Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-17 02:07:40 

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