Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:08 

0001 /**
0002  *  See header file for a description of this class.
0003  *  
0004  *  All the code is under revision
0005  *
0006  *
0007  *  \author A. Vitelli - INFN Torino, V.Palichik
0008  *  \author ported by: R. Bellan - INFN Torino
0009  */
0010 
0011 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedOrcaPatternRecognition.h"
0012 
0013 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
0014 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0015 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0016 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0017 #include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h"
0018 #include "DataFormats/Math/interface/deltaPhi.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 
0021 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0022 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0023 
0024 // Geometry
0025 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0026 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0027 
0028 // Framework
0029 #include "FWCore/Framework/interface/EventSetup.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "FWCore/Framework/interface/ESHandle.h"
0034 
0035 // C++
0036 #include <vector>
0037 
0038 using namespace std;
0039 
0040 const std::string metname = "Muon|RecoMuon|MuonSeedOrcaPatternRecognition";
0041 
0042 // Constructor
0043 MuonSeedOrcaPatternRecognition::MuonSeedOrcaPatternRecognition(const edm::ParameterSet& pset,
0044                                                                edm::ConsumesCollector& iC)
0045     : MuonSeedVPatternRecognition(pset),
0046       theCrackEtas(pset.getParameter<std::vector<double> >("crackEtas")),
0047       theCrackWindow(pset.getParameter<double>("crackWindow")),
0048       theDeltaPhiWindow(
0049           pset.existsAs<double>("deltaPhiSearchWindow") ? pset.getParameter<double>("deltaPhiSearchWindow") : 0.25),
0050       theDeltaEtaWindow(
0051           pset.existsAs<double>("deltaEtaSearchWindow") ? pset.getParameter<double>("deltaEtaSearchWindow") : 0.2),
0052       theDeltaCrackWindow(pset.existsAs<double>("deltaEtaCrackSearchWindow")
0053                               ? pset.getParameter<double>("deltaEtaCrackSearchWindow")
0054                               : 0.25),
0055       muonLayersToken(iC.esConsumes<MuonDetLayerGeometry, MuonRecoGeometryRecord>()) {
0056   muonMeasurements = new MuonDetLayerMeasurements(theDTRecSegmentLabel.label(),
0057                                                   theCSCRecSegmentLabel,
0058                                                   edm::InputTag(),
0059                                                   edm::InputTag(),
0060                                                   theME0RecSegmentLabel,
0061                                                   iC,
0062                                                   enableDTMeasurement,
0063                                                   enableCSCMeasurement,
0064                                                   false,
0065                                                   false,
0066                                                   enableME0Measurement);
0067 }
0068 
0069 // reconstruct muon's seeds
0070 void MuonSeedOrcaPatternRecognition::produce(const edm::Event& event,
0071                                              const edm::EventSetup& eSetup,
0072                                              std::vector<MuonRecHitContainer>& result) {
0073   // divide the RecHits by DetLayer, in order to fill the
0074   // RecHitContainer like it was in ORCA
0075 
0076   // Muon Geometry - DT, CSC, RPC and ME0
0077   edm::ESHandle<MuonDetLayerGeometry> muonLayers = eSetup.getHandle(muonLayersToken);
0078 
0079   // get the DT layers
0080   vector<const DetLayer*> dtLayers = muonLayers->allDTLayers();
0081 
0082   // get the CSC layers
0083   vector<const DetLayer*> cscForwardLayers = muonLayers->forwardCSCLayers();
0084   vector<const DetLayer*> cscBackwardLayers = muonLayers->backwardCSCLayers();
0085 
0086   // get the ME0 layers
0087   vector<const DetLayer*> me0ForwardLayers = muonLayers->forwardME0Layers();
0088   vector<const DetLayer*> me0BackwardLayers = muonLayers->backwardME0Layers();
0089 
0090   // Backward (z<0) EndCap disk
0091   const DetLayer* ME4Bwd = cscBackwardLayers[4];
0092   const DetLayer* ME3Bwd = cscBackwardLayers[3];
0093   const DetLayer* ME2Bwd = cscBackwardLayers[2];
0094   const DetLayer* ME12Bwd = cscBackwardLayers[1];
0095   const DetLayer* ME11Bwd = cscBackwardLayers[0];
0096 
0097   // Forward (z>0) EndCap disk
0098   const DetLayer* ME11Fwd = cscForwardLayers[0];
0099   const DetLayer* ME12Fwd = cscForwardLayers[1];
0100   const DetLayer* ME2Fwd = cscForwardLayers[2];
0101   const DetLayer* ME3Fwd = cscForwardLayers[3];
0102   const DetLayer* ME4Fwd = cscForwardLayers[4];
0103 
0104   // barrel
0105   const DetLayer* MB4DL = dtLayers[3];
0106   const DetLayer* MB3DL = dtLayers[2];
0107   const DetLayer* MB2DL = dtLayers[1];
0108   const DetLayer* MB1DL = dtLayers[0];
0109 
0110   // instantiate the accessor
0111   // Don not use RPC for seeding
0112 
0113   double barreldThetaCut = 0.2;
0114   // still lose good muons to a tighter cut
0115   double endcapdThetaCut = 1.0;
0116 
0117   MuonRecHitContainer list9 = filterSegments(muonMeasurements->recHits(MB4DL, event), barreldThetaCut);
0118   MuonRecHitContainer list6 = filterSegments(muonMeasurements->recHits(MB3DL, event), barreldThetaCut);
0119   MuonRecHitContainer list7 = filterSegments(muonMeasurements->recHits(MB2DL, event), barreldThetaCut);
0120   MuonRecHitContainer list8 = filterSegments(muonMeasurements->recHits(MB1DL, event), barreldThetaCut);
0121 
0122   dumpLayer("MB4 ", list9);
0123   dumpLayer("MB3 ", list6);
0124   dumpLayer("MB2 ", list7);
0125   dumpLayer("MB1 ", list8);
0126 
0127   bool* MB1 = zero(list8.size());
0128   bool* MB2 = zero(list7.size());
0129   bool* MB3 = zero(list6.size());
0130 
0131   MuonRecHitContainer muRH_ME0Fwd, muRH_ME0Bwd;
0132 
0133   if (!me0ForwardLayers.empty()) {  // Forward (z>0) EndCap disk
0134     const DetLayer* ME0Fwd = me0ForwardLayers[0];
0135     muRH_ME0Fwd = filterSegments(muonMeasurements->recHits(ME0Fwd, event), endcapdThetaCut);
0136   }
0137   if (!me0BackwardLayers.empty()) {  // Backward (z<0) EndCap disk
0138     const DetLayer* ME0Bwd = me0BackwardLayers[0];
0139     muRH_ME0Bwd = filterSegments(muonMeasurements->recHits(ME0Bwd, event), endcapdThetaCut);
0140   }
0141 
0142   endcapPatterns(filterSegments(muonMeasurements->recHits(ME11Bwd, event), endcapdThetaCut),
0143                  filterSegments(muonMeasurements->recHits(ME12Bwd, event), endcapdThetaCut),
0144                  filterSegments(muonMeasurements->recHits(ME2Bwd, event), endcapdThetaCut),
0145                  filterSegments(muonMeasurements->recHits(ME3Bwd, event), endcapdThetaCut),
0146                  filterSegments(muonMeasurements->recHits(ME4Bwd, event), endcapdThetaCut),
0147                  muRH_ME0Bwd,
0148                  list8,
0149                  list7,
0150                  list6,
0151                  MB1,
0152                  MB2,
0153                  MB3,
0154                  result);
0155 
0156   endcapPatterns(filterSegments(muonMeasurements->recHits(ME11Fwd, event), endcapdThetaCut),
0157                  filterSegments(muonMeasurements->recHits(ME12Fwd, event), endcapdThetaCut),
0158                  filterSegments(muonMeasurements->recHits(ME2Fwd, event), endcapdThetaCut),
0159                  filterSegments(muonMeasurements->recHits(ME3Fwd, event), endcapdThetaCut),
0160                  filterSegments(muonMeasurements->recHits(ME4Fwd, event), endcapdThetaCut),
0161                  muRH_ME0Fwd,
0162                  list8,
0163                  list7,
0164                  list6,
0165                  MB1,
0166                  MB2,
0167                  MB3,
0168                  result);
0169 
0170   // ----------    Barrel only
0171 
0172   unsigned int counter = 0;
0173   if (list9.size() < 100) {  // +v
0174     for (MuonRecHitContainer::iterator iter = list9.begin(); iter != list9.end(); iter++) {
0175       MuonRecHitContainer seedSegments;
0176       seedSegments.push_back(*iter);
0177       complete(seedSegments, list6, MB3);
0178       complete(seedSegments, list7, MB2);
0179       complete(seedSegments, list8, MB1);
0180       if (check(seedSegments))
0181         result.push_back(seedSegments);
0182     }
0183   }
0184 
0185   if (list6.size() < 100) {  // +v
0186     for (counter = 0; counter < list6.size(); counter++) {
0187       if (!MB3[counter]) {
0188         MuonRecHitContainer seedSegments;
0189         seedSegments.push_back(list6[counter]);
0190         complete(seedSegments, list7, MB2);
0191         complete(seedSegments, list8, MB1);
0192         complete(seedSegments, list9);
0193         if (check(seedSegments))
0194           result.push_back(seedSegments);
0195       }
0196     }
0197   }
0198 
0199   if (list7.size() < 100) {  // +v
0200     for (counter = 0; counter < list7.size(); counter++) {
0201       if (!MB2[counter]) {
0202         MuonRecHitContainer seedSegments;
0203         seedSegments.push_back(list7[counter]);
0204         complete(seedSegments, list8, MB1);
0205         complete(seedSegments, list9);
0206         complete(seedSegments, list6, MB3);
0207         if (seedSegments.size() > 1 || (seedSegments.size() == 1 && seedSegments[0]->dimension() == 4)) {
0208           result.push_back(seedSegments);
0209         }
0210       }
0211     }
0212   }
0213 
0214   if (list8.size() < 100) {  // +v
0215     for (counter = 0; counter < list8.size(); counter++) {
0216       if (!MB1[counter]) {
0217         MuonRecHitContainer seedSegments;
0218         seedSegments.push_back(list8[counter]);
0219         complete(seedSegments, list9);
0220         complete(seedSegments, list6, MB3);
0221         complete(seedSegments, list7, MB2);
0222         if (seedSegments.size() > 1 || (seedSegments.size() == 1 && seedSegments[0]->dimension() == 4)) {
0223           result.push_back(seedSegments);
0224         }
0225       }
0226     }
0227   }
0228 
0229   if (MB3)
0230     delete[] MB3;
0231   if (MB2)
0232     delete[] MB2;
0233   if (MB1)
0234     delete[] MB1;
0235 
0236   if (result.empty()) {
0237     // be a little stricter with single segment seeds
0238     barreldThetaCut = 0.2;
0239     endcapdThetaCut = 0.2;
0240 
0241     MuonRecHitContainer all = muonMeasurements->recHits(ME4Bwd, event);
0242     MuonRecHitContainer tmp = filterSegments(muonMeasurements->recHits(ME3Bwd, event), endcapdThetaCut);
0243     copy(tmp.begin(), tmp.end(), back_inserter(all));
0244 
0245     tmp = filterSegments(muonMeasurements->recHits(ME2Bwd, event), endcapdThetaCut);
0246     copy(tmp.begin(), tmp.end(), back_inserter(all));
0247 
0248     tmp = filterSegments(muonMeasurements->recHits(ME12Bwd, event), endcapdThetaCut);
0249     copy(tmp.begin(), tmp.end(), back_inserter(all));
0250 
0251     tmp = filterSegments(muonMeasurements->recHits(ME11Bwd, event), endcapdThetaCut);
0252     copy(tmp.begin(), tmp.end(), back_inserter(all));
0253 
0254     if (!me0BackwardLayers.empty()) {
0255       const DetLayer* ME0Bwd = me0BackwardLayers[0];
0256       tmp = filterSegments(muonMeasurements->recHits(ME0Bwd, event), endcapdThetaCut);
0257       copy(tmp.begin(), tmp.end(), back_inserter(all));
0258     }
0259     if (!me0ForwardLayers.empty()) {
0260       const DetLayer* ME0Fwd = me0ForwardLayers[0];
0261       tmp = filterSegments(muonMeasurements->recHits(ME0Fwd, event), endcapdThetaCut);
0262       copy(tmp.begin(), tmp.end(), back_inserter(all));
0263     }
0264 
0265     tmp = filterSegments(muonMeasurements->recHits(ME11Fwd, event), endcapdThetaCut);
0266     copy(tmp.begin(), tmp.end(), back_inserter(all));
0267 
0268     tmp = filterSegments(muonMeasurements->recHits(ME12Fwd, event), endcapdThetaCut);
0269     copy(tmp.begin(), tmp.end(), back_inserter(all));
0270 
0271     tmp = filterSegments(muonMeasurements->recHits(ME2Fwd, event), endcapdThetaCut);
0272     copy(tmp.begin(), tmp.end(), back_inserter(all));
0273 
0274     tmp = filterSegments(muonMeasurements->recHits(ME3Fwd, event), endcapdThetaCut);
0275     copy(tmp.begin(), tmp.end(), back_inserter(all));
0276 
0277     tmp = filterSegments(muonMeasurements->recHits(ME4Fwd, event), endcapdThetaCut);
0278     copy(tmp.begin(), tmp.end(), back_inserter(all));
0279 
0280     tmp = filterSegments(muonMeasurements->recHits(MB4DL, event), barreldThetaCut);
0281     copy(tmp.begin(), tmp.end(), back_inserter(all));
0282 
0283     tmp = filterSegments(muonMeasurements->recHits(MB3DL, event), barreldThetaCut);
0284     copy(tmp.begin(), tmp.end(), back_inserter(all));
0285 
0286     tmp = filterSegments(muonMeasurements->recHits(MB2DL, event), barreldThetaCut);
0287     copy(tmp.begin(), tmp.end(), back_inserter(all));
0288 
0289     tmp = filterSegments(muonMeasurements->recHits(MB1DL, event), barreldThetaCut);
0290     copy(tmp.begin(), tmp.end(), back_inserter(all));
0291 
0292     LogTrace(metname) << "Number of segments: " << all.size();
0293 
0294     for (MuonRecHitContainer::const_iterator segmentItr = all.begin(); segmentItr != all.end(); ++segmentItr) {
0295       MuonRecHitContainer singleSegmentContainer;
0296       singleSegmentContainer.push_back(*segmentItr);
0297       result.push_back(singleSegmentContainer);
0298     }
0299   }
0300 }
0301 
0302 bool* MuonSeedOrcaPatternRecognition::zero(unsigned listSize) {
0303   bool* result = nullptr;
0304   if (listSize) {
0305     result = new bool[listSize];
0306     for (size_t i = 0; i < listSize; i++)
0307       result[i] = false;
0308   }
0309   return result;
0310 }
0311 
0312 void MuonSeedOrcaPatternRecognition::endcapPatterns(const MuonRecHitContainer& me11,
0313                                                     const MuonRecHitContainer& me12,
0314                                                     const MuonRecHitContainer& me2,
0315                                                     const MuonRecHitContainer& me3,
0316                                                     const MuonRecHitContainer& me4,
0317                                                     const MuonRecHitContainer& me0,
0318                                                     const MuonRecHitContainer& mb1,
0319                                                     const MuonRecHitContainer& mb2,
0320                                                     const MuonRecHitContainer& mb3,
0321                                                     bool* MB1,
0322                                                     bool* MB2,
0323                                                     bool* MB3,
0324                                                     std::vector<MuonRecHitContainer>& result) {
0325   dumpLayer("ME4 ", me4);
0326   dumpLayer("ME3 ", me3);
0327   dumpLayer("ME2 ", me2);
0328   dumpLayer("ME12 ", me12);
0329   dumpLayer("ME11 ", me11);
0330   dumpLayer("ME0 ", me0);
0331 
0332   std::vector<MuonRecHitContainer> patterns;
0333   MuonRecHitContainer crackSegments;
0334   rememberCrackSegments(me11, crackSegments);
0335   rememberCrackSegments(me12, crackSegments);
0336   rememberCrackSegments(me2, crackSegments);
0337   rememberCrackSegments(me3, crackSegments);
0338   rememberCrackSegments(me4, crackSegments);
0339   rememberCrackSegments(me0, crackSegments);
0340 
0341   const MuonRecHitContainer& list24 = me4;
0342   const MuonRecHitContainer& list23 = me3;
0343 
0344   const MuonRecHitContainer& list12 = me2;
0345 
0346   const MuonRecHitContainer& list22 = me12;
0347   MuonRecHitContainer list21 = me11;
0348   // add ME0 to ME1
0349   list21.reserve(list21.size() + me0.size());
0350   copy(me0.begin(), me0.end(), back_inserter(list21));
0351 
0352   MuonRecHitContainer list11 = list21;
0353   MuonRecHitContainer list5 = list22;
0354   MuonRecHitContainer list13 = list23;
0355   MuonRecHitContainer list4 = list24;
0356 
0357   if (list21.empty()) {
0358     list11 = list22;
0359     list5 = list21;
0360   }
0361 
0362   if (list24.size() < list23.size() && !list24.empty()) {
0363     list13 = list24;
0364     list4 = list23;
0365   }
0366 
0367   if (list23.empty()) {
0368     list13 = list24;
0369     list4 = list23;
0370   }
0371 
0372   MuonRecHitContainer list1 = list11;
0373   MuonRecHitContainer list2 = list12;
0374   MuonRecHitContainer list3 = list13;
0375 
0376   if (list12.empty()) {
0377     list3 = list12;
0378     if (list11.size() <= list13.size() && !list11.empty()) {
0379       list1 = list11;
0380       list2 = list13;
0381     } else {
0382       list1 = list13;
0383       list2 = list11;
0384     }
0385   }
0386 
0387   if (list13.empty()) {
0388     if (list11.size() <= list12.size() && !list11.empty()) {
0389       list1 = list11;
0390       list2 = list12;
0391     } else {
0392       list1 = list12;
0393       list2 = list11;
0394     }
0395   }
0396 
0397   if (!list12.empty() && !list13.empty()) {
0398     if (list11.size() <= list12.size() && list11.size() <= list13.size() && !list11.empty()) {  // ME 1
0399       if (list12.size() > list13.size()) {
0400         list2 = list13;
0401         list3 = list12;
0402       }
0403     } else if (list12.size() <= list13.size()) {  //  start with ME 2
0404       list1 = list12;
0405       if (list11.size() <= list13.size() && !list11.empty()) {
0406         list2 = list11;
0407         list3 = list13;
0408       } else {
0409         list2 = list13;
0410         list3 = list11;
0411       }
0412     } else {  //  start with ME 3
0413       list1 = list13;
0414       if (list11.size() <= list12.size() && !list11.empty()) {
0415         list2 = list11;
0416         list3 = list12;
0417       } else {
0418         list2 = list12;
0419         list3 = list11;
0420       }
0421     }
0422   }
0423 
0424   bool* ME2 = zero(list2.size());
0425   bool* ME3 = zero(list3.size());
0426   bool* ME4 = zero(list4.size());
0427   bool* ME5 = zero(list5.size());
0428 
0429   // creates list of compatible track segments
0430   for (MuonRecHitContainer::iterator iter = list1.begin(); iter != list1.end(); iter++) {
0431     if ((*iter)->recHits().size() < 4 && !list3.empty())
0432       continue;  // 3p.tr-seg. are not so good for starting
0433 
0434     MuonRecHitContainer seedSegments;
0435     seedSegments.push_back(*iter);
0436     complete(seedSegments, list2, ME2);
0437     complete(seedSegments, list3, ME3);
0438     complete(seedSegments, list4, ME4);
0439     complete(seedSegments, list5, ME5);
0440     complete(seedSegments, mb3, MB3);
0441     complete(seedSegments, mb2, MB2);
0442     complete(seedSegments, mb1, MB1);
0443     if (check(seedSegments))
0444       patterns.push_back(seedSegments);
0445   }
0446 
0447   unsigned int counter;
0448 
0449   for (counter = 0; counter < list2.size(); counter++) {
0450     if (!ME2[counter]) {
0451       MuonRecHitContainer seedSegments;
0452       seedSegments.push_back(list2[counter]);
0453       complete(seedSegments, list3, ME3);
0454       complete(seedSegments, list4, ME4);
0455       complete(seedSegments, list5, ME5);
0456       complete(seedSegments, mb3, MB3);
0457       complete(seedSegments, mb2, MB2);
0458       complete(seedSegments, mb1, MB1);
0459       if (check(seedSegments))
0460         patterns.push_back(seedSegments);
0461     }
0462   }
0463 
0464   if (list3.size() < 20) {  // +v
0465     for (counter = 0; counter < list3.size(); counter++) {
0466       if (!ME3[counter]) {
0467         MuonRecHitContainer seedSegments;
0468         seedSegments.push_back(list3[counter]);
0469         complete(seedSegments, list4, ME4);
0470         complete(seedSegments, list5, ME5);
0471         complete(seedSegments, mb3, MB3);
0472         complete(seedSegments, mb2, MB2);
0473         complete(seedSegments, mb1, MB1);
0474         if (check(seedSegments))
0475           patterns.push_back(seedSegments);
0476       }
0477     }
0478   }
0479 
0480   if (list4.size() < 20) {  // +v
0481     for (counter = 0; counter < list4.size(); counter++) {
0482       if (!ME4[counter]) {
0483         MuonRecHitContainer seedSegments;
0484         seedSegments.push_back(list4[counter]);
0485         complete(seedSegments, list5, ME5);
0486         complete(seedSegments, mb3, MB3);
0487         complete(seedSegments, mb2, MB2);
0488         complete(seedSegments, mb1, MB1);
0489         if (check(seedSegments))
0490           patterns.push_back(seedSegments);
0491       }
0492     }
0493   }
0494 
0495   if (ME5)
0496     delete[] ME5;
0497   if (ME4)
0498     delete[] ME4;
0499   if (ME3)
0500     delete[] ME3;
0501   if (ME2)
0502     delete[] ME2;
0503 
0504   if (!patterns.empty()) {
0505     result.insert(result.end(), patterns.begin(), patterns.end());
0506   } else {
0507     if (!crackSegments.empty()) {
0508       // make some single-segment seeds
0509       for (MuonRecHitContainer::const_iterator crackSegmentItr = crackSegments.begin();
0510            crackSegmentItr != crackSegments.end();
0511            ++crackSegmentItr) {
0512         MuonRecHitContainer singleSegmentPattern;
0513         singleSegmentPattern.push_back(*crackSegmentItr);
0514         result.push_back(singleSegmentPattern);
0515       }
0516     }
0517   }
0518 }
0519 
0520 void MuonSeedOrcaPatternRecognition::complete(MuonRecHitContainer& seedSegments,
0521                                               const MuonRecHitContainer& recHits,
0522                                               bool* used) const {
0523   MuonRecHitContainer good_rhit;
0524   MuonPatternRecoDumper theDumper;
0525   //+v get all rhits compatible with the seed on dEta/dPhi Glob.
0526   ConstMuonRecHitPointer first = seedSegments[0];  // first rechit of seed
0527   GlobalPoint ptg2 = first->globalPosition();      // its global pos +v
0528   for (unsigned nr = 0; nr < recHits.size(); ++nr) {
0529     const MuonRecHitPointer& recHit(recHits[nr]);
0530     GlobalPoint ptg1(recHit->globalPosition());
0531     float deta = fabs(ptg1.eta() - ptg2.eta());
0532     // Geom::Phi should keep it in the range [-pi, pi]
0533     float dphi = fabs(deltaPhi(ptg1.barePhi(), ptg2.barePhi()));
0534     // be a little more lenient in cracks
0535     bool crack = isCrack(recHit) || isCrack(first);
0536     //float detaWindow = 0.3;
0537     float detaWindow = crack ? theDeltaCrackWindow : theDeltaEtaWindow;
0538     if (deta > detaWindow || dphi > theDeltaPhiWindow) {
0539       continue;
0540     }  // +vvp!!!
0541 
0542     good_rhit.push_back(recHit);
0543     if (used)
0544       markAsUsed(nr, recHits, used);
0545   }  // recHits iter
0546 
0547   // select the best rhit among the compatible ones (based on Dphi Glob & Dir)
0548   MuonRecHitPointer best = bestMatch(first, good_rhit);
0549   if (best && best->isValid())
0550     seedSegments.push_back(best);
0551 }
0552 
0553 MuonSeedOrcaPatternRecognition::MuonRecHitPointer MuonSeedOrcaPatternRecognition::bestMatch(
0554     const ConstMuonRecHitPointer& first, MuonRecHitContainer& good_rhit) const {
0555   MuonRecHitPointer best = nullptr;
0556   if (good_rhit.size() == 1)
0557     return good_rhit[0];
0558   double bestDiscrim = 10000.;
0559   for (MuonRecHitContainer::iterator iter = good_rhit.begin(); iter != good_rhit.end(); iter++) {
0560     double discrim = discriminator(first, *iter);
0561     if (discrim < bestDiscrim) {
0562       bestDiscrim = discrim;
0563       best = *iter;
0564     }
0565   }
0566   return best;
0567 }
0568 
0569 double MuonSeedOrcaPatternRecognition::discriminator(const ConstMuonRecHitPointer& first,
0570                                                      MuonRecHitPointer& other) const {
0571   GlobalPoint gp1 = first->globalPosition();
0572   GlobalPoint gp2 = other->globalPosition();
0573   GlobalVector gd1 = first->globalDirection();
0574   GlobalVector gd2 = other->globalDirection();
0575   if (first->isDT() || other->isDT()) {
0576     return fabs(deltaPhi(gd1.barePhi(), gd2.barePhi()));
0577   }
0578 
0579   // penalize those 3-hit segments
0580   int nhits = other->recHits().size();
0581   int penalty = std::max(nhits - 2, 1);
0582   float dphig = deltaPhi(gp1.barePhi(), gp2.barePhi());
0583   // ME1A has slanted wires, so matching theta position doesn't work well.
0584   if (isME1A(first) || isME1A(other)) {
0585     return fabs(dphig / penalty);
0586   }
0587 
0588   float dthetag = gp1.theta() - gp2.theta();
0589   float dphid2 = fabs(deltaPhi(gd2.barePhi(), gp2.barePhi()));
0590   if (dphid2 > M_PI * .5)
0591     dphid2 = M_PI - dphid2;  //+v
0592   float dthetad2 = gp2.theta() - gd2.theta();
0593   // for CSC, make a big chi-squared of relevant variables
0594   // FIXME for 100 GeV mnd above muons, this doesn't perform as well as
0595   // previous methods.  needs investigation.
0596   float chisq = ((dphig / 0.02) * (dphig / 0.02) + (dthetag / 0.003) * (dthetag / 0.003) +
0597                  (dphid2 / 0.06) * (dphid2 / 0.06) + (dthetad2 / 0.08) * (dthetad2 / 0.08));
0598   return chisq / penalty;
0599 }
0600 
0601 bool MuonSeedOrcaPatternRecognition::check(const MuonRecHitContainer& segments) { return (segments.size() > 1); }
0602 
0603 void MuonSeedOrcaPatternRecognition::markAsUsed(int nr, const MuonRecHitContainer& recHits, bool* used) const {
0604   used[nr] = true;
0605   // if it's ME1A with two other segments in the container, mark the ghosts as used, too.
0606   if (recHits[nr]->isCSC()) {
0607     CSCDetId detId(recHits[nr]->geographicalId().rawId());
0608     if (detId.ring() == 4) {
0609       std::vector<unsigned> chamberHitNs;
0610       for (unsigned i = 0; i < recHits.size(); ++i) {
0611         if (recHits[i]->geographicalId().rawId() == detId.rawId()) {
0612           chamberHitNs.push_back(i);
0613         }
0614       }
0615       if (chamberHitNs.size() == 3) {
0616         for (unsigned i = 0; i < 3; ++i) {
0617           used[chamberHitNs[i]] = true;
0618         }
0619       }
0620     }
0621   }
0622 }
0623 
0624 bool MuonSeedOrcaPatternRecognition::isCrack(const ConstMuonRecHitPointer& segment) const {
0625   bool result = false;
0626   double absEta = fabs(segment->globalPosition().eta());
0627   for (std::vector<double>::const_iterator crackItr = theCrackEtas.begin(); crackItr != theCrackEtas.end();
0628        ++crackItr) {
0629     if (fabs(absEta - *crackItr) < theCrackWindow) {
0630       result = true;
0631     }
0632   }
0633   return result;
0634 }
0635 
0636 void MuonSeedOrcaPatternRecognition::rememberCrackSegments(const MuonRecHitContainer& segments,
0637                                                            MuonRecHitContainer& crackSegments) const {
0638   for (MuonRecHitContainer::const_iterator segmentItr = segments.begin(); segmentItr != segments.end(); ++segmentItr) {
0639     if ((**segmentItr).hit()->dimension() == 4 && isCrack(*segmentItr)) {
0640       crackSegments.push_back(*segmentItr);
0641     }
0642     // save ME0 segments if eta > 2.4, no other detectors
0643     if ((*segmentItr)->isME0() && std::abs((*segmentItr)->globalPosition().eta()) > 2.4) {
0644       crackSegments.push_back(*segmentItr);
0645     }
0646   }
0647 }
0648 
0649 void MuonSeedOrcaPatternRecognition::dumpLayer(const char* name, const MuonRecHitContainer& segments) const {
0650   MuonPatternRecoDumper theDumper;
0651 
0652   LogTrace(metname) << name << std::endl;
0653   for (MuonRecHitContainer::const_iterator segmentItr = segments.begin(); segmentItr != segments.end(); ++segmentItr) {
0654     LogTrace(metname) << theDumper.dumpMuonId((**segmentItr).geographicalId());
0655   }
0656 }
0657 
0658 MuonSeedOrcaPatternRecognition::MuonRecHitContainer MuonSeedOrcaPatternRecognition::filterSegments(
0659     const MuonRecHitContainer& segments, double dThetaCut) const {
0660   MuonPatternRecoDumper theDumper;
0661   MuonRecHitContainer result;
0662   for (MuonRecHitContainer::const_iterator segmentItr = segments.begin(); segmentItr != segments.end(); ++segmentItr) {
0663     double dtheta = (*segmentItr)->globalDirection().theta() - (*segmentItr)->globalPosition().theta();
0664     if ((*segmentItr)->isDT()) {
0665       // only apply the cut to 4D segments
0666       if ((*segmentItr)->dimension() == 2 || fabs(dtheta) < dThetaCut) {
0667         result.push_back(*segmentItr);
0668       } else {
0669         LogTrace(metname) << "Cutting segment " << theDumper.dumpMuonId((**segmentItr).geographicalId())
0670                           << " because dtheta = " << dtheta;
0671       }
0672 
0673     } else if ((*segmentItr)->isCSC()) {
0674       if (fabs(dtheta) < dThetaCut) {
0675         result.push_back(*segmentItr);
0676       } else {
0677         LogTrace(metname) << "Cutting segment " << theDumper.dumpMuonId((**segmentItr).geographicalId())
0678                           << " because dtheta = " << dtheta;
0679       }
0680     } else if ((*segmentItr)->isME0()) {
0681       if (fabs(dtheta) < dThetaCut) {
0682         result.push_back(*segmentItr);
0683       } else {
0684         LogTrace(metname) << "Cutting segment " << theDumper.dumpMuonId((**segmentItr).geographicalId())
0685                           << " because dtheta = " << dtheta;
0686       }
0687     }
0688   }
0689   filterOverlappingChambers(result);
0690   return result;
0691 }
0692 
0693 void MuonSeedOrcaPatternRecognition::filterOverlappingChambers(MuonRecHitContainer& segments) const {
0694   if (segments.empty())
0695     return;
0696   MuonPatternRecoDumper theDumper;
0697   // need to optimize cuts
0698   double dphiCut = 0.05;
0699   double detaCut = 0.05;
0700   std::vector<unsigned> toKill;
0701   std::vector<unsigned> me1aOverlaps;
0702   // loop over all segment pairs to see if there are two that match up in eta and phi
0703   // but from different chambers
0704   unsigned nseg = segments.size();
0705   for (unsigned i = 0; i < nseg - 1; ++i) {
0706     GlobalPoint pg1 = segments[i]->globalPosition();
0707     for (unsigned j = i + 1; j < nseg; ++j) {
0708       GlobalPoint pg2 = segments[j]->globalPosition();
0709       if (segments[i]->geographicalId().rawId() != segments[j]->geographicalId().rawId() &&
0710           fabs(deltaPhi(pg1.barePhi(), pg2.barePhi())) < dphiCut && fabs(pg1.eta() - pg2.eta()) < detaCut) {
0711         LogTrace(metname) << "OVERLAP " << theDumper.dumpMuonId(segments[i]->geographicalId()) << " "
0712                           << theDumper.dumpMuonId(segments[j]->geographicalId());
0713         // see which one is best
0714         toKill.push_back((countHits(segments[i]) < countHits(segments[j])) ? i : j);
0715         if (isME1A(segments[i])) {
0716           me1aOverlaps.push_back(i);
0717           me1aOverlaps.push_back(j);
0718         }
0719       }
0720     }
0721   }
0722 
0723   if (toKill.empty())
0724     return;
0725 
0726   // try to kill ghosts assigned to overlaps
0727   for (unsigned i = 0; i < me1aOverlaps.size(); ++i) {
0728     DetId detId(segments[me1aOverlaps[i]]->geographicalId());
0729     vector<unsigned> inSameChamber;
0730     for (unsigned j = 0; j < nseg; ++j) {
0731       if (i != j && segments[j]->geographicalId() == detId) {
0732         inSameChamber.push_back(j);
0733       }
0734     }
0735     if (inSameChamber.size() == 2) {
0736       toKill.push_back(inSameChamber[0]);
0737       toKill.push_back(inSameChamber[1]);
0738     }
0739   }
0740   // now kill the killable
0741   MuonRecHitContainer result;
0742   for (unsigned i = 0; i < nseg; ++i) {
0743     if (std::find(toKill.begin(), toKill.end(), i) == toKill.end()) {
0744       result.push_back(segments[i]);
0745     }
0746   }
0747   segments.swap(result);
0748 }
0749 
0750 bool MuonSeedOrcaPatternRecognition::isME1A(const ConstMuonRecHitPointer& segment) const {
0751   return segment->isCSC() && CSCDetId(segment->geographicalId()).ring() == 4;
0752 }
0753 
0754 int MuonSeedOrcaPatternRecognition::countHits(const MuonRecHitPointer& segment) const {
0755   int count = 0;
0756   vector<TrackingRecHit*> components = (*segment).recHits();
0757   for (vector<TrackingRecHit*>::const_iterator component = components.begin(); component != components.end();
0758        ++component) {
0759     int componentSize = (**component).recHits().size();
0760     count += (componentSize == 0) ? 1 : componentSize;
0761   }
0762   return count;
0763 }
[ 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!