Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:48 

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 // Constructors and destructor
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   // Obtention of parameters
0025 
0026   if (debug_)
0027     LogDebug("MuonPathAnalyzerInChamber") << "MuonPathAnalyzer: constructor";
0028 
0029   setChiSquareThreshold(chi2Th_ * 10000.);
0030 
0031   //shift
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 // Main methods (initialise, run, finish)
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   // fit per SL (need to allow for multiple outputs for a single mpath)
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     // Define muonpaths for up/down SL only
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       // UP
0095       if (prim->superLayerId() == 3) {
0096         muonpathUp_ptr->setPrimitive(prim, n);
0097       }
0098       // DOWN
0099       else if (prim->superLayerId() == 1) {
0100         muonpathDown_ptr->setPrimitive(prim, n);
0101       }
0102       // NOT UP NOR DOWN
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 //--- Private method
0126 //------------------------------------------------------------------
0127 void MuonPathAnalyzerInChamber::analyze(MuonPathPtr &inMPath, MuonPathPtrs &outMPath) {
0128   if (debug_)
0129     LogDebug("MuonPathAnalyzerInChamber") << "DTp2:analyze \t\t\t\t starts";
0130 
0131   // Clone the analyzed object
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   // first of all, get info from primitives, so we can reduce the number of latereralities:
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   // LOOP for all lateralities:
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     // Count hits in each SL
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     // Depending on which SL has hits, propagate jm_x to SL1, SL3, or to the center of the chamber
0226     GlobalPoint jm_x_cmssw_global;
0227     if (hits_in_SL1 > 2 && hits_in_SL3 <= 2) {
0228       // Uncorrelated or confirmed with 3 or 4 hits in SL1: propagate to SL1
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       // Uncorrelated or confirmed with 3 or 4 hits in SL3: propagate to SL3
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       // Correlated: stay at chamber center
0239       jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(jm_x, 0., z));
0240     } else {
0241       // Not interesting
0242       continue;
0243     }
0244 
0245     // Protection against non-converged fits
0246     if (edm::isNotFinite(jm_x))
0247       continue;
0248 
0249     // Updating muon-path horizontal position
0250     mPath->setHorizPos(jm_x);
0251 
0252     // Adjusting sector names for MB4
0253     int thisec = MuonPathSLId.sector();
0254     if (thisec == 13)
0255       thisec = 4;
0256     if (thisec == 14)
0257       thisec = 10;
0258 
0259     // Global coordinates from CMSSW
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     // Global coordinates from LUTs (firmware-like)
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     // use SL_for_LUT to decide the shift: x-axis origin for LUTs is left chamber side
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;  // position in cm * precision in JM RF
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   // putting info back into the mpath:
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  * For a combination of up to 8 cells, build all the lateralities to be tested,
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   /* We generate all the possible laterality combinations compatible with the built 
0348      group in the previous step*/
0349   lateralities_.push_back(TLateralities());
0350   for (int ilat = 0; ilat < cmsdt::NUM_LAYERS_2SL; ilat++) {
0351     // Get value from input
0352     LATERAL_CASES lr = (mpath->lateralComb())[ilat];
0353     if (debug_)
0354       LogDebug("MuonPathAnalyzerInChamber") << "[DEBUG_] Input[" << ilat << "]: " << lr;
0355 
0356     // If left/right fill number
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     // both possibilites
0365     else {
0366       // Get the number of possible options now
0367       auto ncurrentoptions = lateralities_.size();
0368 
0369       // Duplicate them
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       // Asign LEFT to first ncurrentoptions and RIGHT to the last
0377       for (unsigned int iall = 0; iall < ncurrentoptions; iall++) {
0378         lateralities_[iall][ilat] = LEFT;
0379         lateralities_[iall + ncurrentoptions][ilat] = RIGHT;
0380       }
0381     }  // else
0382   }    // Iterate over input array
0383 
0384   totalNumValLateralities_ = (int)lateralities_.size();
0385   if (totalNumValLateralities_ > 128) {
0386     // ADD PROTECTION!
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   // Dump values
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;                                                                      //um
0439   float zwire[NUM_LAYERS_2SL] = {-13.7, -12.4, -11.1, -9.8002, 9.79999, 11.1, 12.4, 13.7};  // mm
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);  // in um
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   // Get number of hits in current muonPath
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   // First prepare mpath for fit:
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   //// NOW Start FITTING:
0486 
0487   // fill hit position
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   //Proceed with calculation of fit parameters
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;  //it actually does not depend on laterality
0516     czscal = (-1 + 2 * laterality[lay]) * zwire[lay] + czscal;
0517 
0518     bbscal = b[lay] * b[lay] + bbscal;          //it actually does not depend on laterality
0519     zzscal = zwire[lay] * zwire[lay] + zzscal;  //it actually does not depend on laterality
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   //Calculate results per lat
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   // Results for t0vdrift (BX), slope and position per lat
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   //Get t*v and residuals per layer, and chi2 per laterality
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   // Select the worst hit in order to get rid of it
0598   // Also, check if any hits are close to the wire (rectdriftvdrift[lay] < 0.3 cm)
0599   // --> in that case, try also changing laterality of such hit
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     // If a hit is too close to the wire, set its corresponding "swap" flag to 1
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;  //Changed to 2.11 to account for resolution effects
0626       if (rectdriftvdrift[lay] - 21.1 > maxDif) {
0627         maxDif = rectdriftvdrift[lay] - 21.1;
0628         maxInt = lay;
0629       }
0630     }
0631   }
0632 
0633   // Now consider all possible alternative lateralities and push to lateralities_ those
0634   // we aren't considering yet
0635   if (lat_added == 0) {
0636     std::vector<TLateralities> additional_lateralities;
0637     additional_lateralities.clear();
0638     additional_lateralities.push_back(laterality);
0639     // Everytime the swap flag is 1, duplicate all the current elements
0640     // of additional_lateralities and swap their laterality
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     // Now compare all the additional lateralities with the lateralities we are considering:
0657     // if they are not there, add them
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   // Chi2/ndof
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   //LATERALITY IS NOT VALID
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   // LATERALITY IS VALID...
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!