File indexing completed on 2022-05-12 01:51:37
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() and 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 bool useFitSL1[SL1metaPrimitives.size()];
0123 for (unsigned int i = 0; i < SL1metaPrimitives.size(); i++)
0124 useFitSL1[i] = false;
0125 bool useFitSL3[SL3metaPrimitives.size()];
0126 for (unsigned int i = 0; i < SL3metaPrimitives.size(); i++)
0127 useFitSL3[i] = 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 = -1) {
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(), 3);
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 double x_wire =
0394 shiftinfo_[wireId.rawId()] + ((*digiIt).time() - SL1metaPrimitive->t0) * DRIFT_SPEED / 10.;
0395 double x_wire_left =
0396 shiftinfo_[wireId.rawId()] - ((*digiIt).time() - SL1metaPrimitive->t0) * DRIFT_SPEED / 10.;
0397 lat = 1;
0398 if (std::abs(x_inSL3 - x_wire) > std::abs(x_inSL3 - x_wire_left)) {
0399 x_wire = x_wire_left;
0400 lat = 0;
0401 }
0402 if (std::abs(x_inSL3 - x_wire) < minx) {
0403 min2x = minx;
0404 minx = std::abs(x_inSL3 - x_wire);
0405 next_wire = best_wire;
0406 next_tdc = best_tdc;
0407 next_layer = best_layer;
0408 next_lat = best_lat;
0409
0410 best_wire = (*digiIt).wire();
0411 best_tdc = (*digiIt).time();
0412 best_layer = dtLId.layer();
0413 best_lat = lat;
0414 matched_digis++;
0415 } else if ((std::abs(x_inSL3 - x_wire) >= minx) && (std::abs(x_inSL3 - x_wire) < min2x)) {
0416 min2x = std::abs(x_inSL3 - x_wire);
0417 next_wire = (*digiIt).wire();
0418 next_tdc = (*digiIt).time();
0419 next_layer = dtLId.layer();
0420 next_lat = lat;
0421 matched_digis++;
0422 }
0423 }
0424 }
0425 if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) {
0426 int new_quality = CHIGHQ;
0427 if (SL1metaPrimitive->quality == LOWQ)
0428 new_quality = CLOWQ;
0429
0430 int wi1 = -1;
0431 int tdc1 = -1;
0432 int lat1 = -1;
0433 int wi2 = -1;
0434 int tdc2 = -1;
0435 int lat2 = -1;
0436 int wi3 = -1;
0437 int tdc3 = -1;
0438 int lat3 = -1;
0439 int wi4 = -1;
0440 int tdc4 = -1;
0441 int lat4 = -1;
0442
0443 if (next_layer == 1) {
0444 wi1 = next_wire;
0445 tdc1 = next_tdc;
0446 lat1 = next_lat;
0447 }
0448 if (next_layer == 2) {
0449 wi2 = next_wire;
0450 tdc2 = next_tdc;
0451 lat2 = next_lat;
0452 }
0453 if (next_layer == 3) {
0454 wi3 = next_wire;
0455 tdc3 = next_tdc;
0456 lat3 = next_lat;
0457 }
0458 if (next_layer == 4) {
0459 wi4 = next_wire;
0460 tdc4 = next_tdc;
0461 lat4 = next_lat;
0462 }
0463
0464 if (best_layer == 1) {
0465 wi1 = best_wire;
0466 tdc1 = best_tdc;
0467 lat1 = best_lat;
0468 }
0469 if (best_layer == 2) {
0470 wi2 = best_wire;
0471 tdc2 = best_tdc;
0472 lat2 = best_lat;
0473 }
0474 if (best_layer == 3) {
0475 wi3 = best_wire;
0476 tdc3 = best_tdc;
0477 lat3 = best_lat;
0478 }
0479 if (best_layer == 4) {
0480 wi4 = best_wire;
0481 tdc4 = best_tdc;
0482 lat4 = best_lat;
0483 }
0484
0485 if (!clean_chi2_correlation_)
0486 outMPaths.emplace_back(metaPrimitive({ChId.rawId(),
0487 SL1metaPrimitive->t0,
0488 SL1metaPrimitive->x,
0489 SL1metaPrimitive->tanPhi,
0490 SL1metaPrimitive->phi,
0491 SL1metaPrimitive->phiB,
0492 SL1metaPrimitive->phi_cmssw,
0493 SL1metaPrimitive->phiB_cmssw,
0494 SL1metaPrimitive->chi2,
0495 new_quality,
0496 SL1metaPrimitive->wi1,
0497 SL1metaPrimitive->tdc1,
0498 SL1metaPrimitive->lat1,
0499 SL1metaPrimitive->wi2,
0500 SL1metaPrimitive->tdc2,
0501 SL1metaPrimitive->lat2,
0502 SL1metaPrimitive->wi3,
0503 SL1metaPrimitive->tdc3,
0504 SL1metaPrimitive->lat3,
0505 SL1metaPrimitive->wi4,
0506 SL1metaPrimitive->tdc4,
0507 SL1metaPrimitive->lat4,
0508 wi1,
0509 tdc1,
0510 lat1,
0511 wi2,
0512 tdc2,
0513 lat2,
0514 wi3,
0515 tdc3,
0516 lat3,
0517 wi4,
0518 tdc4,
0519 lat4,
0520 -1}));
0521 else
0522 confirmedMetaPrimitives.emplace_back(metaPrimitive({ChId.rawId(),
0523 SL1metaPrimitive->t0,
0524 SL1metaPrimitive->x,
0525 SL1metaPrimitive->tanPhi,
0526 SL1metaPrimitive->phi,
0527 SL1metaPrimitive->phiB,
0528 SL1metaPrimitive->phi_cmssw,
0529 SL1metaPrimitive->phiB_cmssw,
0530 SL1metaPrimitive->chi2,
0531 new_quality,
0532 SL1metaPrimitive->wi1,
0533 SL1metaPrimitive->tdc1,
0534 SL1metaPrimitive->lat1,
0535 SL1metaPrimitive->wi2,
0536 SL1metaPrimitive->tdc2,
0537 SL1metaPrimitive->lat2,
0538 SL1metaPrimitive->wi3,
0539 SL1metaPrimitive->tdc3,
0540 SL1metaPrimitive->lat3,
0541 SL1metaPrimitive->wi4,
0542 SL1metaPrimitive->tdc4,
0543 SL1metaPrimitive->lat4,
0544 wi1,
0545 tdc1,
0546 lat1,
0547 wi2,
0548 tdc2,
0549 lat2,
0550 wi3,
0551 tdc3,
0552 lat3,
0553 wi4,
0554 tdc4,
0555 lat4,
0556 -1}));
0557 useFitSL1[sl1] = true;
0558 at_least_one_SL1_confirmation = true;
0559 }
0560 }
0561 }
0562
0563
0564
0565
0566 sl3 = 0;
0567 for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end();
0568 ++SL3metaPrimitive, sl3++) {
0569 if (useFitSL3[sl3])
0570 continue;
0571 if ((at_least_one_correlation == false || clean_chi2_correlation_) &&
0572 allow_confirmation_) {
0573
0574 int matched_digis = 0;
0575 double minx = minx_match_2digis_;
0576 double min2x = minx_match_2digis_;
0577 int best_tdc = -1;
0578 int next_tdc = -1;
0579 int best_wire = -1;
0580 int next_wire = -1;
0581 int best_layer = -1;
0582 int next_layer = -1;
0583 int best_lat = -1;
0584 int next_lat = -1;
0585 int lat = -1;
0586
0587 for (const auto &dtLayerId_It : *dtdigis) {
0588 const DTLayerId dtLId = dtLayerId_It.first;
0589
0590 const DTSuperLayerId dtSLId(dtLId.wheel(), dtLId.station(), dtLId.sector(), 1);
0591 if (dtSLId.rawId() != sl1Id.rawId())
0592 continue;
0593 double l_shift = 0;
0594 if (dtLId.layer() == 4)
0595 l_shift = X_POS_L4;
0596 if (dtLId.layer() == 3)
0597 l_shift = X_POS_L3;
0598 if (dtLId.layer() == 2)
0599 l_shift = -1 * X_POS_L3;
0600 if (dtLId.layer() == 1)
0601 l_shift = -1 * X_POS_L4;
0602 double x_inSL1 = SL3metaPrimitive->x + SL3metaPrimitive->tanPhi * (VERT_PHI1_PHI3 - l_shift);
0603 for (auto digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; ++digiIt) {
0604 DTWireId wireId(dtLId, (*digiIt).wire());
0605 double x_wire =
0606 shiftinfo_[wireId.rawId()] + ((*digiIt).time() - SL3metaPrimitive->t0) * DRIFT_SPEED / 10.;
0607 double x_wire_left =
0608 shiftinfo_[wireId.rawId()] - ((*digiIt).time() - SL3metaPrimitive->t0) * DRIFT_SPEED / 10.;
0609 lat = 1;
0610 if (std::abs(x_inSL1 - x_wire) > std::abs(x_inSL1 - x_wire_left)) {
0611 x_wire = x_wire_left;
0612 lat = 0;
0613 }
0614 if (std::abs(x_inSL1 - x_wire) < minx) {
0615 minx = std::abs(x_inSL1 - x_wire);
0616 next_wire = best_wire;
0617 next_tdc = best_tdc;
0618 next_layer = best_layer;
0619 next_lat = best_lat;
0620
0621 best_wire = (*digiIt).wire();
0622 best_tdc = (*digiIt).time();
0623 best_layer = dtLId.layer();
0624 best_lat = lat;
0625 matched_digis++;
0626 } else if ((std::abs(x_inSL1 - x_wire) >= minx) && (std::abs(x_inSL1 - x_wire) < min2x)) {
0627 minx = std::abs(x_inSL1 - x_wire);
0628 next_wire = (*digiIt).wire();
0629 next_tdc = (*digiIt).time();
0630 next_layer = dtLId.layer();
0631 next_lat = lat;
0632 matched_digis++;
0633 }
0634 }
0635 }
0636 if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) {
0637 int new_quality = CHIGHQ;
0638 if (SL3metaPrimitive->quality == LOWQ)
0639 new_quality = CLOWQ;
0640
0641 int wi1 = -1;
0642 int tdc1 = -1;
0643 int lat1 = -1;
0644 int wi2 = -1;
0645 int tdc2 = -1;
0646 int lat2 = -1;
0647 int wi3 = -1;
0648 int tdc3 = -1;
0649 int lat3 = -1;
0650 int wi4 = -1;
0651 int tdc4 = -1;
0652 int lat4 = -1;
0653
0654 if (next_layer == 1) {
0655 wi1 = next_wire;
0656 tdc1 = next_tdc;
0657 lat1 = next_lat;
0658 }
0659 if (next_layer == 2) {
0660 wi2 = next_wire;
0661 tdc2 = next_tdc;
0662 lat2 = next_lat;
0663 }
0664 if (next_layer == 3) {
0665 wi3 = next_wire;
0666 tdc3 = next_tdc;
0667 lat3 = next_lat;
0668 }
0669 if (next_layer == 4) {
0670 wi4 = next_wire;
0671 tdc4 = next_tdc;
0672 lat4 = next_lat;
0673 }
0674
0675 if (best_layer == 1) {
0676 wi1 = best_wire;
0677 tdc1 = best_tdc;
0678 lat1 = best_lat;
0679 }
0680 if (best_layer == 2) {
0681 wi2 = best_wire;
0682 tdc2 = best_tdc;
0683 lat2 = best_lat;
0684 }
0685 if (best_layer == 3) {
0686 wi3 = best_wire;
0687 tdc3 = best_tdc;
0688 lat3 = best_lat;
0689 }
0690 if (best_layer == 4) {
0691 wi4 = best_wire;
0692 tdc4 = best_tdc;
0693 lat4 = best_lat;
0694 }
0695
0696 if (!clean_chi2_correlation_)
0697 outMPaths.push_back(metaPrimitive({ChId.rawId(),
0698 SL3metaPrimitive->t0,
0699 SL3metaPrimitive->x,
0700 SL3metaPrimitive->tanPhi,
0701 SL3metaPrimitive->phi,
0702 SL3metaPrimitive->phiB,
0703 SL3metaPrimitive->phi_cmssw,
0704 SL3metaPrimitive->phiB_cmssw,
0705 SL3metaPrimitive->chi2,
0706 new_quality,
0707 wi1,
0708 tdc1,
0709 lat1,
0710 wi2,
0711 tdc2,
0712 lat2,
0713 wi3,
0714 tdc3,
0715 lat3,
0716 wi4,
0717 tdc4,
0718 lat4,
0719 SL3metaPrimitive->wi1,
0720 SL3metaPrimitive->tdc1,
0721 SL3metaPrimitive->lat1,
0722 SL3metaPrimitive->wi2,
0723 SL3metaPrimitive->tdc2,
0724 SL3metaPrimitive->lat2,
0725 SL3metaPrimitive->wi3,
0726 SL3metaPrimitive->tdc3,
0727 SL3metaPrimitive->lat3,
0728 SL3metaPrimitive->wi4,
0729 SL3metaPrimitive->tdc4,
0730 SL3metaPrimitive->lat4,
0731 -1}));
0732 else
0733 confirmedMetaPrimitives.push_back(metaPrimitive({ChId.rawId(),
0734 SL3metaPrimitive->t0,
0735 SL3metaPrimitive->x,
0736 SL3metaPrimitive->tanPhi,
0737 SL3metaPrimitive->phi,
0738 SL3metaPrimitive->phiB,
0739 SL3metaPrimitive->phi_cmssw,
0740 SL3metaPrimitive->phiB_cmssw,
0741 SL3metaPrimitive->chi2,
0742 new_quality,
0743 wi1,
0744 tdc1,
0745 lat1,
0746 wi2,
0747 tdc2,
0748 lat2,
0749 wi3,
0750 tdc3,
0751 lat3,
0752 wi4,
0753 tdc4,
0754 lat4,
0755 SL3metaPrimitive->wi1,
0756 SL3metaPrimitive->tdc1,
0757 SL3metaPrimitive->lat1,
0758 SL3metaPrimitive->wi2,
0759 SL3metaPrimitive->tdc2,
0760 SL3metaPrimitive->lat2,
0761 SL3metaPrimitive->wi3,
0762 SL3metaPrimitive->tdc3,
0763 SL3metaPrimitive->lat3,
0764 SL3metaPrimitive->wi4,
0765 SL3metaPrimitive->tdc4,
0766 SL3metaPrimitive->lat4,
0767 -1}));
0768 useFitSL3[sl3] = true;
0769 at_least_one_SL3_confirmation = true;
0770 }
0771 }
0772 }
0773
0774 if (clean_chi2_correlation_) {
0775 if (debug_)
0776 LogDebug("MuonPathAssociator") << "Pushing back correlated MPs to the MPs collection";
0777 removeSharingFits(chamberMetaPrimitives, outMPaths);
0778 }
0779 if (clean_chi2_correlation_) {
0780 if (debug_)
0781 LogDebug("MuonPathAssociator") << "Pushing back confirmed MPs to the complete vector";
0782 removeSharingHits(confirmedMetaPrimitives, chamberMetaPrimitives, outMPaths);
0783 }
0784
0785
0786 if (at_least_one_correlation == false || clean_chi2_correlation_) {
0787 if (debug_ && !at_least_one_correlation)
0788 LogDebug("MuonPathAssociator")
0789 << "correlation we found zero correlations, adding both collections as they are to the outMPaths";
0790 if (debug_)
0791 LogDebug("MuonPathAssociator")
0792 << "correlation sizes:" << SL1metaPrimitives.size() << " " << SL3metaPrimitives.size();
0793 if (at_least_one_SL1_confirmation == false || clean_chi2_correlation_) {
0794 sl1 = 0;
0795 for (auto SL1metaPrimitive = SL1metaPrimitives.begin(); SL1metaPrimitive != SL1metaPrimitives.end();
0796 ++SL1metaPrimitive, sl1++) {
0797 if (useFitSL1[sl1])
0798 continue;
0799
0800 DTSuperLayerId SLId(SL1metaPrimitive->rawId);
0801 DTChamberId(SLId.wheel(), SLId.station(), SLId.sector());
0802 metaPrimitive newSL1metaPrimitive = {ChId.rawId(),
0803 SL1metaPrimitive->t0,
0804 SL1metaPrimitive->x,
0805 SL1metaPrimitive->tanPhi,
0806 SL1metaPrimitive->phi,
0807 SL1metaPrimitive->phiB,
0808 SL1metaPrimitive->phi_cmssw,
0809 SL1metaPrimitive->phiB_cmssw,
0810 SL1metaPrimitive->chi2,
0811 SL1metaPrimitive->quality,
0812 SL1metaPrimitive->wi1,
0813 SL1metaPrimitive->tdc1,
0814 SL1metaPrimitive->lat1,
0815 SL1metaPrimitive->wi2,
0816 SL1metaPrimitive->tdc2,
0817 SL1metaPrimitive->lat2,
0818 SL1metaPrimitive->wi3,
0819 SL1metaPrimitive->tdc3,
0820 SL1metaPrimitive->lat3,
0821 SL1metaPrimitive->wi4,
0822 SL1metaPrimitive->tdc4,
0823 SL1metaPrimitive->lat4,
0824 -1,
0825 -1,
0826 -1,
0827 -1,
0828 -1,
0829 -1,
0830 -1,
0831 -1,
0832 -1,
0833 -1,
0834 -1,
0835 -1,
0836 -1};
0837
0838 bool ok = true;
0839 for (auto &metaPrimitive : chamberMetaPrimitives) {
0840 if (!isNotAPrimo(newSL1metaPrimitive, metaPrimitive)) {
0841 ok = false;
0842 break;
0843 }
0844 }
0845 if (!ok)
0846 continue;
0847
0848 if (!clean_chi2_correlation_)
0849 outMPaths.push_back(newSL1metaPrimitive);
0850 else
0851 normalMetaPrimitives.push_back(newSL1metaPrimitive);
0852 }
0853 }
0854 if (at_least_one_SL3_confirmation == false || clean_chi2_correlation_) {
0855 sl3 = 0;
0856 for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end();
0857 ++SL3metaPrimitive, sl3++) {
0858 if (useFitSL3[sl3])
0859 continue;
0860 DTSuperLayerId SLId(SL3metaPrimitive->rawId);
0861 DTChamberId(SLId.wheel(), SLId.station(), SLId.sector());
0862 metaPrimitive newSL3metaPrimitive = {ChId.rawId(),
0863 SL3metaPrimitive->t0,
0864 SL3metaPrimitive->x,
0865 SL3metaPrimitive->tanPhi,
0866 SL3metaPrimitive->phi,
0867 SL3metaPrimitive->phiB,
0868 SL3metaPrimitive->phi_cmssw,
0869 SL3metaPrimitive->phiB_cmssw,
0870 SL3metaPrimitive->chi2,
0871 SL3metaPrimitive->quality,
0872 -1,
0873 -1,
0874 -1,
0875 -1,
0876 -1,
0877 -1,
0878 -1,
0879 -1,
0880 -1,
0881 -1,
0882 -1,
0883 -1,
0884 SL3metaPrimitive->wi1,
0885 SL3metaPrimitive->tdc1,
0886 SL3metaPrimitive->lat1,
0887 SL3metaPrimitive->wi2,
0888 SL3metaPrimitive->tdc2,
0889 SL3metaPrimitive->lat2,
0890 SL3metaPrimitive->wi3,
0891 SL3metaPrimitive->tdc3,
0892 SL3metaPrimitive->lat3,
0893 SL3metaPrimitive->wi4,
0894 SL3metaPrimitive->tdc4,
0895 SL3metaPrimitive->lat4,
0896 -1};
0897
0898 if (!clean_chi2_correlation_)
0899 outMPaths.push_back(newSL3metaPrimitive);
0900 else
0901 normalMetaPrimitives.push_back(newSL3metaPrimitive);
0902 }
0903 }
0904 }
0905
0906 SL1metaPrimitives.clear();
0907 SL1metaPrimitives.erase(SL1metaPrimitives.begin(), SL1metaPrimitives.end());
0908 SL3metaPrimitives.clear();
0909 SL3metaPrimitives.erase(SL3metaPrimitives.begin(), SL3metaPrimitives.end());
0910
0911 vector<metaPrimitive> auxMetaPrimitives;
0912 if (clean_chi2_correlation_) {
0913 if (debug_)
0914 LogDebug("MuonPathAssociator") << "Pushing back normal MPs to the auxiliar vector";
0915 removeSharingHits(normalMetaPrimitives, confirmedMetaPrimitives, auxMetaPrimitives);
0916 }
0917 if (clean_chi2_correlation_) {
0918 if (debug_)
0919 LogDebug("MuonPathAssociator") << "Pushing back normal MPs to the MPs collection";
0920 removeSharingHits(auxMetaPrimitives, chamberMetaPrimitives, outMPaths);
0921 }
0922 }
0923 }
0924 }
0925
0926
0927 std::vector<metaPrimitive> SL2metaPrimitives;
0928
0929 for (int wh = -2; wh <= 2; wh++) {
0930 for (int st = 1; st <= 4; st++) {
0931 for (int se = 1; se <= 14; se++) {
0932 if (se >= 13 && st != 4)
0933 continue;
0934
0935 DTChamberId ChId(wh, st, se);
0936 DTSuperLayerId sl2Id(wh, st, se, 2);
0937
0938
0939 for (auto metaprimitiveIt = inMPaths.begin(); metaprimitiveIt != inMPaths.end(); ++metaprimitiveIt)
0940 if (metaprimitiveIt->rawId == sl2Id.rawId()) {
0941 SL2metaPrimitives.push_back(*metaprimitiveIt);
0942 printmPC(*metaprimitiveIt);
0943 outMPaths.push_back(*metaprimitiveIt);
0944 }
0945 }
0946 }
0947 }
0948
0949 LogDebug("MuonPathAssociator") << "\t etaTP: added " << SL2metaPrimitives.size() << "to outMPaths" << std::endl;
0950
0951 SL2metaPrimitives.clear();
0952 SL2metaPrimitives.erase(SL2metaPrimitives.begin(), SL2metaPrimitives.end());
0953 }
0954
0955 void MuonPathAssociator::removeSharingFits(vector<metaPrimitive> &chamberMPaths, vector<metaPrimitive> &allMPaths) {
0956 bool useFit[chamberMPaths.size()];
0957 for (unsigned int i = 0; i < chamberMPaths.size(); i++) {
0958 useFit[i] = true;
0959 }
0960 for (unsigned int i = 0; i < chamberMPaths.size(); i++) {
0961 if (debug_)
0962 LogDebug("MuonPathAssociator") << "Looking at prim" << i;
0963 if (!useFit[i])
0964 continue;
0965 for (unsigned int j = i + 1; j < chamberMPaths.size(); j++) {
0966 if (debug_)
0967 LogDebug("MuonPathAssociator") << "Comparing with prim " << j;
0968 if (!useFit[j])
0969 continue;
0970 metaPrimitive first = chamberMPaths[i];
0971 metaPrimitive second = chamberMPaths[j];
0972 if (shareFit(first, second)) {
0973 if (first.quality > second.quality)
0974 useFit[j] = false;
0975 else if (first.quality < second.quality)
0976 useFit[i] = false;
0977 else {
0978 if (first.chi2 < second.chi2)
0979 useFit[j] = false;
0980 else {
0981 useFit[i] = false;
0982 break;
0983 }
0984 }
0985 }
0986 }
0987 if (useFit[i]) {
0988 if (debug_)
0989 printmPC(chamberMPaths[i]);
0990 allMPaths.push_back(chamberMPaths[i]);
0991 }
0992 }
0993 if (debug_)
0994 LogDebug("MuonPathAssociator") << "---Swapping chamber---";
0995 }
0996
0997 void MuonPathAssociator::removeSharingHits(std::vector<metaPrimitive> &firstMPaths,
0998 std::vector<metaPrimitive> &secondMPaths,
0999 std::vector<metaPrimitive> &allMPaths) {
1000 for (auto &firstMP : firstMPaths) {
1001 if (debug_)
1002 LogDebug("MuonPathAssociator") << "----------------------------------";
1003 if (debug_)
1004 LogDebug("MuonPathAssociator") << "Turn for ";
1005 if (debug_)
1006 printmPC(firstMP);
1007 bool ok = true;
1008 for (auto &secondMP : secondMPaths) {
1009 if (debug_)
1010 LogDebug("MuonPathAssociator") << "Comparing with ";
1011 if (debug_)
1012 printmPC(secondMP);
1013 if (!isNotAPrimo(firstMP, secondMP)) {
1014 ok = false;
1015 break;
1016 }
1017 }
1018 if (ok) {
1019 allMPaths.push_back(firstMP);
1020 if (debug_)
1021 printmPC(firstMP);
1022 }
1023 if (debug_)
1024 LogDebug("MuonPathAssociator") << "----------------------------------";
1025 }
1026 }
1027
1028 bool MuonPathAssociator::shareFit(metaPrimitive first, metaPrimitive second) {
1029 bool lay1 = (first.wi1 == second.wi1) && (first.tdc1 = second.tdc1);
1030 bool lay2 = (first.wi2 == second.wi2) && (first.tdc2 = second.tdc2);
1031 bool lay3 = (first.wi3 == second.wi3) && (first.tdc3 = second.tdc3);
1032 bool lay4 = (first.wi4 == second.wi4) && (first.tdc4 = second.tdc4);
1033 bool lay5 = (first.wi5 == second.wi5) && (first.tdc5 = second.tdc5);
1034 bool lay6 = (first.wi6 == second.wi6) && (first.tdc6 = second.tdc6);
1035 bool lay7 = (first.wi7 == second.wi7) && (first.tdc7 = second.tdc7);
1036 bool lay8 = (first.wi8 == second.wi8) && (first.tdc8 = second.tdc8);
1037
1038 if (lay1 && lay2 && lay3 && lay4) {
1039 if (lay5 || lay6 || lay7 || lay8)
1040 return true;
1041 else
1042 return false;
1043 } else if (lay5 && lay6 && lay7 && lay8) {
1044 if (lay1 || lay2 || lay3 || lay4)
1045 return true;
1046 else
1047 return false;
1048 } else
1049 return false;
1050 }
1051
1052 bool MuonPathAssociator::isNotAPrimo(metaPrimitive first, metaPrimitive second) {
1053 int hitsSL1 = (first.wi1 != -1) + (first.wi2 != -1) + (first.wi3 != -1) + (first.wi4 != -1);
1054 int hitsSL3 = (first.wi5 != -1) + (first.wi6 != -1) + (first.wi7 != -1) + (first.wi8 != -1);
1055
1056 bool lay1 = (first.wi1 == second.wi1) && (first.tdc1 = second.tdc1) && (first.wi1 != -1);
1057 bool lay2 = (first.wi2 == second.wi2) && (first.tdc2 = second.tdc2) && (first.wi2 != -1);
1058 bool lay3 = (first.wi3 == second.wi3) && (first.tdc3 = second.tdc3) && (first.wi3 != -1);
1059 bool lay4 = (first.wi4 == second.wi4) && (first.tdc4 = second.tdc4) && (first.wi4 != -1);
1060 bool lay5 = (first.wi5 == second.wi5) && (first.tdc5 = second.tdc5) && (first.wi5 != -1);
1061 bool lay6 = (first.wi6 == second.wi6) && (first.tdc6 = second.tdc6) && (first.wi6 != -1);
1062 bool lay7 = (first.wi7 == second.wi7) && (first.tdc7 = second.tdc7) && (first.wi7 != -1);
1063 bool lay8 = (first.wi8 == second.wi8) && (first.tdc8 = second.tdc8) && (first.wi8 != -1);
1064
1065 return (((!lay1 && !lay2 && !lay3 && !lay4) || hitsSL1 < 3) && ((!lay5 && !lay6 && !lay7 && !lay8) || hitsSL3 < 3));
1066 }
1067
1068 void MuonPathAssociator::printmPC(metaPrimitive mP) {
1069 DTChamberId ChId(mP.rawId);
1070 LogDebug("MuonPathAssociator") << ChId << "\t"
1071 << " " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 << " "
1072 << setw(2) << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(2)
1073 << left << mP.wi5 << " " << setw(2) << left << mP.wi6 << " " << setw(2) << left
1074 << mP.wi7 << " " << setw(2) << left << mP.wi8 << " " << setw(5) << left << mP.tdc1
1075 << " " << setw(5) << left << mP.tdc2 << " " << setw(5) << left << mP.tdc3 << " "
1076 << setw(5) << left << mP.tdc4 << " " << setw(5) << left << mP.tdc5 << " " << setw(5)
1077 << left << mP.tdc6 << " " << setw(5) << left << mP.tdc7 << " " << setw(5) << left
1078 << mP.tdc8 << " " << setw(2) << left << mP.lat1 << " " << setw(2) << left << mP.lat2
1079 << " " << setw(2) << left << mP.lat3 << " " << setw(2) << left << mP.lat4 << " "
1080 << setw(2) << left << mP.lat5 << " " << setw(2) << left << mP.lat6 << " " << setw(2)
1081 << left << mP.lat7 << " " << setw(2) << left << mP.lat8 << " " << setw(10) << right
1082 << mP.x << " " << setw(9) << left << mP.tanPhi << " " << setw(5) << left << mP.t0
1083 << " " << setw(13) << left << mP.chi2 << " \n";
1084 }
