File indexing completed on 2024-09-07 04:36:47
0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 #include <cmath>
0004 #include <memory>
0005
0006 using namespace edm;
0007 using namespace std;
0008 using namespace cmsdt;
0009
0010
0011
0012 MuonPathAnalyzerInChamber::MuonPathAnalyzerInChamber(const ParameterSet &pset,
0013 edm::ConsumesCollector &iC,
0014 std::shared_ptr<GlobalCoordsObtainer> &globalcoordsobtainer)
0015 : MuonPathAnalyzer(pset, iC),
0016 debug_(pset.getUntrackedParameter<bool>("debug")),
0017 chi2Th_(pset.getParameter<double>("chi2Th")),
0018 shift_filename_(pset.getParameter<edm::FileInPath>("shift_filename")),
0019 bxTolerance_(30),
0020 minQuality_(LOWQ),
0021 chiSquareThreshold_(50),
0022 minHits4Fit_(pset.getParameter<int>("minHits4Fit")),
0023 splitPathPerSL_(pset.getParameter<bool>("splitPathPerSL")) {
0024
0025
0026 if (debug_)
0027 LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzer: constructor";
0028
0029 setChiSquareThreshold(chi2Th_ * 10000.);
0030
0031
0032 int rawId;
0033 std::ifstream ifin3(shift_filename_.fullPath());
0034 double shift;
0035 if (ifin3.fail()) {
0036 throw cms::Exception("Missing Input File")
0037 << "MuonPathAnalyzerInChamber::MuonPathAnalyzerInChamber() - 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 MuonPathAnalyzerInChamber::~MuonPathAnalyzerInChamber() {
0049 if (debug_)
0050 LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzer: destructor";
0051 }
0052
0053
0054
0055
0056 void MuonPathAnalyzerInChamber::initialise(const edm::EventSetup &iEventSetup) {
0057 if (debug_)
0058 LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzerInChamber::initialiase";
0059
0060 auto geom = iEventSetup.getHandle(dtGeomH);
0061 dtGeo_ = geom.product();
0062 }
0063
0064 void MuonPathAnalyzerInChamber::run(edm::Event &iEvent,
0065 const edm::EventSetup &iEventSetup,
0066 MuonPathPtrs &muonpaths,
0067 MuonPathPtrs &outmuonpaths) {
0068 if (debug_)
0069 LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzerInChamber: run";
0070
0071
0072 int nMuonPath_counter = 0;
0073 for (auto muonpath = muonpaths.begin(); muonpath != muonpaths.end(); ++muonpath) {
0074 if (debug_) {
0075 LogDebug("MuonPathAnalyzerInChamber")
0076 << "Full path: " << nMuonPath_counter << " , " << muonpath->get()->nprimitives() << " , "
0077 << muonpath->get()->nprimitivesUp() << " , " << muonpath->get()->nprimitivesDown();
0078 }
0079 ++nMuonPath_counter;
0080
0081
0082 MuonPathPtr muonpathUp_ptr = std::make_shared<MuonPath>();
0083 muonpathUp_ptr->setNPrimitives(8);
0084 muonpathUp_ptr->setNPrimitivesUp(muonpath->get()->nprimitivesUp());
0085 muonpathUp_ptr->setNPrimitivesDown(0);
0086
0087 MuonPathPtr muonpathDown_ptr = std::make_shared<MuonPath>();
0088 muonpathDown_ptr->setNPrimitives(8);
0089 muonpathDown_ptr->setNPrimitivesUp(0);
0090 muonpathDown_ptr->setNPrimitivesDown(muonpath->get()->nprimitivesDown());
0091
0092 for (int n = 0; n < muonpath->get()->nprimitives(); ++n) {
0093 DTPrimitivePtr prim = muonpath->get()->primitive(n);
0094
0095 if (prim->superLayerId() == 3) {
0096 muonpathUp_ptr->setPrimitive(prim, n);
0097 }
0098
0099 else if (prim->superLayerId() == 1) {
0100 muonpathDown_ptr->setPrimitive(prim, n);
0101 }
0102
0103 else
0104 continue;
0105 }
0106
0107 analyze(*muonpath, outmuonpaths);
0108
0109 if (splitPathPerSL_) {
0110 if (muonpathUp_ptr->nprimitivesUp() > 1 && muonpath->get()->nprimitivesDown() > 0)
0111 analyze(muonpathUp_ptr, outmuonpaths);
0112
0113 if (muonpathDown_ptr->nprimitivesDown() > 1 && muonpath->get()->nprimitivesUp() > 0)
0114 analyze(muonpathDown_ptr, outmuonpaths);
0115 }
0116 }
0117 }
0118
0119 void MuonPathAnalyzerInChamber::finish() {
0120 if (debug_)
0121 LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzer: finish";
0122 };
0123
0124
0125
0126
0127 void MuonPathAnalyzerInChamber::analyze(MuonPathPtr &inMPath, MuonPathPtrs &outMPath) {
0128 if (debug_)
0129 LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analyze \t\t\t\t starts";
0130
0131
0132 if (debug_)
0133 LogDebug("MuonPathAnalyzerInChamber") << inMPath->nprimitives();
0134 auto mPath = std::make_shared<MuonPath>(inMPath);
0135
0136 if (debug_) {
0137 LogDebug("MuonPathAnalyzerInChamber") << "DTp2::analyze, looking at mPath: ";
0138 for (int i = 0; i < mPath->nprimitives(); i++)
0139 LogDebug("MuonPathAnalyzerInChamber")
0140 << mPath->primitive(i)->layerId() << " , " << mPath->primitive(i)->superLayerId() << " , "
0141 << mPath->primitive(i)->channelId() << " , " << mPath->primitive(i)->laterality();
0142 }
0143
0144 if (debug_)
0145 LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analyze \t\t\t\t\t is Analyzable? ";
0146 if (!mPath->isAnalyzable())
0147 return;
0148 if (debug_)
0149 LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analyze \t\t\t\t\t yes it is analyzable " << mPath->isAnalyzable();
0150
0151
0152 buildLateralities(mPath);
0153 setWirePosAndTimeInMP(mPath);
0154
0155 std::shared_ptr<MuonPath> mpAux;
0156 int bestI = -1;
0157 float best_chi2 = 99999.;
0158 int added_lat = 0;
0159
0160
0161 for (int i = 0; i < (int)lateralities_.size(); i++) {
0162 if (debug_)
0163 LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analyze \t\t\t\t\t Start with combination " << i;
0164
0165 int NTotalHits = NUM_LAYERS_2SL;
0166 float xwire[NUM_LAYERS_2SL];
0167 int present_layer[NUM_LAYERS_2SL];
0168 for (int ii = 0; ii < 8; ii++) {
0169 xwire[ii] = mPath->xWirePos(ii);
0170 if (xwire[ii] == 0) {
0171 present_layer[ii] = 0;
0172 NTotalHits--;
0173 } else {
0174 present_layer[ii] = 1;
0175 }
0176 }
0177
0178 while (NTotalHits >= minHits4Fit_) {
0179 mPath->setChiSquare(0);
0180 calculateFitParameters(mPath, lateralities_[i], present_layer, added_lat);
0181 if (mPath->chiSquare() != 0)
0182 break;
0183 NTotalHits--;
0184 }
0185
0186 if (mPath->chiSquare() > chiSquareThreshold_)
0187 continue;
0188
0189 evaluateQuality(mPath);
0190
0191 if (mPath->quality() < minQuality_)
0192 continue;
0193
0194 double z = 0;
0195 double jm_x = (mPath->horizPos());
0196 int selected_Id = 0;
0197 for (int i = 0; i < mPath->nprimitives(); i++) {
0198 if (mPath->primitive(i)->isValidTime()) {
0199 selected_Id = mPath->primitive(i)->cameraId();
0200 mPath->setRawId(selected_Id);
0201 break;
0202 }
0203 }
0204 DTLayerId thisLId(selected_Id);
0205 DTChamberId ChId(thisLId.wheel(), thisLId.station(), thisLId.sector());
0206
0207 if (thisLId.station() >= 3)
0208 z += Z_SHIFT_MB4;
0209
0210 DTSuperLayerId MuonPathSLId(thisLId.wheel(), thisLId.station(), thisLId.sector(), thisLId.superLayer());
0211
0212
0213 int hits_in_SL1 = 0;
0214 int hits_in_SL3 = 0;
0215 for (int i = 0; i < mPath->nprimitives(); i++) {
0216 if (mPath->primitive(i)->isValidTime()) {
0217 if (i <= 3)
0218 ++hits_in_SL1;
0219 else if (i > 3)
0220 ++hits_in_SL3;
0221 }
0222 }
0223
0224 int SL_for_LUT = 0;
0225
0226 GlobalPoint jm_x_cmssw_global;
0227 if (hits_in_SL1 > 2 && hits_in_SL3 <= 2) {
0228
0229 jm_x += mPath->tanPhi() * (11.1 + 0.65);
0230 jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(jm_x, 0., z + 11.75));
0231 SL_for_LUT = 1;
0232 } else if (hits_in_SL1 <= 2 && hits_in_SL3 > 2) {
0233
0234 jm_x -= mPath->tanPhi() * (11.1 + 0.65);
0235 jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(jm_x, 0., z - 11.75));
0236 SL_for_LUT = 3;
0237 } else if (hits_in_SL1 > 2 && hits_in_SL3 > 2) {
0238
0239 jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(jm_x, 0., z));
0240 } else {
0241
0242 continue;
0243 }
0244
0245
0246 if (edm::isNotFinite(jm_x))
0247 continue;
0248
0249
0250 mPath->setHorizPos(jm_x);
0251
0252
0253 int thisec = MuonPathSLId.sector();
0254 if (thisec == 13)
0255 thisec = 4;
0256 if (thisec == 14)
0257 thisec = 10;
0258
0259
0260 double phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1);
0261 double psi = atan(mPath->tanPhi());
0262 mPath->setPhiCMSSW(phi_cmssw);
0263 mPath->setPhiBCMSSW(hasPosRF(MuonPathSLId.wheel(), MuonPathSLId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw);
0264
0265
0266 double phi = -999.;
0267 double phiB = -999.;
0268 double x_lut, slope_lut;
0269 DTSuperLayerId MuonPathSLId1(thisLId.wheel(), thisLId.station(), thisLId.sector(), 1);
0270 DTSuperLayerId MuonPathSLId3(thisLId.wheel(), thisLId.station(), thisLId.sector(), 3);
0271 DTWireId wireId1(MuonPathSLId1, 2, 1);
0272 DTWireId wireId3(MuonPathSLId3, 2, 1);
0273
0274
0275 double shift_for_lut = 0.;
0276 if (SL_for_LUT == 1) {
0277 shift_for_lut = int(10 * shiftinfo_[wireId1.rawId()] * INCREASED_RES_POS_POW);
0278 } else if (SL_for_LUT == 3) {
0279 shift_for_lut = int(10 * shiftinfo_[wireId3.rawId()] * INCREASED_RES_POS_POW);
0280 } else {
0281 int shift_sl1 = int(round(shiftinfo_[wireId1.rawId()] * INCREASED_RES_POS_POW * 10));
0282 int shift_sl3 = int(round(shiftinfo_[wireId3.rawId()] * INCREASED_RES_POS_POW * 10));
0283 if (shift_sl1 < shift_sl3) {
0284 shift_for_lut = shift_sl1;
0285 } else
0286 shift_for_lut = shift_sl3;
0287 }
0288 x_lut = double(jm_x) * 10. * INCREASED_RES_POS_POW - shift_for_lut;
0289 slope_lut = -(double(mPath->tanPhi()) * INCREASED_RES_SLOPE_POW);
0290
0291 auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), SL_for_LUT, x_lut, slope_lut);
0292 phi = global_coords[0];
0293 phiB = global_coords[1];
0294 mPath->setPhi(phi);
0295 mPath->setPhiB(phiB);
0296
0297 if (mPath->chiSquare() < best_chi2 && mPath->chiSquare() > 0) {
0298 mpAux = std::make_shared<MuonPath>(mPath);
0299 bestI = i;
0300 best_chi2 = mPath->chiSquare();
0301 }
0302 }
0303 if (mpAux != nullptr) {
0304 outMPath.push_back(std::move(mpAux));
0305 if (debug_)
0306 LogDebug("MuonPathAnalyzerInChamber")
0307 << "DTp2:analize \t\t\t\t\t Laterality " << bestI << " is the one with smaller chi2";
0308 } else {
0309 if (debug_)
0310 LogDebug("MuonPathAnalyzerInChamber")
0311 << "DTp2:analize \t\t\t\t\t No Laterality found with chi2 smaller than threshold";
0312 }
0313 if (debug_)
0314 LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analize \t\t\t\t\t Ended working with this set of lateralities";
0315 }
0316
0317 void MuonPathAnalyzerInChamber::setCellLayout(MuonPathPtr &mpath) {
0318 for (int i = 0; i <= mpath->nprimitives(); i++) {
0319 if (mpath->primitive(i)->isValidTime())
0320 cellLayout_[i] = mpath->primitive(i)->channelId();
0321 else
0322 cellLayout_[i] = -99;
0323 }
0324
0325
0326 mpath->setCellHorizontalLayout(cellLayout_);
0327 for (int i = 0; i <= mpath->nprimitives(); i++) {
0328 if (cellLayout_[i] >= 0) {
0329 mpath->setBaseChannelId(cellLayout_[i]);
0330 break;
0331 }
0332 }
0333 }
0334
0335
0336
0337
0338 void MuonPathAnalyzerInChamber::buildLateralities(MuonPathPtr &mpath) {
0339 if (debug_)
0340 LogDebug("MuonPathAnalyzerInChamber") << "MPAnalyzer::buildLateralities << setLateralitiesFromPrims ";
0341 mpath->setLateralCombFromPrimitives();
0342
0343 totalNumValLateralities_ = 0;
0344 lateralities_.clear();
0345 latQuality_.clear();
0346
0347
0348
0349 lateralities_.push_back(TLateralities());
0350 for (int ilat = 0; ilat < cmsdt::NUM_LAYERS_2SL; ilat++) {
0351
0352 LATERAL_CASES lr = (mpath->lateralComb())[ilat];
0353 if (debug_)
0354 LogDebug("MuonPathAnalyzerInChamber") << "[DEBUG_] Input[" << ilat << "]: " << lr;
0355
0356
0357 if (lr != NONE) {
0358 if (debug_)
0359 LogDebug("MuonPathAnalyzerInChamber") << "[DEBUG_] - Adding it to " << lateralities_.size() << " lists...";
0360 for (unsigned int iall = 0; iall < lateralities_.size(); iall++) {
0361 lateralities_[iall][ilat] = lr;
0362 }
0363 }
0364
0365 else {
0366
0367 auto ncurrentoptions = lateralities_.size();
0368
0369
0370 if (debug_)
0371 LogDebug("MuonPathAnalyzerInChamber") << "[DEBUG_] - Duplicating " << ncurrentoptions << " lists...";
0372 copy(lateralities_.begin(), lateralities_.end(), back_inserter(lateralities_));
0373 if (debug_)
0374 LogDebug("MuonPathAnalyzerInChamber") << "[DEBUG_] - Now we have " << lateralities_.size() << " lists...";
0375
0376
0377 for (unsigned int iall = 0; iall < ncurrentoptions; iall++) {
0378 lateralities_[iall][ilat] = LEFT;
0379 lateralities_[iall + ncurrentoptions][ilat] = RIGHT;
0380 }
0381 }
0382 }
0383
0384 totalNumValLateralities_ = (int)lateralities_.size();
0385 if (totalNumValLateralities_ > 128) {
0386
0387 LogDebug("MuonPathAnalyzerInChamber") << "[WARNING]: TOO MANY LATERALITIES TO CHECK !!";
0388 LogDebug("MuonPathAnalyzerInChamber") << "[WARNING]: skipping this muon";
0389 lateralities_.clear();
0390 latQuality_.clear();
0391 totalNumValLateralities_ = 0;
0392 }
0393
0394
0395 if (debug_) {
0396 for (unsigned int iall = 0; iall < lateralities_.size(); iall++) {
0397 LogDebug("MuonPathAnalyzerInChamber") << iall << " -> [";
0398 for (int ilat = 0; ilat < cmsdt::NUM_LAYERS_2SL; ilat++) {
0399 if (ilat != 0)
0400 LogDebug("MuonPathAnalyzerInChamber") << ",";
0401 LogDebug("MuonPathAnalyzerInChamber") << lateralities_[iall][ilat];
0402 }
0403 LogDebug("MuonPathAnalyzerInChamber") << "]";
0404 }
0405 }
0406 }
0407
0408 void MuonPathAnalyzerInChamber::setLateralitiesInMP(MuonPathPtr &mpath, TLateralities lat) {
0409 for (int i = 0; i < 8; i++)
0410 mpath->primitive(i)->setLaterality(lat[i]);
0411 }
0412
0413 void MuonPathAnalyzerInChamber::setWirePosAndTimeInMP(MuonPathPtr &mpath) {
0414 int selected_Id = 0;
0415 for (int i = 0; i < mpath->nprimitives(); i++) {
0416 if (mpath->primitive(i)->isValidTime()) {
0417 selected_Id = mpath->primitive(i)->cameraId();
0418 mpath->setRawId(selected_Id);
0419 break;
0420 }
0421 }
0422 DTLayerId thisLId(selected_Id);
0423 DTChamberId chId(thisLId.wheel(), thisLId.station(), thisLId.sector());
0424 if (debug_)
0425 LogDebug("MuonPathAnalyzerInChamber")
0426 << "Id " << chId.rawId() << " Wh " << chId.wheel() << " St " << chId.station() << " Se " << chId.sector();
0427 mpath->setRawId(chId.rawId());
0428
0429 DTSuperLayerId MuonPathSLId1(thisLId.wheel(), thisLId.station(), thisLId.sector(), 1);
0430 DTSuperLayerId MuonPathSLId3(thisLId.wheel(), thisLId.station(), thisLId.sector(), 3);
0431 DTWireId wireId1(MuonPathSLId1, 2, 1);
0432 DTWireId wireId3(MuonPathSLId3, 2, 1);
0433
0434 if (debug_)
0435 LogDebug("MuonPathAnalyzerInChamber")
0436 << "shift1=" << shiftinfo_[wireId1.rawId()] << " shift3=" << shiftinfo_[wireId3.rawId()];
0437
0438 float delta = 42000;
0439 float zwire[NUM_LAYERS_2SL] = {-13.7, -12.4, -11.1, -9.8002, 9.79999, 11.1, 12.4, 13.7};
0440 for (int i = 0; i < mpath->nprimitives(); i++) {
0441 if (mpath->primitive(i)->isValidTime()) {
0442 if (i < 4)
0443 mpath->setXWirePos(10000 * shiftinfo_[wireId1.rawId()] +
0444 (mpath->primitive(i)->channelId() + 0.5 * (double)((i + 1) % 2)) * delta,
0445 i);
0446 if (i >= 4)
0447 mpath->setXWirePos(10000 * shiftinfo_[wireId3.rawId()] +
0448 (mpath->primitive(i)->channelId() + 0.5 * (double)((i + 1) % 2)) * delta,
0449 i);
0450 mpath->setZWirePos(zwire[i] * 1000, i);
0451 mpath->setTWireTDC(mpath->primitive(i)->tdcTimeStamp() * DRIFT_SPEED, i);
0452 } else {
0453 mpath->setXWirePos(0., i);
0454 mpath->setZWirePos(0., i);
0455 mpath->setTWireTDC(-1 * DRIFT_SPEED, i);
0456 }
0457 if (debug_)
0458 LogDebug("MuonPathAnalyzerInChamber") << mpath->primitive(i)->tdcTimeStamp() << " ";
0459 }
0460 if (debug_)
0461 LogDebug("MuonPathAnalyzerInChamber");
0462 }
0463
0464 void MuonPathAnalyzerInChamber::calculateFitParameters(MuonPathPtr &mpath,
0465 TLateralities laterality,
0466 int present_layer[NUM_LAYERS_2SL],
0467 int &lat_added) {
0468
0469 int n_hits = 0;
0470 for (int l = 0; l < 8; ++l) {
0471 if (present_layer[l] == 1)
0472 n_hits++;
0473 }
0474
0475
0476 float xwire[NUM_LAYERS_2SL], zwire[NUM_LAYERS_2SL], tTDCvdrift[NUM_LAYERS_2SL];
0477 double b[NUM_LAYERS_2SL];
0478 for (int i = 0; i < 8; i++) {
0479 xwire[i] = mpath->xWirePos(i);
0480 zwire[i] = mpath->zWirePos(i);
0481 tTDCvdrift[i] = mpath->tWireTDC(i);
0482 b[i] = 1;
0483 }
0484
0485
0486
0487
0488 float xhit[NUM_LAYERS_2SL];
0489 for (int lay = 0; lay < 8; lay++) {
0490 if (debug_)
0491 LogDebug("MuonPathAnalyzerInChamber") << "In fitPerLat " << lay << " xwire " << xwire[lay] << " zwire "
0492 << zwire[lay] << " tTDCvdrift " << tTDCvdrift[lay];
0493 xhit[lay] = xwire[lay] + (-1 + 2 * laterality[lay]) * 1000 * tTDCvdrift[lay];
0494 if (debug_)
0495 LogDebug("MuonPathAnalyzerInChamber") << "In fitPerLat " << lay << " xhit " << xhit[lay];
0496 }
0497
0498
0499 double cbscal = 0.0;
0500 double zbscal = 0.0;
0501 double czscal = 0.0;
0502 double bbscal = 0.0;
0503 double zzscal = 0.0;
0504 double ccscal = 0.0;
0505
0506 for (int lay = 0; lay < 8; lay++) {
0507 if (present_layer[lay] == 0)
0508 continue;
0509 if (debug_)
0510 LogDebug("MuonPathAnalyzerInChamber")
0511 << " For layer " << lay + 1 << " xwire[lay] " << xwire[lay] << " zwire " << zwire[lay] << " b " << b[lay];
0512 if (debug_)
0513 LogDebug("MuonPathAnalyzerInChamber") << " xhit[lat][lay] " << xhit[lay];
0514 cbscal = (-1 + 2 * laterality[lay]) * b[lay] + cbscal;
0515 zbscal = zwire[lay] * b[lay] + zbscal;
0516 czscal = (-1 + 2 * laterality[lay]) * zwire[lay] + czscal;
0517
0518 bbscal = b[lay] * b[lay] + bbscal;
0519 zzscal = zwire[lay] * zwire[lay] + zzscal;
0520 ccscal = (-1 + 2 * laterality[lay]) * (-1 + 2 * laterality[lay]) + ccscal;
0521 }
0522
0523 double cz = 0.0;
0524 double cb = 0.0;
0525 double zb = 0.0;
0526 double zc = 0.0;
0527 double bc = 0.0;
0528 double bz = 0.0;
0529
0530 cz = (cbscal * zbscal - czscal * bbscal) / (zzscal * bbscal - zbscal * zbscal);
0531 cb = (czscal * zbscal - cbscal * zzscal) / (zzscal * bbscal - zbscal * zbscal);
0532
0533 zb = (czscal * cbscal - zbscal * ccscal) / (bbscal * ccscal - cbscal * cbscal);
0534 zc = (zbscal * cbscal - czscal * bbscal) / (bbscal * ccscal - cbscal * cbscal);
0535
0536 bc = (zbscal * czscal - cbscal * zzscal) / (ccscal * zzscal - czscal * czscal);
0537 bz = (cbscal * czscal - zbscal * ccscal) / (ccscal * zzscal - czscal * czscal);
0538
0539 double c_tilde[NUM_LAYERS_2SL];
0540 double z_tilde[NUM_LAYERS_2SL];
0541 double b_tilde[NUM_LAYERS_2SL];
0542
0543 for (int lay = 0; lay < 8; lay++) {
0544 if (present_layer[lay] == 0)
0545 continue;
0546 if (debug_)
0547 LogDebug("MuonPathAnalyzerInChamber")
0548 << " For layer " << lay + 1 << " xwire[lay] " << xwire[lay] << " zwire " << zwire[lay] << " b " << b[lay];
0549 c_tilde[lay] = (-1 + 2 * laterality[lay]) + cz * zwire[lay] + cb * b[lay];
0550 z_tilde[lay] = zwire[lay] + zb * b[lay] + zc * (-1 + 2 * laterality[lay]);
0551 b_tilde[lay] = b[lay] + bc * (-1 + 2 * laterality[lay]) + bz * zwire[lay];
0552 }
0553
0554
0555 double xctilde = 0.0;
0556 double xztilde = 0.0;
0557 double xbtilde = 0.0;
0558 double ctildectilde = 0.0;
0559 double ztildeztilde = 0.0;
0560 double btildebtilde = 0.0;
0561
0562 double rect0vdrift = 0.0;
0563 double recslope = 0.0;
0564 double recpos = 0.0;
0565
0566 for (int lay = 0; lay < 8; lay++) {
0567 if (present_layer[lay] == 0)
0568 continue;
0569 xctilde = xhit[lay] * c_tilde[lay] + xctilde;
0570 ctildectilde = c_tilde[lay] * c_tilde[lay] + ctildectilde;
0571 xztilde = xhit[lay] * z_tilde[lay] + xztilde;
0572 ztildeztilde = z_tilde[lay] * z_tilde[lay] + ztildeztilde;
0573 xbtilde = xhit[lay] * b_tilde[lay] + xbtilde;
0574 btildebtilde = b_tilde[lay] * b_tilde[lay] + btildebtilde;
0575 }
0576
0577
0578 rect0vdrift = xctilde / ctildectilde;
0579 recslope = xztilde / ztildeztilde;
0580 recpos = xbtilde / btildebtilde;
0581 if (debug_) {
0582 LogDebug("MuonPathAnalyzerInChamber") << " In fitPerLat Reconstructed values per lat "
0583 << " rect0vdrift " << rect0vdrift;
0584 LogDebug("MuonPathAnalyzerInChamber")
0585 << "rect0 " << rect0vdrift / DRIFT_SPEED << " recBX " << rect0vdrift / DRIFT_SPEED / 25 << " recslope "
0586 << recslope << " recpos " << recpos;
0587 }
0588
0589
0590 double rectdriftvdrift[NUM_LAYERS_2SL];
0591 double recres[NUM_LAYERS_2SL];
0592 double recchi2 = 0.0;
0593 int sign_tdriftvdrift = {0};
0594 int incell_tdriftvdrift = {0};
0595 int physical_slope = {0};
0596
0597
0598
0599
0600 double maxDif = -1;
0601 int maxInt = -1;
0602 int swap_laterality[8] = {0};
0603
0604 for (int lay = 0; lay < 8; lay++) {
0605 if (present_layer[lay] == 0)
0606 continue;
0607 rectdriftvdrift[lay] = tTDCvdrift[lay] - rect0vdrift / 1000;
0608 if (debug_)
0609 LogDebug("MuonPathAnalyzerInChamber") << rectdriftvdrift[lay];
0610 recres[lay] = xhit[lay] - zwire[lay] * recslope - b[lay] * recpos - (-1 + 2 * laterality[lay]) * rect0vdrift;
0611
0612
0613 if (abs(rectdriftvdrift[lay]) < 3) {
0614 swap_laterality[lay] = 1;
0615 }
0616
0617 if ((present_layer[lay] == 1) && (rectdriftvdrift[lay] < -0.1)) {
0618 sign_tdriftvdrift = -1;
0619 if (-0.1 - rectdriftvdrift[lay] > maxDif) {
0620 maxDif = -0.1 - rectdriftvdrift[lay];
0621 maxInt = lay;
0622 }
0623 }
0624 if ((present_layer[lay] == 1) && (abs(rectdriftvdrift[lay]) > 21.1)) {
0625 incell_tdriftvdrift = -1;
0626 if (rectdriftvdrift[lay] - 21.1 > maxDif) {
0627 maxDif = rectdriftvdrift[lay] - 21.1;
0628 maxInt = lay;
0629 }
0630 }
0631 }
0632
0633
0634
0635 if (lat_added == 0) {
0636 std::vector<TLateralities> additional_lateralities;
0637 additional_lateralities.clear();
0638 additional_lateralities.push_back(laterality);
0639
0640
0641 for (int swap = 0; swap < 8; ++swap) {
0642 if (swap_laterality[swap] == 1) {
0643 int add_lat_size = int(additional_lateralities.size());
0644 for (int ll = 0; ll < add_lat_size; ++ll) {
0645 TLateralities tmp_lat = additional_lateralities[ll];
0646 if (tmp_lat[swap] == LEFT)
0647 tmp_lat[swap] = RIGHT;
0648 else if (tmp_lat[swap] == RIGHT)
0649 tmp_lat[swap] = LEFT;
0650 else
0651 continue;
0652 additional_lateralities.push_back(tmp_lat);
0653 }
0654 }
0655 }
0656
0657
0658 int already_there = 0;
0659 for (int k = 0; k < int(additional_lateralities.size()); ++k) {
0660 already_there = 0;
0661 for (int j = 0; j < int(lateralities_.size()); ++j) {
0662 if (additional_lateralities[k] == lateralities_[j])
0663 already_there = 1;
0664 }
0665 if (already_there == 0) {
0666 lateralities_.push_back(additional_lateralities[k]);
0667 }
0668 }
0669 additional_lateralities.clear();
0670 lat_added = 1;
0671 }
0672
0673 if (fabs(recslope / 10) > 1.3)
0674 physical_slope = -1;
0675
0676 if (physical_slope == -1 && debug_)
0677 LogDebug("MuonPathAnalyzerInChamber") << "Combination with UNPHYSICAL slope ";
0678 if (sign_tdriftvdrift == -1 && debug_)
0679 LogDebug("MuonPathAnalyzerInChamber") << "Combination with negative tdrift-vdrift ";
0680 if (incell_tdriftvdrift == -1 && debug_)
0681 LogDebug("MuonPathAnalyzerInChamber") << "Combination with tdrift-vdrift larger than half cell ";
0682
0683 for (int lay = 0; lay < 8; lay++) {
0684 if (present_layer[lay] == 0)
0685 continue;
0686 recchi2 = recres[lay] * recres[lay] + recchi2;
0687 }
0688
0689 if (n_hits > 4)
0690 recchi2 = recchi2 / (n_hits - 3);
0691
0692 if (debug_)
0693 LogDebug("MuonPathAnalyzerInChamber")
0694 << "In fitPerLat Chi2 " << recchi2 << " with sign " << sign_tdriftvdrift << " within cell "
0695 << incell_tdriftvdrift << " physical_slope " << physical_slope;
0696
0697
0698 if (true && maxInt != -1) {
0699 present_layer[maxInt] = 0;
0700 if (debug_)
0701 LogDebug("MuonPathAnalyzerInChamber") << "We get rid of hit in layer " << maxInt;
0702 }
0703
0704
0705 if (!(sign_tdriftvdrift == -1) && !(incell_tdriftvdrift == -1) && !(physical_slope == -1)) {
0706 mpath->setBxTimeValue((rect0vdrift / DRIFT_SPEED) / 1000);
0707 mpath->setTanPhi(-1 * recslope / 10);
0708 mpath->setHorizPos(recpos / 10000);
0709 mpath->setChiSquare(recchi2 / 100000000);
0710 setLateralitiesInMP(mpath, laterality);
0711 if (debug_)
0712 LogDebug("MuonPathAnalyzerInChamber")
0713 << "In fitPerLat "
0714 << "t0 " << mpath->bxTimeValue() << " slope " << mpath->tanPhi() << " pos " << mpath->horizPos() << " chi2 "
0715 << mpath->chiSquare() << " rawId " << mpath->rawId();
0716 }
0717 }
0718 void MuonPathAnalyzerInChamber::evaluateQuality(MuonPathPtr &mPath) {
0719 mPath->setQuality(NOPATH);
0720
0721 int nPrimsUp(0), nPrimsDown(0);
0722 for (int i = 0; i < NUM_LAYERS_2SL; i++) {
0723 if (mPath->primitive(i)->isValidTime()) {
0724 if (i < 4)
0725 nPrimsDown++;
0726 else if (i >= 4)
0727 nPrimsUp++;
0728 }
0729 }
0730
0731 mPath->setNPrimitivesUp(nPrimsUp);
0732 mPath->setNPrimitivesDown(nPrimsDown);
0733
0734 if (mPath->nprimitivesUp() >= 4 && mPath->nprimitivesDown() >= 4) {
0735 mPath->setQuality(HIGHHIGHQ);
0736 } else if ((mPath->nprimitivesUp() == 4 && mPath->nprimitivesDown() == 3) ||
0737 (mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() == 4)) {
0738 mPath->setQuality(HIGHLOWQ);
0739 } else if ((mPath->nprimitivesUp() == 4 && mPath->nprimitivesDown() <= 2 && mPath->nprimitivesDown() > 0) ||
0740 (mPath->nprimitivesUp() <= 2 && mPath->nprimitivesUp() > 0 && mPath->nprimitivesDown() == 4)) {
0741 mPath->setQuality(CHIGHQ);
0742 } else if ((mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() == 3)) {
0743 mPath->setQuality(LOWLOWQ);
0744 } else if ((mPath->nprimitivesUp() == 3 && mPath->nprimitivesDown() <= 2 && mPath->nprimitivesDown() > 0) ||
0745 (mPath->nprimitivesUp() <= 2 && mPath->nprimitivesUp() > 0 && mPath->nprimitivesDown() == 3) ||
0746 (mPath->nprimitivesUp() == 2 && mPath->nprimitivesDown() == 2)) {
0747 mPath->setQuality(CLOWQ);
0748 } else if (mPath->nprimitivesUp() >= 4 || mPath->nprimitivesDown() >= 4) {
0749 mPath->setQuality(HIGHQ);
0750 } else if (mPath->nprimitivesUp() == 3 || mPath->nprimitivesDown() == 3) {
0751 mPath->setQuality(LOWQ);
0752 }
0753 }