Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /**
0002  *  See header file for a description of this class.
0003  *  
0004  *  \author Shih-Chuan Kao, Dominique Fortin - UCR
0005  */
0006 
0007 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedBuilder.h>
0008 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCreator.h>
0009 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCleaner.h>
0010 
0011 // Data Formats
0012 #include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
0013 #include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
0014 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
0015 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0016 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
0017 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
0018 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
0019 #include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>
0020 
0021 // Geometry
0022 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0023 #include <TrackingTools/DetLayers/interface/DetLayer.h>
0024 #include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
0025 #include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
0026 #include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>
0027 
0028 // muon service
0029 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0030 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
0031 
0032 // Framework
0033 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include <FWCore/Framework/interface/ESHandle.h>
0036 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0037 #include <DataFormats/Common/interface/Handle.h>
0038 
0039 // C++
0040 #include <vector>
0041 #include <deque>
0042 #include <utility>
0043 
0044 //typedef std::pair<double, TrajectorySeed> seedpr ;
0045 //static bool ptDecreasing(const seedpr s1, const seedpr s2) { return ( s1.first > s2.first ); }
0046 //static bool lengthSorting(const TrajectorySeed s1, const TrajectorySeed s2) { return ( s1.nHits() > s2.nHits() ); }
0047 
0048 /*
0049  * Constructor
0050  */
0051 MuonSeedBuilder::MuonSeedBuilder(const edm::ParameterSet& pset, edm::ConsumesCollector& iC) {
0052   // Local Debug flag
0053   debug = pset.getParameter<bool>("DebugMuonSeed");
0054 
0055   // enable the DT chamber
0056   enableDTMeasurement = pset.getParameter<bool>("EnableDTMeasurement");
0057   theDTSegmentLabel = pset.getParameter<edm::InputTag>("DTSegmentLabel");
0058 
0059   // enable the CSC chamber
0060   enableCSCMeasurement = pset.getParameter<bool>("EnableCSCMeasurement");
0061   theCSCSegmentLabel = pset.getParameter<edm::InputTag>("CSCSegmentLabel");
0062 
0063   // Parameters for seed creation in endcap region
0064   minCSCHitsPerSegment = pset.getParameter<int>("minCSCHitsPerSegment");
0065   maxDeltaEtaCSC = pset.getParameter<double>("maxDeltaEtaCSC");
0066   maxDeltaPhiCSC = pset.getParameter<double>("maxDeltaPhiCSC");
0067 
0068   // Parameters for seed creation in overlap region
0069   maxDeltaEtaOverlap = pset.getParameter<double>("maxDeltaEtaOverlap");
0070   maxDeltaPhiOverlap = pset.getParameter<double>("maxDeltaPhiOverlap");
0071 
0072   // Parameters for seed creation in barrel region
0073   minDTHitsPerSegment = pset.getParameter<int>("minDTHitsPerSegment");
0074   maxDeltaEtaDT = pset.getParameter<double>("maxDeltaEtaDT");
0075   maxDeltaPhiDT = pset.getParameter<double>("maxDeltaPhiDT");
0076 
0077   // Parameters to suppress combinatorics (ghosts)
0078   maxEtaResolutionDT = pset.getParameter<double>("maxEtaResolutionDT");
0079   maxPhiResolutionDT = pset.getParameter<double>("maxPhiResolutionDT");
0080   maxEtaResolutionCSC = pset.getParameter<double>("maxEtaResolutionCSC");
0081   maxPhiResolutionCSC = pset.getParameter<double>("maxPhiResolutionCSC");
0082 
0083   // Class for forming muon seeds
0084   muonSeedCreate_ = new MuonSeedCreator(pset);
0085   muonSeedClean_ = new MuonSeedCleaner(pset, edm::ConsumesCollector(iC));
0086 
0087   // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
0088   muonMeasurements = new MuonDetLayerMeasurements(theDTSegmentLabel,
0089                                                   theCSCSegmentLabel,
0090                                                   edm::InputTag(),
0091                                                   edm::InputTag(),
0092                                                   edm::InputTag(),
0093                                                   iC,
0094                                                   enableDTMeasurement,
0095                                                   enableCSCMeasurement,
0096                                                   false,
0097                                                   false,
0098                                                   false);
0099 }
0100 
0101 /*
0102  * Destructor
0103  */
0104 MuonSeedBuilder::~MuonSeedBuilder() {
0105   delete muonSeedCreate_;
0106   delete muonSeedClean_;
0107   if (muonMeasurements)
0108     delete muonMeasurements;
0109 }
0110 
0111 /* 
0112  * build
0113  *
0114  * Where the segments are sorted out to make a protoTrack (vector of matching segments in different 
0115  * stations/layers).  The protoTrack is then passed to the seed creator to create CSC, DT or overlap seeds
0116  *
0117  */
0118 //void MuonSeedBuilder::build( MuonDetLayerMeasurements muonMeasurements, TrajectorySeedCollection& theSeeds ) {
0119 int MuonSeedBuilder::build(edm::Event& event, const edm::EventSetup& eventSetup, TrajectorySeedCollection& theSeeds) {
0120   // Pass the Magnetic Field to where the seed is create
0121   muonSeedCreate_->setBField(BField);
0122 
0123   // Create temporary collection of seeds which will be cleaned up to remove combinatorics
0124   std::vector<TrajectorySeed> rawSeeds;
0125   std::vector<float> etaOfSeed;
0126   std::vector<float> phiOfSeed;
0127   std::vector<int> nSegOnSeed;
0128 
0129   // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
0130   // MuonDetLayerMeasurements muonMeasurements(enableDTMeasurement,enableCSCMeasurement,false,
0131   //                                          theDTSegmentLabel.label(),theCSCSegmentLabel.label());
0132 
0133   // 1) Get the various stations and store segments in containers for each station (layers)
0134 
0135   // 1a. get the DT segments by stations (layers):
0136   std::vector<const DetLayer*> dtLayers = muonLayers->allDTLayers();
0137 
0138   SegmentContainer DTlist4 = muonMeasurements->recHits(dtLayers[3], event);
0139   SegmentContainer DTlist3 = muonMeasurements->recHits(dtLayers[2], event);
0140   SegmentContainer DTlist2 = muonMeasurements->recHits(dtLayers[1], event);
0141   SegmentContainer DTlist1 = muonMeasurements->recHits(dtLayers[0], event);
0142 
0143   // Initialize flags that a given segment has been allocated to a seed
0144   BoolContainer usedDTlist4(DTlist4.size(), false);
0145   BoolContainer usedDTlist3(DTlist3.size(), false);
0146   BoolContainer usedDTlist2(DTlist2.size(), false);
0147   BoolContainer usedDTlist1(DTlist1.size(), false);
0148 
0149   if (debug) {
0150     std::cout << "*** Number of DT segments is: " << DTlist4.size() + DTlist3.size() + DTlist2.size() + DTlist1.size()
0151               << std::endl;
0152     std::cout << "In MB1: " << DTlist1.size() << std::endl;
0153     std::cout << "In MB2: " << DTlist2.size() << std::endl;
0154     std::cout << "In MB3: " << DTlist3.size() << std::endl;
0155     std::cout << "In MB4: " << DTlist4.size() << std::endl;
0156   }
0157 
0158   // 1b. get the CSC segments by stations (layers):
0159   // 1b.1 Global z < 0
0160   std::vector<const DetLayer*> cscBackwardLayers = muonLayers->backwardCSCLayers();
0161   SegmentContainer CSClist4B = muonMeasurements->recHits(cscBackwardLayers[4], event);
0162   SegmentContainer CSClist3B = muonMeasurements->recHits(cscBackwardLayers[3], event);
0163   SegmentContainer CSClist2B = muonMeasurements->recHits(cscBackwardLayers[2], event);
0164   SegmentContainer CSClist1B = muonMeasurements->recHits(cscBackwardLayers[1], event);  // ME1/2 and 1/3
0165   SegmentContainer CSClist0B = muonMeasurements->recHits(cscBackwardLayers[0], event);  // ME11
0166 
0167   BoolContainer usedCSClist4B(CSClist4B.size(), false);
0168   BoolContainer usedCSClist3B(CSClist3B.size(), false);
0169   BoolContainer usedCSClist2B(CSClist2B.size(), false);
0170   BoolContainer usedCSClist1B(CSClist1B.size(), false);
0171   BoolContainer usedCSClist0B(CSClist0B.size(), false);
0172 
0173   // 1b.2 Global z > 0
0174   std::vector<const DetLayer*> cscForwardLayers = muonLayers->forwardCSCLayers();
0175   SegmentContainer CSClist4F = muonMeasurements->recHits(cscForwardLayers[4], event);
0176   SegmentContainer CSClist3F = muonMeasurements->recHits(cscForwardLayers[3], event);
0177   SegmentContainer CSClist2F = muonMeasurements->recHits(cscForwardLayers[2], event);
0178   SegmentContainer CSClist1F = muonMeasurements->recHits(cscForwardLayers[1], event);
0179   SegmentContainer CSClist0F = muonMeasurements->recHits(cscForwardLayers[0], event);
0180 
0181   BoolContainer usedCSClist4F(CSClist4F.size(), false);
0182   BoolContainer usedCSClist3F(CSClist3F.size(), false);
0183   BoolContainer usedCSClist2F(CSClist2F.size(), false);
0184   BoolContainer usedCSClist1F(CSClist1F.size(), false);
0185   BoolContainer usedCSClist0F(CSClist0F.size(), false);
0186 
0187   if (debug) {
0188     std::cout << "*** Number of CSC segments is "
0189               << CSClist4F.size() + CSClist3F.size() + CSClist2F.size() + CSClist1F.size() + CSClist0F.size() +
0190                      CSClist4B.size() + CSClist3B.size() + CSClist2B.size() + CSClist1B.size() + CSClist0B.size()
0191               << std::endl;
0192     std::cout << "In ME11: " << CSClist0F.size() + CSClist0B.size() << std::endl;
0193     std::cout << "In ME12: " << CSClist1F.size() + CSClist1B.size() << std::endl;
0194     std::cout << "In ME2 : " << CSClist2F.size() + CSClist2B.size() << std::endl;
0195     std::cout << "In ME3 : " << CSClist3F.size() + CSClist3B.size() << std::endl;
0196     std::cout << "In ME4 : " << CSClist4F.size() + CSClist4B.size() << std::endl;
0197   }
0198 
0199   /* ******************************************************************************************************************
0200    * Form seeds in barrel region
0201    *
0202    * Proceed from inside -> out
0203    * ******************************************************************************************************************/
0204 
0205   // Loop over all possible MB1 segment to form seeds:
0206   int index = -1;
0207   for (SegmentContainer::iterator it = DTlist1.begin(); it != DTlist1.end(); ++it) {
0208     index++;
0209 
0210     if (usedDTlist1[index] == true)
0211       continue;
0212     if (int((*it)->recHits().size()) < minDTHitsPerSegment)
0213       continue;
0214     if ((*it)->dimension() != 4)
0215       continue;
0216 
0217     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0218     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0219 
0220     // Global position of starting point for protoTrack
0221     GlobalPoint gp = (*it)->globalPosition();
0222     float eta_temp = gp.eta();
0223     float phi_temp = gp.phi();
0224     bool showeringBefore = false;
0225     NShowerSeg = 0;
0226     if (IdentifyShowering(DTlist1, usedDTlist1, eta_temp, phi_temp, -1, NShowerSeg))
0227       showeringBefore = true;
0228     int NShowers = 0;
0229     if (showeringBefore) {
0230       //usedDTlist1[index] = true;
0231       NShowers++;
0232     }
0233 
0234     SegmentContainer protoTrack;
0235     protoTrack.push_back(*it);
0236 
0237     std::vector<int> layers;
0238     layers.push_back(-1);
0239 
0240     // Try adding segment from other stations
0241     if (foundMatchingSegment(
0242             3, protoTrack, DTlist2, usedDTlist2, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0243       layers.push_back(-2);
0244     if (showeringBefore)
0245       NShowers++;
0246     if (foundMatchingSegment(
0247             3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0248       layers.push_back(-3);
0249     if (showeringBefore)
0250       NShowers++;
0251     if (foundMatchingSegment(
0252             3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0253       layers.push_back(-4);
0254     if (showeringBefore)
0255       NShowers++;
0256 
0257     // Check if in overlap region
0258     if (eta_temp > 0.8) {
0259       if (foundMatchingSegment(
0260               2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0261         layers.push_back(1);
0262       if (showeringBefore)
0263         NShowers++;
0264       if (foundMatchingSegment(
0265               2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0266         layers.push_back(2);
0267       if (showeringBefore)
0268         NShowers++;
0269       if (foundMatchingSegment(
0270               2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0271         layers.push_back(3);
0272       if (showeringBefore)
0273         NShowers++;
0274     } else if (eta_temp < -0.8) {
0275       if (foundMatchingSegment(
0276               2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0277         layers.push_back(1);
0278       if (showeringBefore)
0279         NShowers++;
0280       if (foundMatchingSegment(
0281               2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0282         layers.push_back(2);
0283       if (showeringBefore)
0284         NShowers++;
0285       if (foundMatchingSegment(
0286               2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0287         layers.push_back(3);
0288       if (showeringBefore)
0289         NShowers++;
0290     }
0291 
0292     // adding showering information
0293     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0294       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0295         if (ShoweringLayers[i] > 0) {
0296           if (ShoweringLayers[i] <= layers[layers.size() - 1])
0297             continue;
0298           protoTrack.push_back(ShoweringSegments[i]);
0299           layers.push_back(ShoweringLayers[i]);
0300         }
0301         if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
0302           if (ShoweringLayers[i] >= layers[layers.size() - 1])
0303             continue;
0304           protoTrack.push_back(ShoweringSegments[i]);
0305           layers.push_back(ShoweringLayers[i]);
0306         }
0307       }
0308     }
0309     ShoweringSegments.clear();
0310     ShoweringLayers.clear();
0311 
0312     TrajectorySeed thisSeed;
0313 
0314     if (layers.size() < 2) {
0315       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0316     } else {
0317       if (layers[layers.size() - 1] > 0) {
0318         thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
0319       } else {
0320         thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
0321       }
0322     }
0323     // Add the seeds to master collection
0324     rawSeeds.push_back(thisSeed);
0325     etaOfSeed.push_back(eta_temp);
0326     phiOfSeed.push_back(phi_temp);
0327     nSegOnSeed.push_back(protoTrack.size());
0328 
0329     // Marked segment as used
0330     usedDTlist1[index] = true;
0331   }
0332 
0333   // Loop over all possible MB2 segment to form seeds:
0334   index = -1;
0335   for (SegmentContainer::iterator it = DTlist2.begin(); it != DTlist2.end(); ++it) {
0336     index++;
0337 
0338     if (usedDTlist2[index] == true)
0339       continue;
0340     if (int((*it)->recHits().size()) < minDTHitsPerSegment)
0341       continue;
0342     if ((*it)->dimension() != 4)
0343       continue;
0344 
0345     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0346     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0347 
0348     // Global position of starting point for protoTrack
0349     GlobalPoint gp = (*it)->globalPosition();
0350     float eta_temp = gp.eta();
0351     float phi_temp = gp.phi();
0352     bool showeringBefore = false;
0353     NShowerSeg = 0;
0354     if (IdentifyShowering(DTlist2, usedDTlist2, eta_temp, phi_temp, -2, NShowerSeg))
0355       showeringBefore = true;
0356     int NShowers = 0;
0357     if (showeringBefore) {
0358       // usedDTlist2[index] = true;
0359       NShowers++;
0360     }
0361 
0362     SegmentContainer protoTrack;
0363     protoTrack.push_back(*it);
0364 
0365     std::vector<int> layers;
0366     layers.push_back(-2);
0367 
0368     // Try adding segment from other stations
0369     if (foundMatchingSegment(
0370             3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0371       layers.push_back(-3);
0372     if (showeringBefore)
0373       NShowers++;
0374     if (foundMatchingSegment(
0375             3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0376       layers.push_back(-4);
0377     if (showeringBefore)
0378       NShowers++;
0379 
0380     // Check if in overlap region
0381     if (eta_temp > 0.8) {
0382       if (foundMatchingSegment(
0383               2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0384         layers.push_back(1);
0385       if (showeringBefore)
0386         NShowers++;
0387       if (foundMatchingSegment(
0388               2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0389         layers.push_back(2);
0390       if (showeringBefore)
0391         NShowers++;
0392       if (foundMatchingSegment(
0393               2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0394         layers.push_back(3);
0395       if (showeringBefore)
0396         NShowers++;
0397     } else if (eta_temp < -0.8) {
0398       if (foundMatchingSegment(
0399               2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0400         layers.push_back(1);
0401       if (showeringBefore)
0402         NShowers++;
0403       if (foundMatchingSegment(
0404               2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0405         layers.push_back(2);
0406       if (showeringBefore)
0407         NShowers++;
0408       if (foundMatchingSegment(
0409               2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0410         layers.push_back(3);
0411       if (showeringBefore)
0412         NShowers++;
0413     }
0414 
0415     // adding showering information
0416     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0417       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0418         if (ShoweringLayers[i] > 0) {
0419           if (ShoweringLayers[i] <= layers[layers.size() - 1])
0420             continue;
0421           protoTrack.push_back(ShoweringSegments[i]);
0422           layers.push_back(ShoweringLayers[i]);
0423         }
0424         if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
0425           if (ShoweringLayers[i] >= layers[layers.size() - 1])
0426             continue;
0427           protoTrack.push_back(ShoweringSegments[i]);
0428           layers.push_back(ShoweringLayers[i]);
0429         }
0430       }
0431     }
0432     ShoweringSegments.clear();
0433     ShoweringLayers.clear();
0434 
0435     TrajectorySeed thisSeed;
0436 
0437     if (layers.size() < 2) {
0438       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0439     } else {
0440       if (layers[layers.size() - 1] > 0) {
0441         thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
0442       } else {
0443         thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
0444       }
0445     }
0446     // Add the seeds to master collection
0447     rawSeeds.push_back(thisSeed);
0448     etaOfSeed.push_back(eta_temp);
0449     phiOfSeed.push_back(phi_temp);
0450     nSegOnSeed.push_back(protoTrack.size());
0451 
0452     // Marked segment as used
0453     usedDTlist2[index] = true;
0454   }
0455 
0456   // Loop over all possible MB3 segment to form seeds:
0457   index = -1;
0458   for (SegmentContainer::iterator it = DTlist3.begin(); it != DTlist3.end(); ++it) {
0459     index++;
0460 
0461     if (usedDTlist3[index] == true)
0462       continue;
0463     if (int((*it)->recHits().size()) < minDTHitsPerSegment)
0464       continue;
0465     if ((*it)->dimension() != 4)
0466       continue;
0467 
0468     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0469     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0470 
0471     // Global position of starting point for protoTrack
0472     GlobalPoint gp = (*it)->globalPosition();
0473     float eta_temp = gp.eta();
0474     float phi_temp = gp.phi();
0475     bool showeringBefore = false;
0476     NShowerSeg = 0;
0477     if (IdentifyShowering(DTlist3, usedDTlist3, eta_temp, phi_temp, -3, NShowerSeg))
0478       showeringBefore = true;
0479     int NShowers = 0;
0480     if (showeringBefore) {
0481       // usedDTlist3[index] = true;
0482       NShowers++;
0483     }
0484 
0485     SegmentContainer protoTrack;
0486     protoTrack.push_back(*it);
0487 
0488     std::vector<int> layers;
0489     layers.push_back(-3);
0490 
0491     // Try adding segment from other stations
0492     if (foundMatchingSegment(
0493             3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0494       layers.push_back(-4);
0495     if (showeringBefore)
0496       NShowers++;
0497 
0498     // Check if in overlap region
0499     if (eta_temp > 0.8) {
0500       if (foundMatchingSegment(
0501               2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0502         layers.push_back(1);
0503       if (showeringBefore)
0504         NShowers++;
0505       if (foundMatchingSegment(
0506               2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0507         layers.push_back(2);
0508       if (showeringBefore)
0509         NShowers++;
0510       if (foundMatchingSegment(
0511               2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0512         layers.push_back(3);
0513       if (showeringBefore)
0514         NShowers++;
0515     } else if (eta_temp < -0.8) {
0516       if (foundMatchingSegment(
0517               2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0518         layers.push_back(1);
0519       if (showeringBefore)
0520         NShowers++;
0521       if (foundMatchingSegment(
0522               2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0523         layers.push_back(2);
0524       if (showeringBefore)
0525         NShowers++;
0526       if (foundMatchingSegment(
0527               2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0528         layers.push_back(3);
0529       if (showeringBefore)
0530         NShowers++;
0531     }
0532 
0533     // adding showering information
0534     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0535       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0536         if (ShoweringLayers[i] > 0) {
0537           if (ShoweringLayers[i] <= layers[layers.size() - 1])
0538             continue;
0539           protoTrack.push_back(ShoweringSegments[i]);
0540           layers.push_back(ShoweringLayers[i]);
0541         }
0542         if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
0543           if (ShoweringLayers[i] >= layers[layers.size() - 1])
0544             continue;
0545           protoTrack.push_back(ShoweringSegments[i]);
0546           layers.push_back(ShoweringLayers[i]);
0547         }
0548       }
0549     }
0550     ShoweringSegments.clear();
0551     ShoweringLayers.clear();
0552 
0553     TrajectorySeed thisSeed;
0554     if (layers.size() < 2) {
0555       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0556     } else {
0557       if (layers[layers.size() - 1] > 0) {
0558         thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
0559       } else {
0560         thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
0561       }
0562     }
0563 
0564     // Add the seeds to master collection
0565     rawSeeds.push_back(thisSeed);
0566     etaOfSeed.push_back(eta_temp);
0567     phiOfSeed.push_back(phi_temp);
0568     nSegOnSeed.push_back(protoTrack.size());
0569 
0570     // Marked segment as used
0571     usedDTlist3[index] = true;
0572   }
0573 
0574   /* *********************************************************************************************************************
0575    * Form seeds from backward endcap
0576    *
0577    * Proceed from inside -> out
0578    * *********************************************************************************************************************/
0579 
0580   // Loop over all possible ME11 segment to form seeds:
0581   index = -1;
0582 
0583   for (SegmentContainer::iterator it = CSClist0B.begin(); it != CSClist0B.end(); ++it) {
0584     index++;
0585 
0586     if (usedCSClist0B[index] == true)
0587       continue;
0588     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0589       continue;
0590 
0591     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0592     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0593 
0594     // Global position of starting point for protoTrack
0595     GlobalPoint gp = (*it)->globalPosition();
0596     float eta_temp = gp.eta();
0597     float phi_temp = gp.phi();
0598     bool showeringBefore = false;
0599     NShowerSeg = 0;
0600     if (IdentifyShowering(CSClist0B, usedCSClist0B, eta_temp, phi_temp, 0, NShowerSeg))
0601       showeringBefore = true;
0602     int NShowers = 0;
0603     if (showeringBefore) {
0604       // usedCSClist0B[index] = true;
0605       NShowers++;
0606     }
0607 
0608     SegmentContainer protoTrack;
0609     protoTrack.push_back(*it);
0610 
0611     std::vector<int> layers;
0612     layers.push_back(0);
0613 
0614     // Try adding segment from other station
0615     if (foundMatchingSegment(
0616             1, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0617       layers.push_back(1);
0618     if (showeringBefore)
0619       NShowers++;
0620     if (foundMatchingSegment(
0621             1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0622       layers.push_back(2);
0623     if (showeringBefore)
0624       NShowers++;
0625     if (foundMatchingSegment(
0626             1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0627       layers.push_back(3);
0628     if (showeringBefore)
0629       NShowers++;
0630     if (foundMatchingSegment(
0631             1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0632       layers.push_back(4);
0633     if (showeringBefore)
0634       NShowers++;
0635 
0636     // adding showering information
0637     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0638       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0639         if (ShoweringLayers[i] <= layers[layers.size() - 1])
0640           continue;
0641         protoTrack.push_back(ShoweringSegments[i]);
0642         layers.push_back(ShoweringLayers[i]);
0643       }
0644     }
0645     ShoweringSegments.clear();
0646     ShoweringLayers.clear();
0647 
0648     TrajectorySeed thisSeed;
0649     if (layers.size() < 2) {
0650       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0651     } else {
0652       if (fabs(gp.eta()) > 1.7) {
0653         thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
0654       } else {
0655         thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
0656       }
0657     }
0658 
0659     // Add the seeds to master collection
0660     rawSeeds.push_back(thisSeed);
0661     etaOfSeed.push_back(eta_temp);
0662     phiOfSeed.push_back(phi_temp);
0663     nSegOnSeed.push_back(protoTrack.size());
0664 
0665     // mark this segment as used
0666     usedCSClist0B[index] = true;
0667   }
0668 
0669   // Loop over all possible ME1/2 ME1/3 segment to form seeds:
0670   index = -1;
0671   for (SegmentContainer::iterator it = CSClist1B.begin(); it != CSClist1B.end(); ++it) {
0672     index++;
0673 
0674     if (usedCSClist1B[index] == true)
0675       continue;
0676     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0677       continue;
0678 
0679     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0680     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0681 
0682     // Global position of starting point for protoTrack
0683     GlobalPoint gp = (*it)->globalPosition();
0684     float eta_temp = gp.eta();
0685     float phi_temp = gp.phi();
0686     bool showeringBefore = false;
0687     NShowerSeg = 0;
0688     if (IdentifyShowering(CSClist1B, usedCSClist1B, eta_temp, phi_temp, 1, NShowerSeg))
0689       showeringBefore = true;
0690     int NShowers = 0;
0691     if (showeringBefore) {
0692       //    usedCSClist1B[index] = true;
0693       NShowers++;
0694     }
0695 
0696     SegmentContainer protoTrack;
0697     protoTrack.push_back(*it);
0698 
0699     std::vector<int> layers;
0700     layers.push_back(1);
0701 
0702     // Try adding segment from other stations
0703     if (foundMatchingSegment(
0704             1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0705       layers.push_back(2);
0706     if (showeringBefore)
0707       NShowers++;
0708     if (foundMatchingSegment(
0709             1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0710       layers.push_back(3);
0711     if (showeringBefore)
0712       NShowers++;
0713     if (foundMatchingSegment(
0714             1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0715       layers.push_back(4);
0716     if (showeringBefore)
0717       NShowers++;
0718 
0719     // adding showering information
0720     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0721       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0722         if (ShoweringLayers[i] <= layers[layers.size() - 1])
0723           continue;
0724         protoTrack.push_back(ShoweringSegments[i]);
0725         layers.push_back(ShoweringLayers[i]);
0726       }
0727     }
0728     ShoweringSegments.clear();
0729     ShoweringLayers.clear();
0730 
0731     TrajectorySeed thisSeed;
0732     if (layers.size() < 2) {
0733       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0734     } else {
0735       thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
0736     }
0737     // Add the seeds to master collection
0738     rawSeeds.push_back(thisSeed);
0739     etaOfSeed.push_back(eta_temp);
0740     phiOfSeed.push_back(phi_temp);
0741     nSegOnSeed.push_back(protoTrack.size());
0742 
0743     // mark this segment as used
0744     usedCSClist1B[index] = true;
0745   }
0746 
0747   // Loop over all possible ME2 segment to form seeds:
0748   index = -1;
0749   for (SegmentContainer::iterator it = CSClist2B.begin(); it != CSClist2B.end(); ++it) {
0750     index++;
0751 
0752     if (usedCSClist2B[index] == true)
0753       continue;
0754     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0755       continue;
0756 
0757     double dof = static_cast<double>((*it)->degreesOfFreedom());
0758     if (((*it)->chi2() / dof) > 20000.0)
0759       continue;
0760 
0761     // Global position of starting point for protoTrack
0762     GlobalPoint gp = (*it)->globalPosition();
0763     float eta_temp = gp.eta();
0764     float phi_temp = gp.phi();
0765     bool showeringBefore = false;
0766     NShowerSeg = 0;
0767     if (IdentifyShowering(CSClist2B, usedCSClist2B, eta_temp, phi_temp, 2, NShowerSeg))
0768       showeringBefore = true;
0769     int NShowers = 0;
0770     if (showeringBefore) {
0771       // usedCSClist2B[index] = true;
0772       NShowers++;
0773     }
0774 
0775     SegmentContainer protoTrack;
0776     protoTrack.push_back(*it);
0777 
0778     std::vector<int> layers;
0779     layers.push_back(2);
0780 
0781     // Try adding segment from other stations
0782     if (foundMatchingSegment(
0783             1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0784       layers.push_back(3);
0785     if (showeringBefore)
0786       NShowers++;
0787     if (foundMatchingSegment(
0788             1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0789       layers.push_back(4);
0790     if (showeringBefore)
0791       NShowers++;
0792 
0793     // adding showering information
0794     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0795       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0796         if (ShoweringLayers[i] <= layers[layers.size() - 1])
0797           continue;
0798         protoTrack.push_back(ShoweringSegments[i]);
0799         layers.push_back(ShoweringLayers[i]);
0800       }
0801     }
0802     ShoweringSegments.clear();
0803     ShoweringLayers.clear();
0804 
0805     TrajectorySeed thisSeed;
0806     if (layers.size() < 2) {
0807       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0808     } else {
0809       thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
0810     }
0811 
0812     // Add the seeds to master collection
0813     rawSeeds.push_back(thisSeed);
0814     etaOfSeed.push_back(eta_temp);
0815     phiOfSeed.push_back(phi_temp);
0816     nSegOnSeed.push_back(protoTrack.size());
0817     // mark this segment as used
0818     usedCSClist2B[index] = true;
0819   }
0820 
0821   // Loop over all possible ME3 segment to form seeds:
0822   index = -1;
0823   for (SegmentContainer::iterator it = CSClist3B.begin(); it != CSClist3B.end(); ++it) {
0824     index++;
0825 
0826     if (usedCSClist3B[index] == true)
0827       continue;
0828     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0829       continue;
0830 
0831     double dof = static_cast<double>((*it)->degreesOfFreedom());
0832     if (((*it)->chi2() / dof) > 20000.0)
0833       continue;
0834 
0835     // Global position of starting point for protoTrack
0836     GlobalPoint gp = (*it)->globalPosition();
0837     float eta_temp = gp.eta();
0838     float phi_temp = gp.phi();
0839     bool showeringBefore = false;
0840     NShowerSeg = 0;
0841     if (IdentifyShowering(CSClist3B, usedCSClist3B, eta_temp, phi_temp, 3, NShowerSeg))
0842       showeringBefore = true;
0843     int NShowers = 0;
0844     if (showeringBefore) {
0845       //    usedCSClist3B[index] = true;
0846       NShowers++;
0847     }
0848 
0849     SegmentContainer protoTrack;
0850     protoTrack.push_back(*it);
0851 
0852     std::vector<int> layers;
0853     layers.push_back(2);
0854 
0855     // Try adding segment from other stations
0856     if (foundMatchingSegment(
0857             1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0858       layers.push_back(4);
0859     if (showeringBefore)
0860       NShowers++;
0861 
0862     // adding showering information
0863     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0864       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0865         if (ShoweringLayers[i] <= layers[layers.size() - 1])
0866           continue;
0867         protoTrack.push_back(ShoweringSegments[i]);
0868         layers.push_back(ShoweringLayers[i]);
0869       }
0870     }
0871     ShoweringSegments.clear();
0872     ShoweringLayers.clear();
0873 
0874     // mark this segment as used
0875     usedCSClist3B[index] = true;
0876 
0877     if (layers.size() < 2)
0878       continue;
0879     TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
0880 
0881     // Add the seeds to master collection
0882     rawSeeds.push_back(thisSeed);
0883     etaOfSeed.push_back(eta_temp);
0884     phiOfSeed.push_back(phi_temp);
0885     nSegOnSeed.push_back(protoTrack.size());
0886   }
0887 
0888   /* *****************************************************************************************************************
0889    * Form seeds from forward endcap
0890    *
0891    * Proceed from inside -> out
0892    * *****************************************************************************************************************/
0893 
0894   // Loop over all possible ME11 segment to form seeds:
0895   index = -1;
0896   for (SegmentContainer::iterator it = CSClist0F.begin(); it != CSClist0F.end(); ++it) {
0897     index++;
0898 
0899     if (usedCSClist0F[index] == true)
0900       continue;
0901     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0902       continue;
0903 
0904     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0905     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0906 
0907     // Global position of starting point for protoTrack
0908     GlobalPoint gp = (*it)->globalPosition();
0909     float eta_temp = gp.eta();
0910     float phi_temp = gp.phi();
0911     bool showeringBefore = false;
0912     NShowerSeg = 0;
0913     if (IdentifyShowering(CSClist0F, usedCSClist0F, eta_temp, phi_temp, 0, NShowerSeg))
0914       showeringBefore = true;
0915     int NShowers = 0;
0916     if (showeringBefore) {
0917       // usedCSClist0F[index] = true;
0918       NShowers++;
0919     }
0920 
0921     SegmentContainer protoTrack;
0922     protoTrack.push_back(*it);
0923 
0924     std::vector<int> layers;
0925     layers.push_back(0);
0926 
0927     // Try adding segment from other station
0928     if (foundMatchingSegment(
0929             1, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0930       layers.push_back(1);
0931     if (showeringBefore)
0932       NShowers++;
0933     if (foundMatchingSegment(
0934             1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0935       layers.push_back(2);
0936     if (showeringBefore)
0937       NShowers++;
0938     if (foundMatchingSegment(
0939             1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0940       layers.push_back(3);
0941     if (showeringBefore)
0942       NShowers++;
0943     if (foundMatchingSegment(
0944             1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
0945       layers.push_back(4);
0946     if (showeringBefore)
0947       NShowers++;
0948 
0949     // adding showering information
0950     if (layers.size() < 2 && !ShoweringSegments.empty()) {
0951       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
0952         if (ShoweringLayers[i] <= layers[layers.size() - 1])
0953           continue;
0954         protoTrack.push_back(ShoweringSegments[i]);
0955         layers.push_back(ShoweringLayers[i]);
0956       }
0957     }
0958     ShoweringSegments.clear();
0959     ShoweringLayers.clear();
0960 
0961     TrajectorySeed thisSeed;
0962     if (layers.size() < 2) {
0963       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
0964     } else {
0965       if (fabs(gp.eta()) > 1.7) {
0966         thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
0967       } else {
0968         thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
0969       }
0970     }
0971     // Add the seeds to master collection
0972     rawSeeds.push_back(thisSeed);
0973     etaOfSeed.push_back(eta_temp);
0974     phiOfSeed.push_back(phi_temp);
0975     nSegOnSeed.push_back(protoTrack.size());
0976 
0977     // mark this segment as used
0978     usedCSClist0F[index] = true;
0979   }
0980 
0981   // Loop over all possible ME1/2 ME1/3 segment to form seeds:
0982   index = -1;
0983   for (SegmentContainer::iterator it = CSClist1F.begin(); it != CSClist1F.end(); ++it) {
0984     index++;
0985 
0986     if (usedCSClist1F[index] == true)
0987       continue;
0988     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
0989       continue;
0990 
0991     //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
0992     //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
0993 
0994     // Global position of starting point for protoTrack
0995     GlobalPoint gp = (*it)->globalPosition();
0996     float eta_temp = gp.eta();
0997     float phi_temp = gp.phi();
0998     bool showeringBefore = false;
0999     NShowerSeg = 0;
1000     if (IdentifyShowering(CSClist1F, usedCSClist1F, eta_temp, phi_temp, 1, NShowerSeg))
1001       showeringBefore = true;
1002     int NShowers = 0;
1003     if (showeringBefore) {
1004       //   usedCSClist1F[index] = true;
1005       NShowers++;
1006     }
1007 
1008     SegmentContainer protoTrack;
1009     protoTrack.push_back(*it);
1010 
1011     std::vector<int> layers;
1012     layers.push_back(1);
1013 
1014     // Try adding segment from other stations
1015     if (foundMatchingSegment(
1016             1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1017       layers.push_back(2);
1018     if (showeringBefore)
1019       NShowers++;
1020     if (foundMatchingSegment(
1021             1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1022       layers.push_back(3);
1023     if (showeringBefore)
1024       NShowers++;
1025     if (foundMatchingSegment(
1026             1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1027       layers.push_back(4);
1028     if (showeringBefore)
1029       NShowers++;
1030 
1031     // adding showering information
1032     if (layers.size() < 2 && !ShoweringSegments.empty()) {
1033       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1034         if (ShoweringLayers[i] <= layers[layers.size() - 1])
1035           continue;
1036         protoTrack.push_back(ShoweringSegments[i]);
1037         layers.push_back(ShoweringLayers[i]);
1038       }
1039     }
1040     ShoweringSegments.clear();
1041     ShoweringLayers.clear();
1042 
1043     TrajectorySeed thisSeed;
1044     if (layers.size() < 2) {
1045       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
1046     } else {
1047       thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1048     }
1049 
1050     // Add the seeds to master collection
1051     rawSeeds.push_back(thisSeed);
1052     etaOfSeed.push_back(eta_temp);
1053     phiOfSeed.push_back(phi_temp);
1054     nSegOnSeed.push_back(protoTrack.size());
1055 
1056     // mark this segment as used
1057     usedCSClist1F[index] = true;
1058   }
1059 
1060   // Loop over all possible ME2 segment to form seeds:
1061   index = -1;
1062   for (SegmentContainer::iterator it = CSClist2F.begin(); it != CSClist2F.end(); ++it) {
1063     index++;
1064 
1065     if (usedCSClist2F[index] == true)
1066       continue;
1067     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
1068       continue;
1069 
1070     double dof = static_cast<double>((*it)->degreesOfFreedom());
1071     if (((*it)->chi2() / dof) > 20000.0)
1072       continue;
1073 
1074     // Global position of starting point for protoTrack
1075     GlobalPoint gp = (*it)->globalPosition();
1076     float eta_temp = gp.eta();
1077     float phi_temp = gp.phi();
1078     bool showeringBefore = false;
1079     NShowerSeg = 0;
1080     if (IdentifyShowering(CSClist2F, usedCSClist2F, eta_temp, phi_temp, 2, NShowerSeg))
1081       showeringBefore = true;
1082     int NShowers = 0;
1083     if (showeringBefore) {
1084       //    usedCSClist2F[index] = true;
1085       NShowers++;
1086     }
1087 
1088     SegmentContainer protoTrack;
1089     protoTrack.push_back(*it);
1090 
1091     std::vector<int> layers;
1092     layers.push_back(2);
1093 
1094     // Try adding segment from other stations
1095     if (foundMatchingSegment(
1096             1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1097       layers.push_back(3);
1098     if (showeringBefore)
1099       NShowers++;
1100     if (foundMatchingSegment(
1101             1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1102       layers.push_back(4);
1103     if (showeringBefore)
1104       NShowers++;
1105 
1106     // adding showering information
1107     if (layers.size() < 2 && !ShoweringSegments.empty()) {
1108       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1109         if (ShoweringLayers[i] <= layers[layers.size() - 1])
1110           continue;
1111         protoTrack.push_back(ShoweringSegments[i]);
1112         layers.push_back(ShoweringLayers[i]);
1113       }
1114     }
1115     ShoweringSegments.clear();
1116     ShoweringLayers.clear();
1117 
1118     TrajectorySeed thisSeed;
1119     if (layers.size() < 2) {
1120       thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
1121     } else {
1122       thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1123     }
1124 
1125     // Add the seeds to master collection
1126     rawSeeds.push_back(thisSeed);
1127     etaOfSeed.push_back(eta_temp);
1128     phiOfSeed.push_back(phi_temp);
1129     nSegOnSeed.push_back(protoTrack.size());
1130 
1131     // mark this segment as used
1132     usedCSClist2F[index] = true;
1133   }
1134 
1135   // Loop over all possible ME3 segment to form seeds:
1136   index = -1;
1137   for (SegmentContainer::iterator it = CSClist3F.begin(); it != CSClist3F.end(); ++it) {
1138     index++;
1139 
1140     if (usedCSClist3F[index] == true)
1141       continue;
1142     if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
1143       continue;
1144 
1145     double dof = static_cast<double>((*it)->degreesOfFreedom());
1146     if (((*it)->chi2() / dof) > 20000.0)
1147       continue;
1148 
1149     // Global position of starting point for protoTrack
1150     GlobalPoint gp = (*it)->globalPosition();
1151     float eta_temp = gp.eta();
1152     float phi_temp = gp.phi();
1153     bool showeringBefore = false;
1154     NShowerSeg = 0;
1155     if (IdentifyShowering(CSClist3F, usedCSClist3F, eta_temp, phi_temp, 3, NShowerSeg))
1156       showeringBefore = true;
1157     int NShowers = 0;
1158     if (showeringBefore) {
1159       //   usedCSClist3F[index] = true;
1160       NShowers++;
1161     }
1162 
1163     SegmentContainer protoTrack;
1164     protoTrack.push_back(*it);
1165 
1166     std::vector<int> layers;
1167     layers.push_back(2);
1168 
1169     // Try adding segment from other stations
1170     if (foundMatchingSegment(
1171             1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1172       layers.push_back(4);
1173     if (showeringBefore)
1174       NShowers++;
1175 
1176     // adding showering information
1177     if (layers.size() < 2 && !ShoweringSegments.empty()) {
1178       for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1179         if (ShoweringLayers[i] <= layers[layers.size() - 1])
1180           continue;
1181         protoTrack.push_back(ShoweringSegments[i]);
1182         layers.push_back(ShoweringLayers[i]);
1183       }
1184     }
1185     ShoweringSegments.clear();
1186     ShoweringLayers.clear();
1187 
1188     // mark this segment as used
1189     usedCSClist3F[index] = true;
1190 
1191     if (layers.size() < 2)
1192       continue;
1193 
1194     TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1195 
1196     // Add the seeds to master collection
1197     rawSeeds.push_back(thisSeed);
1198     etaOfSeed.push_back(eta_temp);
1199     phiOfSeed.push_back(phi_temp);
1200     nSegOnSeed.push_back(protoTrack.size());
1201   }
1202 
1203   /* *********************************************************************************************************************
1204    * Clean up raw seed collection and pass to master collection
1205    * *********************************************************************************************************************/
1206 
1207   if (debug)
1208     std::cout << "*** CLEAN UP " << std::endl;
1209   if (debug)
1210     std::cout << "Number of seeds BEFORE " << rawSeeds.size() << std::endl;
1211 
1212   int goodSeeds = 0;
1213 
1214   theSeeds = muonSeedClean_->seedCleaner(eventSetup, rawSeeds);
1215   goodSeeds = theSeeds.size();
1216 
1217   //std::cout << " === Before Cleaning : " << rawSeeds.size() <<std::endl;
1218   //std::cout << " => AFTER :" << goodSeeds << std::endl;
1219 
1220   if (debug)
1221     std::cout << "Number of seeds AFTER " << goodSeeds << std::endl;
1222 
1223   return goodSeeds;
1224 }
1225 
1226 /* *********************************************************************************************************************
1227  * Try to match segment from different station (layer)
1228  *
1229  * Note type = 1 --> endcap
1230  *           = 2 --> overlap
1231  *           = 3 --> barrel
1232  * *********************************************************************************************************************/
1233 
1234 ///                                                   segment for seeding         , segments collection
1235 bool MuonSeedBuilder::foundMatchingSegment(int type,
1236                                            SegmentContainer& protoTrack,
1237                                            SegmentContainer& segs,
1238                                            BoolContainer& usedSeg,
1239                                            float& eta_last,
1240                                            float& phi_last,
1241                                            int& lastLayer,
1242                                            bool& showeringBefore) {
1243   bool ok = false;
1244   int scanlayer = (lastLayer < 0) ? (lastLayer - 1) : (lastLayer + 1);
1245 
1246   if (IdentifyShowering(segs, usedSeg, eta_last, phi_last, scanlayer, NShowerSeg)) {
1247     showeringBefore = true;
1248     return ok;
1249   }
1250 
1251   // Setup the searching cone-size
1252   // searching cone-size should changes with bending power
1253   double maxdEta;
1254   double maxdPhi;
1255   if (type == 1) {
1256     // CSC
1257     maxdEta = maxDeltaEtaCSC;
1258     if (lastLayer == 0 || lastLayer == 1) {
1259       if (fabs(eta_last) < 2.1) {
1260         maxdPhi = maxDeltaPhiCSC;
1261       } else {
1262         maxdPhi = 0.06;
1263       }
1264     } else if (lastLayer == 2) {
1265       maxdPhi = 0.5 * maxDeltaPhiCSC;
1266     } else {
1267       maxdPhi = 0.2 * maxDeltaPhiCSC;
1268     }
1269 
1270   } else if (type == 2) {
1271     // Overlap
1272     maxdEta = maxDeltaEtaOverlap;
1273     if (lastLayer == -1) {
1274       maxdPhi = maxDeltaPhiDT;
1275     } else {
1276       maxdPhi = maxDeltaPhiOverlap;
1277     }
1278 
1279   } else {
1280     // DT
1281     maxdEta = maxDeltaEtaDT;
1282     if (lastLayer == -1) {
1283       maxdPhi = maxDeltaPhiDT;
1284     } else if (lastLayer == -2) {
1285       maxdPhi = 0.8 * maxDeltaPhiDT;
1286     } else {
1287       maxdPhi = 0.4 * maxDeltaPhiDT;
1288     }
1289   }
1290 
1291   // if previous layer showers, limite the maxdPhi < 0.06
1292   if (showeringBefore && maxdPhi > 0.03)
1293     maxdPhi = 0.03;
1294   // reset the showering flag
1295   showeringBefore = false;
1296 
1297   // global phi/eta from previous segment
1298   float eta_temp = eta_last;
1299   float phi_temp = phi_last;
1300 
1301   // Index counter to keep track of element used in segs
1302   int index = -1;
1303   int best_match = index;
1304   float best_R = sqrt((maxdEta * maxdEta) + (maxdPhi * maxdPhi));
1305   float best_chi2 = 200;
1306   int best_dimension = 2;
1307   int best_nhits = minDTHitsPerSegment;
1308   if (type == 1)
1309     best_nhits = minCSCHitsPerSegment;
1310   // Loop over segments in other station (layer) and find candidate match
1311   for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1312     index++;
1313 
1314     // Not to get confused:  eta_last is from the previous layer.
1315     // This is only to find the best set of segments by comparing at the distance layer-by-layer
1316     GlobalPoint gp2 = (*it)->globalPosition();
1317     double dh = fabs(gp2.eta() - eta_temp);
1318     double df = fabs(gp2.phi() - phi_temp);
1319     double dR = sqrt((dh * dh) + (df * df));
1320 
1321     // dEta and dPhi should be within certain range
1322     bool case1 = (dh < maxdEta && df < maxdPhi) ? true : false;
1323     // for DT station 4 ; CSCSegment is always 4D
1324     bool case2 = (((*it)->dimension() != 4) && (dh < 0.5) && (df < maxdPhi)) ? true : false;
1325     if (!case1 && !case2)
1326       continue;
1327 
1328     int NRechits = muonSeedClean_->NRecHitsFromSegment(&*(*it));
1329 
1330     if (NRechits < best_nhits)
1331       continue;
1332     best_nhits = NRechits;
1333 
1334     // reject 2D segments if 4D segments are available
1335     if ((*it)->dimension() < best_dimension)
1336       continue;
1337     best_dimension = (*it)->dimension();
1338 
1339     // pick the segment with best chi2/dof within a fixed cone size
1340     if (dR > best_R)
1341       continue;
1342 
1343     // select smaller chi2/dof
1344     double dof = static_cast<double>((*it)->degreesOfFreedom());
1345     /// reject possible edge segments
1346     if ((*it)->chi2() / dof < 0.001 && NRechits < 6 && type == 1)
1347       continue;
1348     if ((*it)->chi2() / dof > best_chi2)
1349       continue;
1350     best_chi2 = (*it)->chi2() / dof;
1351     best_match = index;
1352     // propagate the eta and phi to next layer
1353     if ((*it)->dimension() != 4) {
1354       // Self assignment, is this a bug??
1355       // should this have been (phi/eta)_last = (phi/eta)_temp to make it reset?
1356       //phi_last = phi_last;
1357       //eta_last = eta_last;
1358     } else {
1359       phi_last = gp2.phi();
1360       eta_last = gp2.eta();
1361     }
1362   }
1363 
1364   if (best_match < 0)
1365     return ok;
1366 
1367   // Add best matching segment to protoTrack:
1368   index = -1;
1369   for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1370     index++;
1371     if (index != best_match)
1372       continue;
1373     protoTrack.push_back(*it);
1374     usedSeg[best_match] = true;
1375     ok = true;
1376   }
1377   return ok;
1378 }
1379 
1380 bool MuonSeedBuilder::IdentifyShowering(SegmentContainer& segs,
1381                                         BoolContainer& usedSeg,
1382                                         float& eta_last,
1383                                         float& phi_last,
1384                                         int layer,
1385                                         int& NShoweringSegments) {
1386   bool showering = false;
1387 
1388   int nSeg = 0;
1389   int nRhits = 0;
1390   double nChi2 = 9999.;
1391   int theOrigin = -1;
1392   std::vector<int> badtag;
1393   int index = -1;
1394   double aveEta = 0.0;
1395   for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1396     index++;
1397     GlobalPoint gp = (*it)->globalPosition();
1398     double dh = gp.eta() - eta_last;
1399     double df = gp.phi() - phi_last;
1400     double dR = sqrt((dh * dh) + (df * df));
1401 
1402     double dof = static_cast<double>((*it)->degreesOfFreedom());
1403     double nX2 = (*it)->chi2() / dof;
1404 
1405     bool isDT = false;
1406     DetId geoId = (*it)->geographicalId();
1407     if (geoId.subdetId() == MuonSubdetId::DT)
1408       isDT = true;
1409 
1410     if (dR < 0.3) {
1411       nSeg++;
1412       badtag.push_back(index);
1413       aveEta += fabs(gp.eta());
1414       // pick up the best segment from showering chamber
1415       int rh = muonSeedClean_->NRecHitsFromSegment(&*(*it));
1416       if (rh < 6 && !isDT)
1417         continue;
1418       if (rh < 12 && isDT)
1419         continue;
1420       if (rh > nRhits) {
1421         nRhits = rh;
1422         if (nX2 > nChi2)
1423           continue;
1424         if (layer != 0 && layer != 1 && layer != -1) {
1425           theOrigin = index;
1426         }
1427       }
1428     }
1429   }
1430   aveEta = aveEta / static_cast<double>(nSeg);
1431   bool isME11A = (aveEta >= 2.1 && layer == 0) ? true : false;
1432   bool isME12 = (aveEta > 1.2 && aveEta <= 1.65 && layer == 1) ? true : false;
1433   bool isME11 = (aveEta > 1.65 && aveEta <= 2.1 && layer == 0) ? true : false;
1434   bool is1stLayer = (layer == -1 || layer == 0 || isME12 || isME11 || isME11A) ? true : false;
1435 
1436   NShoweringSegments += nSeg;
1437 
1438   if (nSeg > 3 && !isME11A)
1439     showering = true;
1440   if (nSeg > 6 && isME11A)
1441     showering = true;
1442 
1443   // if showering, flag all segments in order to skip this layer for pt estimation except 1st layer
1444   //std::cout<<" from Showering "<<std::endl;
1445   if (showering && !is1stLayer) {
1446     for (std::vector<int>::iterator it = badtag.begin(); it != badtag.end(); ++it) {
1447       usedSeg[*it] = true;
1448       if ((*it) != theOrigin)
1449         continue;
1450       ShoweringSegments.push_back(segs[*it]);
1451       ShoweringLayers.push_back(layer);
1452     }
1453   }
1454   return showering;
1455 }
1456 
1457 double MuonSeedBuilder::etaError(const GlobalPoint gp, double rErr) {
1458   double dHdTheta = 0.0;
1459   double dThetadR = 0.0;
1460   double etaErr = 1.0;
1461 
1462   if (gp.perp() != 0) {
1463     dHdTheta = (gp.mag() + gp.z()) / gp.perp();
1464     dThetadR = gp.z() / gp.perp2();
1465     etaErr = 0.25 * (dHdTheta * dThetadR) * (dHdTheta * dThetadR) * rErr;
1466   }
1467 
1468   return etaErr;
1469 }
[ 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!