13 debug_(
pset.getUntrackedParameter<
bool>(
"debug")),
14 chi2Th_(
pset.getUntrackedParameter<double>(
"chi2Th")),
18 chiSquareThreshold_(50),
19 minHits4Fit_(
pset.getUntrackedParameter<
int>(
"minHits4Fit")) {
23 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MuonPathAnalyzer: constructor";
35 while (ifin3.good()) {
36 ifin3 >> rawId >>
shift;
45 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MuonPathAnalyzer: destructor";
53 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MuonPathAnalyzerInChamber::initialiase";
64 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MuonPathAnalyzerInChamber: run";
67 for (
auto muonpath = muonpaths.begin(); muonpath != muonpaths.end(); ++muonpath) {
68 analyze(*muonpath, outmuonpaths);
74 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MuonPathAnalyzer: finish";
82 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2:analyze \t\t\t\t starts";
86 LogDebug(
"MuonPathAnalyzerInChamber") << inMPath->nprimitives();
87 auto mPath = std::make_shared<MuonPath>(inMPath);
90 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2::analyze, looking at mPath: ";
91 for (
int i = 0;
i < mPath->nprimitives();
i++)
92 LogDebug(
"MuonPathAnalyzerInChamber")
93 << mPath->primitive(
i)->layerId() <<
" , " << mPath->primitive(
i)->superLayerId() <<
" , "
94 << mPath->primitive(
i)->channelId() <<
" , " << mPath->primitive(
i)->laterality();
98 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2:analyze \t\t\t\t\t is Analyzable? ";
99 if (!mPath->isAnalyzable())
102 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2:analyze \t\t\t\t\t yes it is analyzable " << mPath->isAnalyzable();
108 std::shared_ptr<MuonPath> mpAux;
110 float best_chi2 = 99999.;
113 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2:analyze \t\t\t\t\t Start with combination " <<
i;
117 for (
int ii = 0;
ii < 8;
ii++) {
118 xwire[
ii] = mPath->xWirePos(
ii);
119 if (xwire[
ii] == 0) {
120 present_layer[
ii] = 0;
123 present_layer[
ii] = 1;
128 mPath->setChiSquare(0);
130 if (mPath->chiSquare() != 0)
143 double jm_x = (mPath->horizPos());
145 for (
int i = 0;
i < mPath->nprimitives();
i++) {
146 if (mPath->primitive(
i)->isValidTime()) {
147 selected_Id = mPath->primitive(
i)->cameraId();
148 mPath->setRawId(selected_Id);
159 int thisec = MuonPathSLId.sector();
164 double phi = jm_x_cmssw_global.phi() -
PHI_CONV * (thisec - 1);
165 double psi = atan(mPath->tanPhi());
166 mPath->setPhi(jm_x_cmssw_global.phi() -
PHI_CONV * (thisec - 1));
169 if (mPath->chiSquare() < best_chi2 && mPath->chiSquare() > 0) {
170 mpAux = std::make_shared<MuonPath>(mPath);
172 best_chi2 = mPath->chiSquare();
175 if (mpAux !=
nullptr) {
178 LogDebug(
"MuonPathAnalyzerInChamber")
179 <<
"DTp2:analize \t\t\t\t\t Laterality " << bestI <<
" is the one with smaller chi2";
182 LogDebug(
"MuonPathAnalyzerInChamber")
183 <<
"DTp2:analize \t\t\t\t\t No Laterality found with chi2 smaller than threshold";
186 LogDebug(
"MuonPathAnalyzerInChamber") <<
"DTp2:analize \t\t\t\t\t Ended working with this set of lateralities";
190 for (
int i = 0;
i <= mpath->nprimitives();
i++) {
191 if (mpath->primitive(
i)->isValidTime())
199 for (
int i = 0;
i <= mpath->nprimitives();
i++) {
212 LogDebug(
"MuonPathAnalyzerInChamber") <<
"MPAnalyzer::buildLateralities << setLateralitiesFromPrims ";
213 mpath->setLateralCombFromPrimitives();
222 for (
int ilat = 0; ilat < NLayers; ilat++) {
226 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[DEBUG_] Input[" << ilat <<
"]: " << lr;
231 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[DEBUG_] - Adding it to " <<
lateralities_.size() <<
" lists...";
232 for (
unsigned int iall = 0; iall <
lateralities_.size(); iall++) {
243 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[DEBUG_] - Duplicating " << ncurrentoptions <<
" lists...";
246 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[DEBUG_] - Now we have " <<
lateralities_.size() <<
" lists...";
249 for (
unsigned int iall = 0; iall < ncurrentoptions; iall++) {
259 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[WARNING]: TOO MANY LATERALITIES TO CHECK !!";
260 LogDebug(
"MuonPathAnalyzerInChamber") <<
"[WARNING]: skipping this muon";
268 for (
unsigned int iall = 0; iall <
lateralities_.size(); iall++) {
269 LogDebug(
"MuonPathAnalyzerInChamber") << iall <<
" -> [";
270 for (
int ilat = 0; ilat < NLayers; ilat++) {
272 LogDebug(
"MuonPathAnalyzerInChamber") <<
",";
275 LogDebug(
"MuonPathAnalyzerInChamber") <<
"]";
281 for (
int i = 0;
i < 8;
i++)
284 mpath->setLateralComb(
tmp);
288 for (
int i = 0;
i < mpath->nprimitives();
i++) {
289 if (mpath->primitive(
i)->isValidTime()) {
290 selected_Id = mpath->primitive(
i)->cameraId();
291 mpath->setRawId(selected_Id);
298 LogDebug(
"MuonPathAnalyzerInChamber")
299 <<
"Id " << chId.rawId() <<
" Wh " << chId.wheel() <<
" St " << chId.station() <<
" Se " << chId.sector();
300 mpath->setRawId(chId.rawId());
304 DTWireId wireId1(MuonPathSLId1, 2, 1);
305 DTWireId wireId3(MuonPathSLId3, 2, 1);
308 LogDebug(
"MuonPathAnalyzerInChamber")
312 float zwire[
NUM_LAYERS_2SL] = {-13.7, -12.4, -11.1, -9.8002, 9.79999, 11.1, 12.4, 13.7};
313 for (
int i = 0;
i < mpath->nprimitives();
i++) {
314 if (mpath->primitive(
i)->isValidTime()) {
316 mpath->setXWirePos(10000 *
shiftinfo_[wireId1.rawId()] +
317 (mpath->primitive(
i)->channelId() + 0.5 * (double)((
i + 1) % 2)) *
delta,
321 (mpath->primitive(
i)->channelId() + 0.5 * (double)((
i + 1) % 2)) *
delta,
323 mpath->setZWirePos(zwire[
i] * 1000,
i);
324 mpath->setTWireTDC(mpath->primitive(
i)->tdcTimeStamp() *
DRIFT_SPEED,
i);
326 mpath->setXWirePos(0.,
i);
327 mpath->setZWirePos(0.,
i);
331 LogDebug(
"MuonPathAnalyzerInChamber") << mpath->primitive(
i)->tdcTimeStamp() <<
" ";
334 LogDebug(
"MuonPathAnalyzerInChamber");
337 TLateralities laterality,
342 for (
int i = 0;
i < 8;
i++) {
343 xwire[
i] = mpath->xWirePos(
i);
344 zwire[
i] = mpath->zWirePos(
i);
345 tTDCvdrift[
i] = mpath->tWireTDC(
i);
353 for (
int lay = 0; lay < 8; lay++) {
355 LogDebug(
"MuonPathAnalyzerInChamber") <<
"In fitPerLat " << lay <<
" xwire " << xwire[lay] <<
" zwire "
356 << zwire[lay] <<
" tTDCvdrift " << tTDCvdrift[lay];
357 xhit[lay] = xwire[lay] + (-1 + 2 * laterality[lay]) * 1000 * tTDCvdrift[lay];
359 LogDebug(
"MuonPathAnalyzerInChamber") <<
"In fitPerLat " << lay <<
" xhit " << xhit[lay];
370 for (
int lay = 0; lay < 8; lay++) {
371 if (present_layer[lay] == 0)
374 LogDebug(
"MuonPathAnalyzerInChamber")
375 <<
" For layer " << lay + 1 <<
" xwire[lay] " << xwire[lay] <<
" zwire " << zwire[lay] <<
" b " <<
b[lay];
377 LogDebug(
"MuonPathAnalyzerInChamber") <<
" xhit[lat][lay] " << xhit[lay];
378 cbscal = (-1 + 2 * laterality[lay]) *
b[lay] + cbscal;
379 zbscal = zwire[lay] *
b[lay] + zbscal;
380 czscal = (-1 + 2 * laterality[lay]) * zwire[lay] + czscal;
382 bbscal =
b[lay] *
b[lay] + bbscal;
383 zzscal = zwire[lay] * zwire[lay] + zzscal;
384 ccscal = (-1 + 2 * laterality[lay]) * (-1 + 2 * laterality[lay]) + ccscal;
394 cz = (cbscal * zbscal - czscal * bbscal) / (zzscal * bbscal - zbscal * zbscal);
395 cb = (czscal * zbscal - cbscal * zzscal) / (zzscal * bbscal - zbscal * zbscal);
397 zb = (czscal * cbscal - zbscal * ccscal) / (bbscal * ccscal - cbscal * cbscal);
398 zc = (zbscal * cbscal - czscal * bbscal) / (bbscal * ccscal - cbscal * cbscal);
400 bc = (zbscal * czscal - cbscal * zzscal) / (ccscal * zzscal - czscal * czscal);
401 bz = (cbscal * czscal - zbscal * ccscal) / (ccscal * zzscal - czscal * czscal);
407 for (
int lay = 0; lay < 8; lay++) {
408 if (present_layer[lay] == 0)
411 LogDebug(
"MuonPathAnalyzerInChamber")
412 <<
" For layer " << lay + 1 <<
" xwire[lay] " << xwire[lay] <<
" zwire " << zwire[lay] <<
" b " <<
b[lay];
413 c_tilde[lay] = (-1 + 2 * laterality[lay]) + cz * zwire[lay] + cb *
b[lay];
414 z_tilde[lay] = zwire[lay] + zb *
b[lay] + zc * (-1 + 2 * laterality[lay]);
415 b_tilde[lay] =
b[lay] + bc * (-1 + 2 * laterality[lay]) + bz * zwire[lay];
419 double xctilde = 0.0;
420 double xztilde = 0.0;
421 double xbtilde = 0.0;
422 double ctildectilde = 0.0;
423 double ztildeztilde = 0.0;
424 double btildebtilde = 0.0;
426 double rect0vdrift = 0.0;
427 double recslope = 0.0;
430 for (
int lay = 0; lay < 8; lay++) {
431 if (present_layer[lay] == 0)
433 xctilde = xhit[lay] * c_tilde[lay] + xctilde;
434 ctildectilde = c_tilde[lay] * c_tilde[lay] + ctildectilde;
435 xztilde = xhit[lay] * z_tilde[lay] + xztilde;
436 ztildeztilde = z_tilde[lay] * z_tilde[lay] + ztildeztilde;
437 xbtilde = xhit[lay] * b_tilde[lay] + xbtilde;
438 btildebtilde = b_tilde[lay] * b_tilde[lay] + btildebtilde;
442 rect0vdrift = xctilde / ctildectilde;
443 recslope = xztilde / ztildeztilde;
444 recpos = xbtilde / btildebtilde;
446 LogDebug(
"MuonPathAnalyzerInChamber") <<
" In fitPerLat Reconstructed values per lat "
447 <<
" rect0vdrift " << rect0vdrift;
448 LogDebug(
"MuonPathAnalyzerInChamber")
450 << recslope <<
" recpos " << recpos;
456 double recchi2 = 0.0;
457 int sign_tdriftvdrift = {0};
458 int incell_tdriftvdrift = {0};
459 int physical_slope = {0};
465 for (
int lay = 0; lay < 8; lay++) {
466 if (present_layer[lay] == 0)
468 rectdriftvdrift[lay] = tTDCvdrift[lay] - rect0vdrift / 1000;
470 LogDebug(
"MuonPathAnalyzerInChamber") << rectdriftvdrift[lay];
471 recres[lay] = xhit[lay] - zwire[lay] * recslope -
b[lay] * recpos - (-1 + 2 * laterality[lay]) * rect0vdrift;
472 if ((present_layer[lay] == 1) && (rectdriftvdrift[lay] < -0.1)) {
473 sign_tdriftvdrift = -1;
474 if (-0.1 - rectdriftvdrift[lay] > maxDif) {
475 maxDif = -0.1 - rectdriftvdrift[lay];
479 if ((present_layer[lay] == 1) && (
abs(rectdriftvdrift[lay]) > 21.1)) {
480 incell_tdriftvdrift = -1;
481 if (rectdriftvdrift[lay] - 21.1 > maxDif) {
482 maxDif = rectdriftvdrift[lay] - 21.1;
488 if (fabs(recslope / 10) > 1.3)
491 if (physical_slope == -1 &&
debug_)
492 LogDebug(
"MuonPathAnalyzerInChamber") <<
"Combination with UNPHYSICAL slope ";
493 if (sign_tdriftvdrift == -1 &&
debug_)
494 LogDebug(
"MuonPathAnalyzerInChamber") <<
"Combination with negative tdrift-vdrift ";
495 if (incell_tdriftvdrift == -1 &&
debug_)
496 LogDebug(
"MuonPathAnalyzerInChamber") <<
"Combination with tdrift-vdrift larger than half cell ";
498 for (
int lay = 0; lay < 8; lay++) {
499 if (present_layer[lay] == 0)
501 recchi2 = recres[lay] * recres[lay] + recchi2;
504 LogDebug(
"MuonPathAnalyzerInChamber")
505 <<
"In fitPerLat Chi2 " << recchi2 <<
" with sign " << sign_tdriftvdrift <<
" within cell "
506 << incell_tdriftvdrift <<
" physical_slope " << physical_slope;
509 if (
true && maxInt != -1) {
510 present_layer[maxInt] = 0;
512 LogDebug(
"MuonPathAnalyzerInChamber") <<
"We get rid of hit in layer " << maxInt;
516 if (!(sign_tdriftvdrift == -1) && !(incell_tdriftvdrift == -1) && !(physical_slope == -1)) {
517 mpath->setBxTimeValue((rect0vdrift /
DRIFT_SPEED) / 1000);
518 mpath->setTanPhi(-1 * recslope / 10);
519 mpath->setHorizPos(recpos / 10000);
520 mpath->setChiSquare(recchi2 / 100000000);
523 LogDebug(
"MuonPathAnalyzerInChamber")
525 <<
"t0 " << mpath->bxTimeValue() <<
" slope " << mpath->tanPhi() <<
" pos " << mpath->horizPos() <<
" chi2 "
526 << mpath->chiSquare() <<
" rawId " << mpath->rawId();
530 mPath->setQuality(
NOPATH);
532 if (mPath->nprimitivesUp() >= 4 && mPath->nprimitivesDown() >= 4) {
534 }
else if ((mPath->nprimitivesUp() == 4 && mPath->nprimitivesDown() == 3) ||
535 (mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() == 4)) {
537 }
else if ((mPath->nprimitivesUp() == 4 && mPath->nprimitivesDown() <= 2 && mPath->nprimitivesDown() > 0) ||
538 (mPath->nprimitivesUp() <= 2 && mPath->nprimitivesUp() > 0 && mPath->nprimitivesDown() == 4)) {
539 mPath->setQuality(
CHIGHQ);
540 }
else if ((mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() == 3)) {
542 }
else if ((mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() <= 2 && mPath->nprimitivesDown() > 0) ||
543 (mPath->nprimitivesUp() <= 2 && mPath->nprimitivesUp() > 0 && mPath->nprimitivesDown() == 3) ||
544 (mPath->nprimitivesUp() == 2 && mPath->nprimitivesDown() == 2)) {
545 mPath->setQuality(
CLOWQ);
546 }
else if (mPath->nprimitivesUp() >= 4 || mPath->nprimitivesDown() >= 4) {
547 mPath->setQuality(
HIGHQ);
548 }
else if (mPath->nprimitivesUp() == 3 || mPath->nprimitivesDown() == 3) {
549 mPath->setQuality(
LOWQ);