File indexing completed on 2023-06-01 00:41:32
0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAssociator.h"
0002 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h"
0003 #include "L1Trigger/DTTriggerPhase2/interface/constants.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 using namespace edm;
0007 using namespace std;
0008 using namespace cmsdt;
0009
0010
0011
0012
0013 MuonPathAssociator::MuonPathAssociator(const ParameterSet &pset,
0014 edm::ConsumesCollector &iC,
0015 std::shared_ptr<GlobalCoordsObtainer> &globalcoordsobtainer)
0016 : debug_(pset.getUntrackedParameter<bool>("debug")),
0017 clean_chi2_correlation_(pset.getParameter<bool>("clean_chi2_correlation")),
0018 useBX_correlation_(pset.getParameter<bool>("useBX_correlation")),
0019 allow_confirmation_(pset.getParameter<bool>("allow_confirmation")),
0020 dT0_correlate_TP_(pset.getParameter<double>("dT0_correlate_TP")),
0021 dBX_correlate_TP_(pset.getParameter<int>("dBX_correlate_TP")),
0022 dTanPsi_correlate_TP_(pset.getParameter<double>("dTanPsi_correlate_TP")),
0023 minx_match_2digis_(pset.getParameter<double>("minx_match_2digis")),
0024 chi2corTh_(pset.getParameter<double>("chi2corTh")) {
0025
0026
0027 if (debug_)
0028 LogDebug("MuonPathAssociator") << "MuonPathAssociator: constructor";
0029
0030
0031 int rawId;
0032 shift_filename_ = pset.getParameter<edm::FileInPath>("shift_filename");
0033 std::ifstream ifin3(shift_filename_.fullPath());
0034 double shift;
0035 if (ifin3.fail()) {
0036 throw cms::Exception("Missing Input File")
0037 << "MuonPathAnalyzerPerSL::MuonPathAnalyzerPerSL() - Cannot find " << shift_filename_.fullPath();
0038 }
0039 while (ifin3.good()) {
0040 ifin3 >> rawId >> shift;
0041 shiftinfo_[rawId] = shift;
0042 }
0043
0044 dtGeomH_ = iC.esConsumes<DTGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0045 globalcoordsobtainer_ = globalcoordsobtainer;
0046 }
0047
0048 MuonPathAssociator::~MuonPathAssociator() {
0049 if (debug_)
0050 LogDebug("MuonPathAssociator") << "MuonPathAssociator: destructor";
0051 }
0052
0053
0054
0055
0056 void MuonPathAssociator::initialise(const edm::EventSetup &iEventSetup) {
0057 if (debug_)
0058 LogDebug("MuonPathAssociator") << "MuonPathAssociator::initialiase";
0059
0060 auto geom = iEventSetup.getHandle(dtGeomH_);
0061 dtGeo_ = &(*geom);
0062 }
0063
0064 void MuonPathAssociator::run(edm::Event &iEvent,
0065 const edm::EventSetup &iEventSetup,
0066 edm::Handle<DTDigiCollection> digis,
0067 std::vector<metaPrimitive> &inMPaths,
0068 std::vector<metaPrimitive> &outMPaths) {
0069 if (dT0_correlate_TP_)
0070 correlateMPaths(digis, inMPaths, outMPaths);
0071 else {
0072 outMPaths.insert(outMPaths.end(), inMPaths.begin(), inMPaths.end());
0073 }
0074 }
0075
0076 void MuonPathAssociator::finish() {
0077 if (debug_)
0078 LogDebug("MuonPathAssociator") << "MuonPathAssociator: finish";
0079 };
0080
0081 void MuonPathAssociator::correlateMPaths(edm::Handle<DTDigiCollection> dtdigis,
0082 std::vector<metaPrimitive> &inMPaths,
0083 std::vector<metaPrimitive> &outMPaths) {
0084 if (debug_)
0085 LogDebug("MuonPathAssociator") << "starting correlation";
0086
0087 for (int wh = -2; wh <= 2; wh++) {
0088 for (int st = 1; st <= 4; st++) {
0089 for (int se = 1; se <= 14; se++) {
0090 if (se >= 13 && st != 4)
0091 continue;
0092
0093 DTChamberId ChId(wh, st, se);
0094 DTSuperLayerId sl1Id(wh, st, se, 1);
0095 DTSuperLayerId sl3Id(wh, st, se, 3);
0096
0097
0098 std::vector<metaPrimitive> SL1metaPrimitives;
0099 for (const auto &metaprimitiveIt : inMPaths) {
0100 if (metaprimitiveIt.rawId == sl1Id.rawId())
0101 SL1metaPrimitives.push_back(metaprimitiveIt);
0102 }
0103
0104
0105 std::vector<metaPrimitive> SL3metaPrimitives;
0106 for (const auto &metaprimitiveIt : inMPaths) {
0107 if (metaprimitiveIt.rawId == sl3Id.rawId())
0108 SL3metaPrimitives.push_back(metaprimitiveIt);
0109 }
0110
0111 if (SL1metaPrimitives.empty() or SL3metaPrimitives.empty())
0112 continue;
0113
0114 if (debug_)
0115 LogDebug("MuonPathAssociator") << "correlating " << SL1metaPrimitives.size() << " metaPrim in SL1 and "
0116 << SL3metaPrimitives.size() << " in SL3 for " << sl3Id;
0117
0118 bool at_least_one_correlation = false;
0119 bool at_least_one_SL1_confirmation = false;
0120 bool at_least_one_SL3_confirmation = false;
0121
0122 vector<bool> useFitSL1;
0123 for (unsigned int i = 0; i < SL1metaPrimitives.size(); i++)
0124 useFitSL1.push_back(false);
0125 vector<bool> useFitSL3;
0126 for (unsigned int i = 0; i < SL3metaPrimitives.size(); i++)
0127 useFitSL3.push_back(false);
0128
0129
0130 vector<metaPrimitive> chamberMetaPrimitives;
0131 vector<metaPrimitive> confirmedMetaPrimitives;
0132 vector<metaPrimitive> normalMetaPrimitives;
0133 int sl1 = 0;
0134 int sl3 = 0;
0135 for (auto SL1metaPrimitive = SL1metaPrimitives.begin(); SL1metaPrimitive != SL1metaPrimitives.end();
0136 ++SL1metaPrimitive, sl1++, sl3 = 0) {
0137 if (clean_chi2_correlation_)
0138 at_least_one_correlation = false;
0139 for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end();
0140 ++SL3metaPrimitive, sl3++) {
0141 if (std::abs(SL1metaPrimitive->tanPhi - SL3metaPrimitive->tanPhi) > dTanPsi_correlate_TP_)
0142 continue;
0143 if (useBX_correlation_) {
0144 if (abs(round(SL1metaPrimitive->t0 / (float)LHC_CLK_FREQ) -
0145 round(SL3metaPrimitive->t0 / (float)LHC_CLK_FREQ)) > dBX_correlate_TP_)
0146 continue;
0147 } else {
0148 if (std::abs(SL1metaPrimitive->t0 - SL3metaPrimitive->t0) >= dT0_correlate_TP_)
0149 continue;
0150 }
0151 long int PosSL1 = (int)round(INCREASED_RES_POS_POW * 10 * SL1metaPrimitive->x);
0152 long int PosSL3 = (int)round(INCREASED_RES_POS_POW * 10 * SL3metaPrimitive->x);
0153 double NewSlope = -999.;
0154
0155 long int pos = (PosSL3 + PosSL1) / 2;
0156
0157
0158 if (((PosSL3 + PosSL1) % 2 != 0) && (pos < 0)) {
0159 pos--;
0160 }
0161
0162 long int difPos_mm_x4 = PosSL3 - PosSL1;
0163 long int tanPsi_x4096_x128 = (difPos_mm_x4)*VERT_PHI1_PHI3_INV;
0164 long int tanpsi = tanPsi_x4096_x128 / ((long int)pow(2, 5 + INCREASED_RES_POS));
0165 if (tanpsi < 0 && tanPsi_x4096_x128 % ((long int)pow(2, 5 + INCREASED_RES_POS)) != 0)
0166 tanpsi--;
0167 NewSlope = -tanpsi / (double)INCREASED_RES_SLOPE_POW;
0168 double MeanT0 = (SL1metaPrimitive->t0 + SL3metaPrimitive->t0) / 2;
0169 double MeanPos = (PosSL3 + PosSL1) / (2. * INCREASED_RES_POS_POW * 10);
0170
0171 DTSuperLayerId SLId1(SL1metaPrimitive->rawId);
0172 DTSuperLayerId SLId3(SL3metaPrimitive->rawId);
0173 DTWireId wireId1(SLId1, 2, 1);
0174 DTWireId wireId3(SLId3, 2, 1);
0175
0176 int shift_sl1 = int(round(shiftinfo_[wireId1.rawId()] * INCREASED_RES_POS_POW * 10));
0177 int shift_sl3 = int(round(shiftinfo_[wireId3.rawId()] * INCREASED_RES_POS_POW * 10));
0178 if (shift_sl1 < shift_sl3) {
0179 pos -= shift_sl1;
0180 } else
0181 pos -= shift_sl3;
0182
0183 int wi[8], tdc[8], lat[8];
0184 wi[0] = SL1metaPrimitive->wi1;
0185 tdc[0] = SL1metaPrimitive->tdc1;
0186 lat[0] = SL1metaPrimitive->lat1;
0187 wi[1] = SL1metaPrimitive->wi2;
0188 tdc[1] = SL1metaPrimitive->tdc2;
0189 lat[1] = SL1metaPrimitive->lat2;
0190 wi[2] = SL1metaPrimitive->wi3;
0191 tdc[2] = SL1metaPrimitive->tdc3;
0192 lat[2] = SL1metaPrimitive->lat3;
0193 wi[3] = SL1metaPrimitive->wi4;
0194 tdc[3] = SL1metaPrimitive->tdc4;
0195 lat[3] = SL1metaPrimitive->lat4;
0196 wi[4] = SL3metaPrimitive->wi1;
0197 tdc[4] = SL3metaPrimitive->tdc1;
0198 lat[4] = SL3metaPrimitive->lat1;
0199 wi[5] = SL3metaPrimitive->wi2;
0200 tdc[5] = SL3metaPrimitive->tdc2;
0201 lat[5] = SL3metaPrimitive->lat2;
0202 wi[6] = SL3metaPrimitive->wi3;
0203 tdc[6] = SL3metaPrimitive->tdc3;
0204 lat[6] = SL3metaPrimitive->lat3;
0205 wi[7] = SL3metaPrimitive->wi4;
0206 tdc[7] = SL3metaPrimitive->tdc4;
0207 lat[7] = SL3metaPrimitive->lat4;
0208
0209 long int chi2 = 0;
0210
0211 long int Z_FACTOR_CORR[8] = {-6, -2, 2, 6, -6, -2, 2, 6};
0212
0213 for (int i = 0; i < 8; i++) {
0214 int sign = 2 * (i / 4) - 1;
0215 Z_FACTOR_CORR[i] = Z_FACTOR_CORR[i] * CELL_HEIGHT + CH_CENTER_TO_MID_SL_X2 * sign;
0216 }
0217 long int sum_A, sum_B;
0218 for (int i = 0; i < NUM_LAYERS_2SL; i++) {
0219 long int shift, slTime;
0220 if (i / NUM_LAYERS == 0) {
0221 shift = shift_sl1;
0222 slTime = SL1metaPrimitive->t0;
0223 } else {
0224 shift = shift_sl3;
0225 slTime = SL3metaPrimitive->t0;
0226 }
0227 if (wi[i] != -1) {
0228 long int drift_dist_um_x4 = DRIFT_SPEED_X4 * (((long int)tdc[i]) - slTime);
0229 long int wireHorizPos_x4 =
0230 (CELL_LENGTH * wi[i] + ((i + 1) % 2) * CELL_SEMILENGTH) * INCREASED_RES_POS_POW;
0231 long int pos_mm_x4;
0232
0233 if (lat[i] == 0) {
0234 pos_mm_x4 = wireHorizPos_x4 - (drift_dist_um_x4 >> 10);
0235 } else {
0236 pos_mm_x4 = wireHorizPos_x4 + (drift_dist_um_x4 >> 10);
0237 }
0238 sum_A = shift + pos_mm_x4 - (long int)round(MeanPos * 10 * INCREASED_RES_POS_POW);
0239 sum_A = sum_A << (14 - INCREASED_RES_POS);
0240 sum_B = Z_FACTOR_CORR[i] * (long int)round(-NewSlope * INCREASED_RES_SLOPE_POW);
0241 chi2 += ((sum_A - sum_B) * (sum_A - sum_B)) >> 2;
0242 }
0243 }
0244
0245 double newChi2 = (double)(chi2 >> INCREASED_RES_POS_POW) / (1024. * 100.);
0246
0247 if (newChi2 > chi2corTh_)
0248 continue;
0249
0250
0251 useFitSL1[sl1] = true;
0252 useFitSL3[sl3] = true;
0253
0254 int quality = 0;
0255 if (SL3metaPrimitive->quality == LOWQ and SL1metaPrimitive->quality == LOWQ)
0256 quality = LOWLOWQ;
0257
0258 if ((SL3metaPrimitive->quality == HIGHQ && SL1metaPrimitive->quality == LOWQ) or
0259 (SL1metaPrimitive->quality == HIGHQ && SL3metaPrimitive->quality == LOWQ))
0260 quality = HIGHLOWQ;
0261
0262 if (SL3metaPrimitive->quality == HIGHQ && SL1metaPrimitive->quality == HIGHQ)
0263 quality = HIGHHIGHQ;
0264
0265 double phi = -999.;
0266 double phiB = -999.;
0267 double phi_cmssw = -999.;
0268 double phiB_cmssw = -999.;
0269 double z = 0;
0270 if (ChId.station() >= 3)
0271 z = Z_SHIFT_MB4;
0272 GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(
0273 LocalPoint(MeanPos, 0., z));
0274 int thisec = ChId.sector();
0275 if (se == 13)
0276 thisec = 4;
0277 if (se == 14)
0278 thisec = 10;
0279 phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1);
0280 double psi = atan(NewSlope);
0281 phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw;
0282
0283 auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), 0, pos, tanpsi);
0284 phi = global_coords[0];
0285 phiB = global_coords[1];
0286
0287 if (!clean_chi2_correlation_)
0288 outMPaths.emplace_back(ChId.rawId(),
0289 MeanT0,
0290 MeanPos,
0291 NewSlope,
0292 phi,
0293 phiB,
0294 phi_cmssw,
0295 phiB_cmssw,
0296 newChi2,
0297 quality,
0298 SL1metaPrimitive->wi1,
0299 SL1metaPrimitive->tdc1,
0300 SL1metaPrimitive->lat1,
0301 SL1metaPrimitive->wi2,
0302 SL1metaPrimitive->tdc2,
0303 SL1metaPrimitive->lat2,
0304 SL1metaPrimitive->wi3,
0305 SL1metaPrimitive->tdc3,
0306 SL1metaPrimitive->lat3,
0307 SL1metaPrimitive->wi4,
0308 SL1metaPrimitive->tdc4,
0309 SL1metaPrimitive->lat4,
0310 SL3metaPrimitive->wi1,
0311 SL3metaPrimitive->tdc1,
0312 SL3metaPrimitive->lat1,
0313 SL3metaPrimitive->wi2,
0314 SL3metaPrimitive->tdc2,
0315 SL3metaPrimitive->lat2,
0316 SL3metaPrimitive->wi3,
0317 SL3metaPrimitive->tdc3,
0318 SL3metaPrimitive->lat3,
0319 SL3metaPrimitive->wi4,
0320 SL3metaPrimitive->tdc4,
0321 SL3metaPrimitive->lat4);
0322 else
0323 chamberMetaPrimitives.emplace_back(ChId.rawId(),
0324 MeanT0,
0325 MeanPos,
0326 NewSlope,
0327 phi,
0328 phiB,
0329 phi_cmssw,
0330 phiB_cmssw,
0331 newChi2,
0332 quality,
0333 SL1metaPrimitive->wi1,
0334 SL1metaPrimitive->tdc1,
0335 SL1metaPrimitive->lat1,
0336 SL1metaPrimitive->wi2,
0337 SL1metaPrimitive->tdc2,
0338 SL1metaPrimitive->lat2,
0339 SL1metaPrimitive->wi3,
0340 SL1metaPrimitive->tdc3,
0341 SL1metaPrimitive->lat3,
0342 SL1metaPrimitive->wi4,
0343 SL1metaPrimitive->tdc4,
0344 SL1metaPrimitive->lat4,
0345 SL3metaPrimitive->wi1,
0346 SL3metaPrimitive->tdc1,
0347 SL3metaPrimitive->lat1,
0348 SL3metaPrimitive->wi2,
0349 SL3metaPrimitive->tdc2,
0350 SL3metaPrimitive->lat2,
0351 SL3metaPrimitive->wi3,
0352 SL3metaPrimitive->tdc3,
0353 SL3metaPrimitive->lat3,
0354 SL3metaPrimitive->wi4,
0355 SL3metaPrimitive->tdc4,
0356 SL3metaPrimitive->lat4);
0357
0358 at_least_one_correlation = true;
0359 }
0360
0361 if (at_least_one_correlation == false &&
0362 allow_confirmation_ == true) {
0363 int matched_digis = 0;
0364 double minx = minx_match_2digis_;
0365 double min2x = minx_match_2digis_;
0366 int best_tdc = -1;
0367 int next_tdc = -1;
0368 int best_wire = -1;
0369 int next_wire = -1;
0370 int best_layer = -1;
0371 int next_layer = -1;
0372 int best_lat = -1;
0373 int next_lat = -1;
0374 int lat = -1;
0375 for (const auto &dtLayerId_It : *dtdigis) {
0376 const DTLayerId dtLId = dtLayerId_It.first;
0377
0378 const DTSuperLayerId dtSLId(dtLId.wheel(), dtLId.station(), dtLId.sector(), dtLId.superLayer());
0379 if (dtSLId.rawId() != sl3Id.rawId())
0380 continue;
0381 double l_shift = 0;
0382 if (dtLId.layer() == 4)
0383 l_shift = X_POS_L4;
0384 else if (dtLId.layer() == 3)
0385 l_shift = X_POS_L3;
0386 else if (dtLId.layer() == 2)
0387 l_shift = -1 * X_POS_L3;
0388 else if (dtLId.layer() == 1)
0389 l_shift = -1 * X_POS_L4;
0390 double x_inSL3 = SL1metaPrimitive->x - SL1metaPrimitive->tanPhi * (VERT_PHI1_PHI3 + l_shift);
0391 for (auto digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; ++digiIt) {
0392 DTWireId wireId(dtLId, (*digiIt).wire());
0393 if ((*digiIt).time() < SL1metaPrimitive->t0)
0394 continue;
0395 double x_wire =
0396 shiftinfo_[wireId.rawId()] + ((*digiIt).time() - SL1metaPrimitive->t0) * DRIFT_SPEED / 10.;
0397 double x_wire_left =
0398 shiftinfo_[wireId.rawId()] - ((*digiIt).time() - SL1metaPrimitive->t0) * DRIFT_SPEED / 10.;
0399 lat = 1;
0400 if (std::abs(x_inSL3 - x_wire) > std::abs(x_inSL3 - x_wire_left)) {
0401 x_wire = x_wire_left;
0402 lat = 0;
0403 }
0404 if (std::abs(x_inSL3 - x_wire) < minx) {
0405
0406
0407
0408
0409 if (dtLId.layer() != best_layer) {
0410 minx = std::abs(x_inSL3 - x_wire);
0411 next_wire = best_wire;
0412 next_tdc = best_tdc;
0413 next_layer = best_layer;
0414 next_lat = best_lat;
0415 matched_digis++;
0416 }
0417 best_wire = (*digiIt).wire();
0418 best_tdc = (*digiIt).time();
0419 best_layer = dtLId.layer();
0420 best_lat = lat;
0421
0422 } else if ((std::abs(x_inSL3 - x_wire) >= minx) && (std::abs(x_inSL3 - x_wire) < min2x)) {
0423
0424 if (dtLId.layer() == best_layer)
0425 continue;
0426
0427
0428
0429
0430 matched_digis++;
0431
0432
0433 min2x = std::abs(x_inSL3 - x_wire);
0434 next_wire = (*digiIt).wire();
0435 next_tdc = (*digiIt).time();
0436 next_layer = dtLId.layer();
0437 next_lat = lat;
0438 }
0439 }
0440 }
0441 if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) {
0442 int new_quality = CHIGHQ;
0443 if (SL1metaPrimitive->quality == LOWQ)
0444 new_quality = CLOWQ;
0445
0446 int wi1 = -1;
0447 int tdc1 = -1;
0448 int lat1 = -1;
0449 int wi2 = -1;
0450 int tdc2 = -1;
0451 int lat2 = -1;
0452 int wi3 = -1;
0453 int tdc3 = -1;
0454 int lat3 = -1;
0455 int wi4 = -1;
0456 int tdc4 = -1;
0457 int lat4 = -1;
0458
0459 if (next_layer == 1) {
0460 wi1 = next_wire;
0461 tdc1 = next_tdc;
0462 lat1 = next_lat;
0463 }
0464 if (next_layer == 2) {
0465 wi2 = next_wire;
0466 tdc2 = next_tdc;
0467 lat2 = next_lat;
0468 }
0469 if (next_layer == 3) {
0470 wi3 = next_wire;
0471 tdc3 = next_tdc;
0472 lat3 = next_lat;
0473 }
0474 if (next_layer == 4) {
0475 wi4 = next_wire;
0476 tdc4 = next_tdc;
0477 lat4 = next_lat;
0478 }
0479
0480 if (best_layer == 1) {
0481 wi1 = best_wire;
0482 tdc1 = best_tdc;
0483 lat1 = best_lat;
0484 }
0485 if (best_layer == 2) {
0486 wi2 = best_wire;
0487 tdc2 = best_tdc;
0488 lat2 = best_lat;
0489 }
0490 if (best_layer == 3) {
0491 wi3 = best_wire;
0492 tdc3 = best_tdc;
0493 lat3 = best_lat;
0494 }
0495 if (best_layer == 4) {
0496 wi4 = best_wire;
0497 tdc4 = best_tdc;
0498 lat4 = best_lat;
0499 }
0500
0501 if (!clean_chi2_correlation_)
0502 outMPaths.emplace_back(metaPrimitive({ChId.rawId(),
0503 SL1metaPrimitive->t0,
0504 SL1metaPrimitive->x,
0505 SL1metaPrimitive->tanPhi,
0506 SL1metaPrimitive->phi,
0507 SL1metaPrimitive->phiB,
0508 SL1metaPrimitive->phi_cmssw,
0509 SL1metaPrimitive->phiB_cmssw,
0510 SL1metaPrimitive->chi2,
0511 new_quality,
0512 SL1metaPrimitive->wi1,
0513 SL1metaPrimitive->tdc1,
0514 SL1metaPrimitive->lat1,
0515 SL1metaPrimitive->wi2,
0516 SL1metaPrimitive->tdc2,
0517 SL1metaPrimitive->lat2,
0518 SL1metaPrimitive->wi3,
0519 SL1metaPrimitive->tdc3,
0520 SL1metaPrimitive->lat3,
0521 SL1metaPrimitive->wi4,
0522 SL1metaPrimitive->tdc4,
0523 SL1metaPrimitive->lat4,
0524 wi1,
0525 tdc1,
0526 lat1,
0527 wi2,
0528 tdc2,
0529 lat2,
0530 wi3,
0531 tdc3,
0532 lat3,
0533 wi4,
0534 tdc4,
0535 lat4,
0536 -1}));
0537 else
0538 confirmedMetaPrimitives.emplace_back(metaPrimitive({ChId.rawId(),
0539 SL1metaPrimitive->t0,
0540 SL1metaPrimitive->x,
0541 SL1metaPrimitive->tanPhi,
0542 SL1metaPrimitive->phi,
0543 SL1metaPrimitive->phiB,
0544 SL1metaPrimitive->phi_cmssw,
0545 SL1metaPrimitive->phiB_cmssw,
0546 SL1metaPrimitive->chi2,
0547 new_quality,
0548 SL1metaPrimitive->wi1,
0549 SL1metaPrimitive->tdc1,
0550 SL1metaPrimitive->lat1,
0551 SL1metaPrimitive->wi2,
0552 SL1metaPrimitive->tdc2,
0553 SL1metaPrimitive->lat2,
0554 SL1metaPrimitive->wi3,
0555 SL1metaPrimitive->tdc3,
0556 SL1metaPrimitive->lat3,
0557 SL1metaPrimitive->wi4,
0558 SL1metaPrimitive->tdc4,
0559 SL1metaPrimitive->lat4,
0560 wi1,
0561 tdc1,
0562 lat1,
0563 wi2,
0564 tdc2,
0565 lat2,
0566 wi3,
0567 tdc3,
0568 lat3,
0569 wi4,
0570 tdc4,
0571 lat4,
0572 -1}));
0573 useFitSL1[sl1] = true;
0574 at_least_one_SL1_confirmation = true;
0575 }
0576 }
0577 }
0578
0579
0580
0581
0582 sl3 = 0;
0583 for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end();
0584 ++SL3metaPrimitive, sl3++) {
0585 if (useFitSL3[sl3])
0586 continue;
0587 if ((at_least_one_correlation == false || clean_chi2_correlation_) &&
0588 allow_confirmation_) {
0589
0590 int matched_digis = 0;
0591 double minx = minx_match_2digis_;
0592 double min2x = minx_match_2digis_;
0593 int best_tdc = -1;
0594 int next_tdc = -1;
0595 int best_wire = -1;
0596 int next_wire = -1;
0597 int best_layer = -1;
0598 int next_layer = -1;
0599 int best_lat = -1;
0600 int next_lat = -1;
0601 int lat = -1;
0602
0603 for (const auto &dtLayerId_It : *dtdigis) {
0604 const DTLayerId dtLId = dtLayerId_It.first;
0605
0606 const DTSuperLayerId dtSLId(dtLId.wheel(), dtLId.station(), dtLId.sector(), dtLId.superLayer());
0607 if (dtSLId.rawId() != sl1Id.rawId())
0608 continue;
0609 double l_shift = 0;
0610 if (dtLId.layer() == 4)
0611 l_shift = X_POS_L4;
0612 if (dtLId.layer() == 3)
0613 l_shift = X_POS_L3;
0614 if (dtLId.layer() == 2)
0615 l_shift = -1 * X_POS_L3;
0616 if (dtLId.layer() == 1)
0617 l_shift = -1 * X_POS_L4;
0618 double x_inSL1 = SL3metaPrimitive->x + SL3metaPrimitive->tanPhi * (VERT_PHI1_PHI3 - l_shift);
0619 for (auto digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; ++digiIt) {
0620 DTWireId wireId(dtLId, (*digiIt).wire());
0621 if ((*digiIt).time() < SL3metaPrimitive->t0)
0622 continue;
0623 double x_wire =
0624 shiftinfo_[wireId.rawId()] + ((*digiIt).time() - SL3metaPrimitive->t0) * DRIFT_SPEED / 10.;
0625 double x_wire_left =
0626 shiftinfo_[wireId.rawId()] - ((*digiIt).time() - SL3metaPrimitive->t0) * DRIFT_SPEED / 10.;
0627 lat = 1;
0628 if (std::abs(x_inSL1 - x_wire) > std::abs(x_inSL1 - x_wire_left)) {
0629 x_wire = x_wire_left;
0630 lat = 0;
0631 }
0632 if (std::abs(x_inSL1 - x_wire) < minx) {
0633
0634
0635
0636
0637 if (dtLId.layer() != best_layer) {
0638 minx = std::abs(x_inSL1 - x_wire);
0639 next_wire = best_wire;
0640 next_tdc = best_tdc;
0641 next_layer = best_layer;
0642 next_lat = best_lat;
0643 matched_digis++;
0644 }
0645 best_wire = (*digiIt).wire();
0646 best_tdc = (*digiIt).time();
0647 best_layer = dtLId.layer();
0648 best_lat = lat;
0649 } else if ((std::abs(x_inSL1 - x_wire) >= minx) && (std::abs(x_inSL1 - x_wire) < min2x)) {
0650
0651 if (dtLId.layer() == best_layer)
0652 continue;
0653
0654
0655
0656
0657 matched_digis++;
0658
0659
0660 min2x = std::abs(x_inSL1 - x_wire);
0661 next_wire = (*digiIt).wire();
0662 next_tdc = (*digiIt).time();
0663 next_layer = dtLId.layer();
0664 next_lat = lat;
0665 }
0666 }
0667 }
0668 if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) {
0669 int new_quality = CHIGHQ;
0670 if (SL3metaPrimitive->quality == LOWQ)
0671 new_quality = CLOWQ;
0672
0673 int wi1 = -1;
0674 int tdc1 = -1;
0675 int lat1 = -1;
0676 int wi2 = -1;
0677 int tdc2 = -1;
0678 int lat2 = -1;
0679 int wi3 = -1;
0680 int tdc3 = -1;
0681 int lat3 = -1;
0682 int wi4 = -1;
0683 int tdc4 = -1;
0684 int lat4 = -1;
0685
0686 if (next_layer == 1) {
0687 wi1 = next_wire;
0688 tdc1 = next_tdc;
0689 lat1 = next_lat;
0690 }
0691 if (next_layer == 2) {
0692 wi2 = next_wire;
0693 tdc2 = next_tdc;
0694 lat2 = next_lat;
0695 }
0696 if (next_layer == 3) {
0697 wi3 = next_wire;
0698 tdc3 = next_tdc;
0699 lat3 = next_lat;
0700 }
0701 if (next_layer == 4) {
0702 wi4 = next_wire;
0703 tdc4 = next_tdc;
0704 lat4 = next_lat;
0705 }
0706
0707 if (best_layer == 1) {
0708 wi1 = best_wire;
0709 tdc1 = best_tdc;
0710 lat1 = best_lat;
0711 }
0712 if (best_layer == 2) {
0713 wi2 = best_wire;
0714 tdc2 = best_tdc;
0715 lat2 = best_lat;
0716 }
0717 if (best_layer == 3) {
0718 wi3 = best_wire;
0719 tdc3 = best_tdc;
0720 lat3 = best_lat;
0721 }
0722 if (best_layer == 4) {
0723 wi4 = best_wire;
0724 tdc4 = best_tdc;
0725 lat4 = best_lat;
0726 }
0727
0728 if (!clean_chi2_correlation_)
0729 outMPaths.push_back(metaPrimitive({ChId.rawId(),
0730 SL3metaPrimitive->t0,
0731 SL3metaPrimitive->x,
0732 SL3metaPrimitive->tanPhi,
0733 SL3metaPrimitive->phi,
0734 SL3metaPrimitive->phiB,
0735 SL3metaPrimitive->phi_cmssw,
0736 SL3metaPrimitive->phiB_cmssw,
0737 SL3metaPrimitive->chi2,
0738 new_quality,
0739 wi1,
0740 tdc1,
0741 lat1,
0742 wi2,
0743 tdc2,
0744 lat2,
0745 wi3,
0746 tdc3,
0747 lat3,
0748 wi4,
0749 tdc4,
0750 lat4,
0751 SL3metaPrimitive->wi1,
0752 SL3metaPrimitive->tdc1,
0753 SL3metaPrimitive->lat1,
0754 SL3metaPrimitive->wi2,
0755 SL3metaPrimitive->tdc2,
0756 SL3metaPrimitive->lat2,
0757 SL3metaPrimitive->wi3,
0758 SL3metaPrimitive->tdc3,
0759 SL3metaPrimitive->lat3,
0760 SL3metaPrimitive->wi4,
0761 SL3metaPrimitive->tdc4,
0762 SL3metaPrimitive->lat4,
0763 -1}));
0764 else
0765 confirmedMetaPrimitives.push_back(metaPrimitive({ChId.rawId(),
0766 SL3metaPrimitive->t0,
0767 SL3metaPrimitive->x,
0768 SL3metaPrimitive->tanPhi,
0769 SL3metaPrimitive->phi,
0770 SL3metaPrimitive->phiB,
0771 SL3metaPrimitive->phi_cmssw,
0772 SL3metaPrimitive->phiB_cmssw,
0773 SL3metaPrimitive->chi2,
0774 new_quality,
0775 wi1,
0776 tdc1,
0777 lat1,
0778 wi2,
0779 tdc2,
0780 lat2,
0781 wi3,
0782 tdc3,
0783 lat3,
0784 wi4,
0785 tdc4,
0786 lat4,
0787 SL3metaPrimitive->wi1,
0788 SL3metaPrimitive->tdc1,
0789 SL3metaPrimitive->lat1,
0790 SL3metaPrimitive->wi2,
0791 SL3metaPrimitive->tdc2,
0792 SL3metaPrimitive->lat2,
0793 SL3metaPrimitive->wi3,
0794 SL3metaPrimitive->tdc3,
0795 SL3metaPrimitive->lat3,
0796 SL3metaPrimitive->wi4,
0797 SL3metaPrimitive->tdc4,
0798 SL3metaPrimitive->lat4,
0799 -1}));
0800 useFitSL3[sl3] = true;
0801 at_least_one_SL3_confirmation = true;
0802 }
0803 }
0804 }
0805
0806 if (clean_chi2_correlation_) {
0807 if (debug_)
0808 LogDebug("MuonPathAssociator") << "Pushing back correlated MPs to the MPs collection";
0809 removeSharingFits(chamberMetaPrimitives, outMPaths);
0810 }
0811 if (clean_chi2_correlation_) {
0812 if (debug_)
0813 LogDebug("MuonPathAssociator") << "Pushing back confirmed MPs to the complete vector";
0814 removeSharingHits(confirmedMetaPrimitives, chamberMetaPrimitives, outMPaths);
0815 }
0816
0817
0818 if (at_least_one_correlation == false || clean_chi2_correlation_) {
0819 if (debug_ && !at_least_one_correlation)
0820 LogDebug("MuonPathAssociator")
0821 << "correlation we found zero correlations, adding both collections as they are to the outMPaths";
0822 if (debug_)
0823 LogDebug("MuonPathAssociator")
0824 << "correlation sizes:" << SL1metaPrimitives.size() << " " << SL3metaPrimitives.size();
0825 if (at_least_one_SL1_confirmation == false || clean_chi2_correlation_) {
0826 sl1 = 0;
0827 for (auto SL1metaPrimitive = SL1metaPrimitives.begin(); SL1metaPrimitive != SL1metaPrimitives.end();
0828 ++SL1metaPrimitive, sl1++) {
0829 if (useFitSL1[sl1])
0830 continue;
0831
0832 DTSuperLayerId SLId(SL1metaPrimitive->rawId);
0833 DTChamberId(SLId.wheel(), SLId.station(), SLId.sector());
0834 metaPrimitive newSL1metaPrimitive = {ChId.rawId(),
0835 SL1metaPrimitive->t0,
0836 SL1metaPrimitive->x,
0837 SL1metaPrimitive->tanPhi,
0838 SL1metaPrimitive->phi,
0839 SL1metaPrimitive->phiB,
0840 SL1metaPrimitive->phi_cmssw,
0841 SL1metaPrimitive->phiB_cmssw,
0842 SL1metaPrimitive->chi2,
0843 SL1metaPrimitive->quality,
0844 SL1metaPrimitive->wi1,
0845 SL1metaPrimitive->tdc1,
0846 SL1metaPrimitive->lat1,
0847 SL1metaPrimitive->wi2,
0848 SL1metaPrimitive->tdc2,
0849 SL1metaPrimitive->lat2,
0850 SL1metaPrimitive->wi3,
0851 SL1metaPrimitive->tdc3,
0852 SL1metaPrimitive->lat3,
0853 SL1metaPrimitive->wi4,
0854 SL1metaPrimitive->tdc4,
0855 SL1metaPrimitive->lat4,
0856 -1,
0857 -1,
0858 -1,
0859 -1,
0860 -1,
0861 -1,
0862 -1,
0863 -1,
0864 -1,
0865 -1,
0866 -1,
0867 -1,
0868 -1};
0869
0870 bool ok = true;
0871 for (auto &metaPrimitive : chamberMetaPrimitives) {
0872 if (!isNotAPrimo(newSL1metaPrimitive, metaPrimitive)) {
0873 ok = false;
0874 break;
0875 }
0876 }
0877 if (!ok)
0878 continue;
0879
0880 if (!clean_chi2_correlation_)
0881 outMPaths.push_back(newSL1metaPrimitive);
0882 else
0883 normalMetaPrimitives.push_back(newSL1metaPrimitive);
0884 }
0885 }
0886 if (at_least_one_SL3_confirmation == false || clean_chi2_correlation_) {
0887 sl3 = 0;
0888 for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end();
0889 ++SL3metaPrimitive, sl3++) {
0890 if (useFitSL3[sl3])
0891 continue;
0892 DTSuperLayerId SLId(SL3metaPrimitive->rawId);
0893 DTChamberId(SLId.wheel(), SLId.station(), SLId.sector());
0894 metaPrimitive newSL3metaPrimitive = {ChId.rawId(),
0895 SL3metaPrimitive->t0,
0896 SL3metaPrimitive->x,
0897 SL3metaPrimitive->tanPhi,
0898 SL3metaPrimitive->phi,
0899 SL3metaPrimitive->phiB,
0900 SL3metaPrimitive->phi_cmssw,
0901 SL3metaPrimitive->phiB_cmssw,
0902 SL3metaPrimitive->chi2,
0903 SL3metaPrimitive->quality,
0904 -1,
0905 -1,
0906 -1,
0907 -1,
0908 -1,
0909 -1,
0910 -1,
0911 -1,
0912 -1,
0913 -1,
0914 -1,
0915 -1,
0916 SL3metaPrimitive->wi1,
0917 SL3metaPrimitive->tdc1,
0918 SL3metaPrimitive->lat1,
0919 SL3metaPrimitive->wi2,
0920 SL3metaPrimitive->tdc2,
0921 SL3metaPrimitive->lat2,
0922 SL3metaPrimitive->wi3,
0923 SL3metaPrimitive->tdc3,
0924 SL3metaPrimitive->lat3,
0925 SL3metaPrimitive->wi4,
0926 SL3metaPrimitive->tdc4,
0927 SL3metaPrimitive->lat4,
0928 -1};
0929
0930 if (!clean_chi2_correlation_)
0931 outMPaths.push_back(newSL3metaPrimitive);
0932 else
0933 normalMetaPrimitives.push_back(newSL3metaPrimitive);
0934 }
0935 }
0936 }
0937
0938 SL1metaPrimitives.clear();
0939 SL1metaPrimitives.erase(SL1metaPrimitives.begin(), SL1metaPrimitives.end());
0940 SL3metaPrimitives.clear();
0941 SL3metaPrimitives.erase(SL3metaPrimitives.begin(), SL3metaPrimitives.end());
0942
0943 vector<metaPrimitive> auxMetaPrimitives;
0944 if (clean_chi2_correlation_) {
0945 if (debug_)
0946 LogDebug("MuonPathAssociator") << "Pushing back normal MPs to the auxiliar vector";
0947 removeSharingHits(normalMetaPrimitives, confirmedMetaPrimitives, auxMetaPrimitives);
0948 }
0949 if (clean_chi2_correlation_) {
0950 if (debug_)
0951 LogDebug("MuonPathAssociator") << "Pushing back normal MPs to the MPs collection";
0952 removeSharingHits(auxMetaPrimitives, chamberMetaPrimitives, outMPaths);
0953 }
0954 }
0955 }
0956 }
0957
0958
0959 std::vector<metaPrimitive> SL2metaPrimitives;
0960
0961 for (int wh = -2; wh <= 2; wh++) {
0962 for (int st = 1; st <= 4; st++) {
0963 for (int se = 1; se <= 14; se++) {
0964 if (se >= 13 && st != 4)
0965 continue;
0966
0967 DTChamberId ChId(wh, st, se);
0968 DTSuperLayerId sl2Id(wh, st, se, 2);
0969
0970
0971 for (auto metaprimitiveIt = inMPaths.begin(); metaprimitiveIt != inMPaths.end(); ++metaprimitiveIt)
0972 if (metaprimitiveIt->rawId == sl2Id.rawId()) {
0973 SL2metaPrimitives.push_back(*metaprimitiveIt);
0974 if (debug_)
0975 printmPC(*metaprimitiveIt);
0976 outMPaths.push_back(*metaprimitiveIt);
0977 }
0978 }
0979 }
0980 }
0981
0982 LogDebug("MuonPathAssociator") << "\t etaTP: added " << SL2metaPrimitives.size() << "to outMPaths" << std::endl;
0983
0984 SL2metaPrimitives.clear();
0985 SL2metaPrimitives.erase(SL2metaPrimitives.begin(), SL2metaPrimitives.end());
0986 }
0987
0988 void MuonPathAssociator::removeSharingFits(vector<metaPrimitive> &chamberMPaths, vector<metaPrimitive> &allMPaths) {
0989 vector<bool> useFit;
0990 for (unsigned int i = 0; i < chamberMPaths.size(); i++) {
0991 useFit.push_back(true);
0992 }
0993 for (unsigned int i = 0; i < chamberMPaths.size(); i++) {
0994 if (debug_)
0995 LogDebug("MuonPathAssociator") << "Looking at prim" << i;
0996 if (!useFit[i])
0997 continue;
0998 for (unsigned int j = i + 1; j < chamberMPaths.size(); j++) {
0999 if (debug_)
1000 LogDebug("MuonPathAssociator") << "Comparing with prim " << j;
1001 if (!useFit[j])
1002 continue;
1003 metaPrimitive first = chamberMPaths[i];
1004 metaPrimitive second = chamberMPaths[j];
1005 if (shareFit(first, second)) {
1006 if (first.quality > second.quality)
1007 useFit[j] = false;
1008 else if (first.quality < second.quality)
1009 useFit[i] = false;
1010 else {
1011 if (first.chi2 < second.chi2)
1012 useFit[j] = false;
1013 else {
1014 useFit[i] = false;
1015 break;
1016 }
1017 }
1018 }
1019 }
1020 if (useFit[i]) {
1021 if (debug_)
1022 printmPC(chamberMPaths[i]);
1023 allMPaths.push_back(chamberMPaths[i]);
1024 }
1025 }
1026 if (debug_)
1027 LogDebug("MuonPathAssociator") << "---Swapping chamber---";
1028 }
1029
1030 void MuonPathAssociator::removeSharingHits(std::vector<metaPrimitive> &firstMPaths,
1031 std::vector<metaPrimitive> &secondMPaths,
1032 std::vector<metaPrimitive> &allMPaths) {
1033 for (auto &firstMP : firstMPaths) {
1034 if (debug_)
1035 LogDebug("MuonPathAssociator") << "----------------------------------";
1036 if (debug_)
1037 LogDebug("MuonPathAssociator") << "Turn for ";
1038 if (debug_)
1039 printmPC(firstMP);
1040 bool ok = true;
1041 for (auto &secondMP : secondMPaths) {
1042 if (debug_)
1043 LogDebug("MuonPathAssociator") << "Comparing with ";
1044 if (debug_)
1045 printmPC(secondMP);
1046 if (!isNotAPrimo(firstMP, secondMP)) {
1047 ok = false;
1048 break;
1049 }
1050 }
1051 if (ok) {
1052 allMPaths.push_back(firstMP);
1053 if (debug_)
1054 printmPC(firstMP);
1055 }
1056 if (debug_)
1057 LogDebug("MuonPathAssociator") << "----------------------------------";
1058 }
1059 }
1060
1061 bool MuonPathAssociator::shareFit(metaPrimitive first, metaPrimitive second) {
1062 bool lay1 = (first.wi1 == second.wi1) && (first.tdc1 = second.tdc1);
1063 bool lay2 = (first.wi2 == second.wi2) && (first.tdc2 = second.tdc2);
1064 bool lay3 = (first.wi3 == second.wi3) && (first.tdc3 = second.tdc3);
1065 bool lay4 = (first.wi4 == second.wi4) && (first.tdc4 = second.tdc4);
1066 bool lay5 = (first.wi5 == second.wi5) && (first.tdc5 = second.tdc5);
1067 bool lay6 = (first.wi6 == second.wi6) && (first.tdc6 = second.tdc6);
1068 bool lay7 = (first.wi7 == second.wi7) && (first.tdc7 = second.tdc7);
1069 bool lay8 = (first.wi8 == second.wi8) && (first.tdc8 = second.tdc8);
1070
1071 if (lay1 && lay2 && lay3 && lay4) {
1072 if (lay5 || lay6 || lay7 || lay8)
1073 return true;
1074 else
1075 return false;
1076 } else if (lay5 && lay6 && lay7 && lay8) {
1077 if (lay1 || lay2 || lay3 || lay4)
1078 return true;
1079 else
1080 return false;
1081 } else
1082 return false;
1083 }
1084
1085 bool MuonPathAssociator::isNotAPrimo(metaPrimitive first, metaPrimitive second) {
1086 int hitsSL1 = (first.wi1 != -1) + (first.wi2 != -1) + (first.wi3 != -1) + (first.wi4 != -1);
1087 int hitsSL3 = (first.wi5 != -1) + (first.wi6 != -1) + (first.wi7 != -1) + (first.wi8 != -1);
1088
1089 bool lay1 = (first.wi1 == second.wi1) && (first.tdc1 = second.tdc1) && (first.wi1 != -1);
1090 bool lay2 = (first.wi2 == second.wi2) && (first.tdc2 = second.tdc2) && (first.wi2 != -1);
1091 bool lay3 = (first.wi3 == second.wi3) && (first.tdc3 = second.tdc3) && (first.wi3 != -1);
1092 bool lay4 = (first.wi4 == second.wi4) && (first.tdc4 = second.tdc4) && (first.wi4 != -1);
1093 bool lay5 = (first.wi5 == second.wi5) && (first.tdc5 = second.tdc5) && (first.wi5 != -1);
1094 bool lay6 = (first.wi6 == second.wi6) && (first.tdc6 = second.tdc6) && (first.wi6 != -1);
1095 bool lay7 = (first.wi7 == second.wi7) && (first.tdc7 = second.tdc7) && (first.wi7 != -1);
1096 bool lay8 = (first.wi8 == second.wi8) && (first.tdc8 = second.tdc8) && (first.wi8 != -1);
1097
1098 return (((!lay1 && !lay2 && !lay3 && !lay4) || hitsSL1 < 3) && ((!lay5 && !lay6 && !lay7 && !lay8) || hitsSL3 < 3));
1099 }
1100
1101 void MuonPathAssociator::printmPC(metaPrimitive mP) {
1102 DTChamberId ChId(mP.rawId);
1103 LogDebug("MuonPathAssociator") << ChId << "\t"
1104 << " " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 << " "
1105 << setw(2) << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(2)
1106 << left << mP.wi5 << " " << setw(2) << left << mP.wi6 << " " << setw(2) << left
1107 << mP.wi7 << " " << setw(2) << left << mP.wi8 << " " << setw(5) << left << mP.tdc1
1108 << " " << setw(5) << left << mP.tdc2 << " " << setw(5) << left << mP.tdc3 << " "
1109 << setw(5) << left << mP.tdc4 << " " << setw(5) << left << mP.tdc5 << " " << setw(5)
1110 << left << mP.tdc6 << " " << setw(5) << left << mP.tdc7 << " " << setw(5) << left
1111 << mP.tdc8 << " " << setw(2) << left << mP.lat1 << " " << setw(2) << left << mP.lat2
1112 << " " << setw(2) << left << mP.lat3 << " " << setw(2) << left << mP.lat4 << " "
1113 << setw(2) << left << mP.lat5 << " " << setw(2) << left << mP.lat6 << " " << setw(2)
1114 << left << mP.lat7 << " " << setw(2) << left << mP.lat8 << " " << setw(10) << right
1115 << mP.x << " " << setw(9) << left << mP.tanPhi << " " << setw(5) << left << mP.t0
1116 << " " << setw(13) << left << mP.chi2 << " \n";
1117 }