CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonSeedBuilder.cc
Go to the documentation of this file.
1 
8 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedBuilder.h>
9 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCreator.h>
10 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCleaner.h>
11 
12 // Data Formats
13 #include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
14 #include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
15 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
16 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
17 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
18 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
19 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
20 #include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>
21 
22 // Geometry
23 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
24 #include <TrackingTools/DetLayers/interface/DetLayer.h>
25 #include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
26 #include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
27 #include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>
28 
29 // muon service
30 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
31 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
32 
33 // Framework
34 #include <FWCore/ParameterSet/interface/ParameterSet.h>
35 #include "FWCore/Framework/interface/Event.h"
36 #include <FWCore/Framework/interface/ESHandle.h>
37 #include <FWCore/MessageLogger/interface/MessageLogger.h>
38 #include <DataFormats/Common/interface/Handle.h>
39 
40 // C++
41 #include <vector>
42 #include <deque>
43 #include <utility>
44 
45 //typedef std::pair<double, TrajectorySeed> seedpr ;
46 //static bool ptDecreasing(const seedpr s1, const seedpr s2) { return ( s1.first > s2.first ); }
47 //static bool lengthSorting(const TrajectorySeed s1, const TrajectorySeed s2) { return ( s1.nHits() > s2.nHits() ); }
48 
49 /*
50  * Constructor
51  */
52 MuonSeedBuilder::MuonSeedBuilder(const edm::ParameterSet& pset,edm::ConsumesCollector& iC){
53 
54  // Local Debug flag
55  debug = pset.getParameter<bool>("DebugMuonSeed");
56 
57  // enable the DT chamber
58  enableDTMeasurement = pset.getParameter<bool>("EnableDTMeasurement");
59  theDTSegmentLabel = pset.getParameter<edm::InputTag>("DTSegmentLabel");
60 
61  // enable the CSC chamber
62  enableCSCMeasurement = pset.getParameter<bool>("EnableCSCMeasurement");
63  theCSCSegmentLabel = pset.getParameter<edm::InputTag>("CSCSegmentLabel");
64 
65  // Parameters for seed creation in endcap region
66  minCSCHitsPerSegment = pset.getParameter<int>("minCSCHitsPerSegment");
67  maxDeltaEtaCSC = pset.getParameter<double>("maxDeltaEtaCSC");
68  maxDeltaPhiCSC = pset.getParameter<double>("maxDeltaPhiCSC");
69 
70  // Parameters for seed creation in overlap region
71  maxDeltaEtaOverlap = pset.getParameter<double>("maxDeltaEtaOverlap");
72  maxDeltaPhiOverlap = pset.getParameter<double>("maxDeltaPhiOverlap");
73 
74  // Parameters for seed creation in barrel region
75  minDTHitsPerSegment = pset.getParameter<int>("minDTHitsPerSegment");
76  maxDeltaEtaDT = pset.getParameter<double>("maxDeltaEtaDT");
77  maxDeltaPhiDT = pset.getParameter<double>("maxDeltaPhiDT");
78 
79  // Parameters to suppress combinatorics (ghosts)
80  maxEtaResolutionDT = pset.getParameter<double>("maxEtaResolutionDT");
81  maxPhiResolutionDT = pset.getParameter<double>("maxPhiResolutionDT");
82  maxEtaResolutionCSC = pset.getParameter<double>("maxEtaResolutionCSC");
83  maxPhiResolutionCSC = pset.getParameter<double>("maxPhiResolutionCSC");
84 
85  // muon service
86  edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
87  theService = new MuonServiceProxy(serviceParameters);
88 
89  // Class for forming muon seeds
90  muonSeedCreate_ = new MuonSeedCreator( pset );
91  muonSeedClean_ = new MuonSeedCleaner( pset );
92 
93  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
94  muonMeasurements = new MuonDetLayerMeasurements(theDTSegmentLabel,theCSCSegmentLabel,edm::InputTag(),iC,
95  enableDTMeasurement,enableCSCMeasurement,false);
96 
97 
98 
99 }
100 
101 /*
102  * Destructor
103  */
104 MuonSeedBuilder::~MuonSeedBuilder(){
105 
106  delete muonSeedCreate_;
107  delete muonSeedClean_;
108  if (theService) delete theService;
109  if( muonMeasurements) delete muonMeasurements;
110 
111 }
112 
113 /*
114  * build
115  *
116  * Where the segments are sorted out to make a protoTrack (vector of matching segments in different
117  * stations/layers). The protoTrack is then passed to the seed creator to create CSC, DT or overlap seeds
118  *
119  */
120 //void MuonSeedBuilder::build( MuonDetLayerMeasurements muonMeasurements, TrajectorySeedCollection& theSeeds ) {
121 int MuonSeedBuilder::build( edm::Event& event, const edm::EventSetup& eventSetup, TrajectorySeedCollection& theSeeds ) {
122 
123 
124  // Pass the Magnetic Field to where the seed is create
125  muonSeedCreate_->setBField(BField);
126 
127  // Create temporary collection of seeds which will be cleaned up to remove combinatorics
128  std::vector<TrajectorySeed> rawSeeds;
129  std::vector<float> etaOfSeed;
130  std::vector<float> phiOfSeed;
131  std::vector<int> nSegOnSeed;
132 
133 
134  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
135  // MuonDetLayerMeasurements muonMeasurements(enableDTMeasurement,enableCSCMeasurement,false,
136  // theDTSegmentLabel.label(),theCSCSegmentLabel.label());
137 
138 
139  // 1) Get the various stations and store segments in containers for each station (layers)
140 
141  // 1a. get the DT segments by stations (layers):
142  std::vector<const DetLayer*> dtLayers = muonLayers->allDTLayers();
143 
144  SegmentContainer DTlist4 = muonMeasurements->recHits( dtLayers[3], event );
145  SegmentContainer DTlist3 = muonMeasurements->recHits( dtLayers[2], event );
146  SegmentContainer DTlist2 = muonMeasurements->recHits( dtLayers[1], event );
147  SegmentContainer DTlist1 = muonMeasurements->recHits( dtLayers[0], event );
148 
149  // Initialize flags that a given segment has been allocated to a seed
150  BoolContainer usedDTlist4(DTlist4.size(), false);
151  BoolContainer usedDTlist3(DTlist3.size(), false);
152  BoolContainer usedDTlist2(DTlist2.size(), false);
153  BoolContainer usedDTlist1(DTlist1.size(), false);
154 
155  if (debug) {
156  std::cout << "*** Number of DT segments is: " << DTlist4.size()+DTlist3.size()+DTlist2.size()+DTlist1.size() << std::endl;
157  std::cout << "In MB1: " << DTlist1.size() << std::endl;
158  std::cout << "In MB2: " << DTlist2.size() << std::endl;
159  std::cout << "In MB3: " << DTlist3.size() << std::endl;
160  std::cout << "In MB4: " << DTlist4.size() << std::endl;
161  }
162 
163  // 1b. get the CSC segments by stations (layers):
164  // 1b.1 Global z < 0
165  std::vector<const DetLayer*> cscBackwardLayers = muonLayers->backwardCSCLayers();
166  SegmentContainer CSClist4B = muonMeasurements->recHits( cscBackwardLayers[4], event );
167  SegmentContainer CSClist3B = muonMeasurements->recHits( cscBackwardLayers[3], event );
168  SegmentContainer CSClist2B = muonMeasurements->recHits( cscBackwardLayers[2], event );
169  SegmentContainer CSClist1B = muonMeasurements->recHits( cscBackwardLayers[1], event ); // ME1/2 and 1/3
170  SegmentContainer CSClist0B = muonMeasurements->recHits( cscBackwardLayers[0], event ); // ME11
171 
172  BoolContainer usedCSClist4B(CSClist4B.size(), false);
173  BoolContainer usedCSClist3B(CSClist3B.size(), false);
174  BoolContainer usedCSClist2B(CSClist2B.size(), false);
175  BoolContainer usedCSClist1B(CSClist1B.size(), false);
176  BoolContainer usedCSClist0B(CSClist0B.size(), false);
177 
178  // 1b.2 Global z > 0
179  std::vector<const DetLayer*> cscForwardLayers = muonLayers->forwardCSCLayers();
180  SegmentContainer CSClist4F = muonMeasurements->recHits( cscForwardLayers[4], event );
181  SegmentContainer CSClist3F = muonMeasurements->recHits( cscForwardLayers[3], event );
182  SegmentContainer CSClist2F = muonMeasurements->recHits( cscForwardLayers[2], event );
183  SegmentContainer CSClist1F = muonMeasurements->recHits( cscForwardLayers[1], event );
184  SegmentContainer CSClist0F = muonMeasurements->recHits( cscForwardLayers[0], event );
185 
186  BoolContainer usedCSClist4F(CSClist4F.size(), false);
187  BoolContainer usedCSClist3F(CSClist3F.size(), false);
188  BoolContainer usedCSClist2F(CSClist2F.size(), false);
189  BoolContainer usedCSClist1F(CSClist1F.size(), false);
190  BoolContainer usedCSClist0F(CSClist0F.size(), false);
191 
192  if (debug) {
193  std::cout << "*** Number of CSC segments is " << CSClist4F.size()+CSClist3F.size()+CSClist2F.size()+CSClist1F.size()+CSClist0F.size()+CSClist4B.size()+CSClist3B.size()+CSClist2B.size()+CSClist1B.size()+CSClist0B.size()<< std::endl;
194  std::cout << "In ME11: " << CSClist0F.size()+CSClist0B.size() << std::endl;
195  std::cout << "In ME12: " << CSClist1F.size()+CSClist1B.size() << std::endl;
196  std::cout << "In ME2 : " << CSClist2F.size()+CSClist2B.size() << std::endl;
197  std::cout << "In ME3 : " << CSClist3F.size()+CSClist3B.size() << std::endl;
198  std::cout << "In ME4 : " << CSClist4F.size()+CSClist4B.size() << std::endl;
199  }
200 
201 
202  /* ******************************************************************************************************************
203  * Form seeds in barrel region
204  *
205  * Proceed from inside -> out
206  * ******************************************************************************************************************/
207 
208 
209  // Loop over all possible MB1 segment to form seeds:
210  int index = -1;
211  for (SegmentContainer::iterator it = DTlist1.begin(); it != DTlist1.end(); ++it ){
212 
213  index++;
214 
215  if (usedDTlist1[index] == true) continue;
216  if ( int ((*it)->recHits().size()) < minDTHitsPerSegment ) continue;
217  if ((*it)->dimension() != 4) continue;
218 
219  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
220  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
221 
222  // Global position of starting point for protoTrack
223  GlobalPoint gp = (*it)->globalPosition();
224  float eta_temp = gp.eta();
225  float phi_temp = gp.phi();
226  bool showeringBefore = false;
227  NShowerSeg = 0;
228  if ( IdentifyShowering( DTlist1, usedDTlist1, eta_temp, phi_temp, -1, NShowerSeg ) ) showeringBefore = true ;
229  int NShowers = 0;
230  if ( showeringBefore ) {
231  //usedDTlist1[index] = true;
232  NShowers++ ;
233  }
234 
235  SegmentContainer protoTrack;
236  protoTrack.push_back(*it);
237 
238  std::vector<int> layers;
239  layers.push_back(-1);
240 
241  // Try adding segment from other stations
242  if (foundMatchingSegment(3, protoTrack, DTlist2, usedDTlist2, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-2);
243  if ( showeringBefore ) NShowers++ ;
244  if (foundMatchingSegment(3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-3);
245  if ( showeringBefore ) NShowers++ ;
246  if (foundMatchingSegment(3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-4);
247  if ( showeringBefore ) NShowers++ ;
248 
249  // Check if in overlap region
250  if (eta_temp > 0.8) {
251  if (foundMatchingSegment(2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
252  if ( showeringBefore ) NShowers++ ;
253  if (foundMatchingSegment(2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
254  if ( showeringBefore ) NShowers++ ;
255  if (foundMatchingSegment(2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
256  if ( showeringBefore ) NShowers++ ;
257  }
258  else if (eta_temp < -0.8) {
259  if (foundMatchingSegment(2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
260  if ( showeringBefore ) NShowers++ ;
261  if (foundMatchingSegment(2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
262  if ( showeringBefore ) NShowers++ ;
263  if (foundMatchingSegment(2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
264  if ( showeringBefore ) NShowers++ ;
265  }
266 
267 
268  // adding showering information
269  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
270  for (size_t i=0; i< ShoweringSegments.size(); i++) {
271  if ( ShoweringLayers[i] > 0 ) {
272  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
273  protoTrack.push_back( ShoweringSegments[i] );
274  layers.push_back( ShoweringLayers[i] );
275  }
276  if ( ShoweringLayers[i] < 0 && layers[ layers.size()-1] < 0 ) {
277  if ( ShoweringLayers[i] >= layers[ layers.size()-1] ) continue;
278  protoTrack.push_back( ShoweringSegments[i] );
279  layers.push_back( ShoweringLayers[i] );
280  }
281  }
282  }
283  ShoweringSegments.clear() ;
284  ShoweringLayers.clear() ;
285 
286  TrajectorySeed thisSeed;
287 
288  if ( layers.size() < 2 ) {
289  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
290  } else {
291  if ( layers[ layers.size()-1] > 0 ) {
292  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg );
293  } else {
294  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg );
295  }
296  }
297  // Add the seeds to master collection
298  rawSeeds.push_back(thisSeed);
299  etaOfSeed.push_back(eta_temp);
300  phiOfSeed.push_back(phi_temp);
301  nSegOnSeed.push_back( protoTrack.size() );
302 
303  // Marked segment as used
304  usedDTlist1[index] = true;
305  }
306 
307 
308  // Loop over all possible MB2 segment to form seeds:
309  index = -1;
310  for (SegmentContainer::iterator it = DTlist2.begin(); it != DTlist2.end(); ++it ){
311 
312  index++;
313 
314  if (usedDTlist2[index] == true) continue;
315  if ( int ((*it)->recHits().size()) < minDTHitsPerSegment ) continue;
316  if ((*it)->dimension() != 4) continue;
317 
318  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
319  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
320 
321  // Global position of starting point for protoTrack
322  GlobalPoint gp = (*it)->globalPosition();
323  float eta_temp = gp.eta();
324  float phi_temp = gp.phi();
325  bool showeringBefore = false;
326  NShowerSeg = 0;
327  if ( IdentifyShowering( DTlist2, usedDTlist2, eta_temp, phi_temp, -2, NShowerSeg) ) showeringBefore = true ;
328  int NShowers = 0;
329  if ( showeringBefore ) {
330  // usedDTlist2[index] = true;
331  NShowers++ ;
332  }
333 
334  SegmentContainer protoTrack;
335  protoTrack.push_back(*it);
336 
337  std::vector<int> layers;
338  layers.push_back(-2);
339 
340 
341  // Try adding segment from other stations
342  if (foundMatchingSegment(3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-3);
343  if ( showeringBefore ) NShowers++ ;
344  if (foundMatchingSegment(3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-4);
345  if ( showeringBefore ) NShowers++ ;
346 
347  // Check if in overlap region
348  if (eta_temp > 0.8) {
349  if (foundMatchingSegment(2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
350  if ( showeringBefore ) NShowers++ ;
351  if (foundMatchingSegment(2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
352  if ( showeringBefore ) NShowers++ ;
353  if (foundMatchingSegment(2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
354  if ( showeringBefore ) NShowers++ ;
355  }
356  else if (eta_temp < -0.8) {
357  if (foundMatchingSegment(2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
358  if ( showeringBefore ) NShowers++ ;
359  if (foundMatchingSegment(2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
360  if ( showeringBefore ) NShowers++ ;
361  if (foundMatchingSegment(2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
362  if ( showeringBefore ) NShowers++ ;
363  }
364 
365  // adding showering information
366  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
367  for (size_t i=0; i< ShoweringSegments.size(); i++) {
368  if ( ShoweringLayers[i] > 0 ) {
369  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
370  protoTrack.push_back( ShoweringSegments[i] );
371  layers.push_back( ShoweringLayers[i] );
372  }
373  if ( ShoweringLayers[i] < 0 && layers[ layers.size()-1] < 0 ) {
374  if ( ShoweringLayers[i] >= layers[ layers.size()-1] ) continue;
375  protoTrack.push_back( ShoweringSegments[i] );
376  layers.push_back( ShoweringLayers[i] );
377  }
378  }
379  }
380  ShoweringSegments.clear() ;
381  ShoweringLayers.clear() ;
382 
383 
384  TrajectorySeed thisSeed;
385 
386  if ( layers.size() < 2 ) {
387  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
388  } else {
389  if ( layers[ layers.size()-1] > 0 ) {
390  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg );
391  } else {
392  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg );
393  }
394  }
395  // Add the seeds to master collection
396  rawSeeds.push_back(thisSeed);
397  etaOfSeed.push_back(eta_temp);
398  phiOfSeed.push_back(phi_temp);
399  nSegOnSeed.push_back( protoTrack.size() );
400 
401  // Marked segment as used
402  usedDTlist2[index] = true;
403  }
404 
405 
406  // Loop over all possible MB3 segment to form seeds:
407  index = -1;
408  for (SegmentContainer::iterator it = DTlist3.begin(); it != DTlist3.end(); ++it ){
409 
410  index++;
411 
412  if (usedDTlist3[index] == true) continue;
413  if ( int ((*it)->recHits().size()) < minDTHitsPerSegment ) continue;
414  if ((*it)->dimension() != 4) continue;
415 
416  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
417  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
418 
419  // Global position of starting point for protoTrack
420  GlobalPoint gp = (*it)->globalPosition();
421  float eta_temp = gp.eta();
422  float phi_temp = gp.phi();
423  bool showeringBefore = false;
424  NShowerSeg = 0;
425  if ( IdentifyShowering( DTlist3, usedDTlist3, eta_temp, phi_temp, -3, NShowerSeg) ) showeringBefore = true ;
426  int NShowers = 0;
427  if ( showeringBefore ) {
428  // usedDTlist3[index] = true;
429  NShowers++ ;
430  }
431 
432  SegmentContainer protoTrack;
433  protoTrack.push_back(*it);
434 
435  std::vector<int> layers;
436  layers.push_back(-3);
437 
438  // Try adding segment from other stations
439  if (foundMatchingSegment(3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(-4);
440  if ( showeringBefore ) NShowers++ ;
441 
442  // Check if in overlap region
443  if (eta_temp > 0.8) {
444  if (foundMatchingSegment(2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
445  if ( showeringBefore ) NShowers++ ;
446  if (foundMatchingSegment(2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
447  if ( showeringBefore ) NShowers++ ;
448  if (foundMatchingSegment(2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
449  if ( showeringBefore ) NShowers++ ;
450  }
451  else if (eta_temp < -0.8) {
452  if (foundMatchingSegment(2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
453  if ( showeringBefore ) NShowers++ ;
454  if (foundMatchingSegment(2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
455  if ( showeringBefore ) NShowers++ ;
456  if (foundMatchingSegment(2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
457  if ( showeringBefore ) NShowers++ ;
458  }
459 
460  // adding showering information
461  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
462  for (size_t i=0; i< ShoweringSegments.size(); i++) {
463  if ( ShoweringLayers[i] > 0 ) {
464  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
465  protoTrack.push_back( ShoweringSegments[i] );
466  layers.push_back( ShoweringLayers[i] );
467  }
468  if ( ShoweringLayers[i] < 0 && layers[ layers.size()-1] < 0 ) {
469  if ( ShoweringLayers[i] >= layers[ layers.size()-1] ) continue;
470  protoTrack.push_back( ShoweringSegments[i] );
471  layers.push_back( ShoweringLayers[i] );
472  }
473  }
474  }
475  ShoweringSegments.clear() ;
476  ShoweringLayers.clear() ;
477 
478 
479  TrajectorySeed thisSeed;
480  if ( layers.size() < 2 ) {
481  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
482  }else {
483  if ( layers[ layers.size()-1] > 0 ) {
484  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg );
485  } else {
486  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg );
487  }
488  }
489 
490  // Add the seeds to master collection
491  rawSeeds.push_back(thisSeed);
492  etaOfSeed.push_back(eta_temp);
493  phiOfSeed.push_back(phi_temp);
494  nSegOnSeed.push_back( protoTrack.size() );
495 
496  // Marked segment as used
497  usedDTlist3[index] = true;
498  }
499 
500  /* *********************************************************************************************************************
501  * Form seeds from backward endcap
502  *
503  * Proceed from inside -> out
504  * *********************************************************************************************************************/
505 
506  // Loop over all possible ME11 segment to form seeds:
507  index = -1;
508 
509  for (SegmentContainer::iterator it = CSClist0B.begin(); it != CSClist0B.end(); ++it ){
510 
511  index++;
512 
513  if (usedCSClist0B[index] == true) continue;
514  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
515 
516  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
517  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
518 
519  // Global position of starting point for protoTrack
520  GlobalPoint gp = (*it)->globalPosition();
521  float eta_temp = gp.eta();
522  float phi_temp = gp.phi();
523  bool showeringBefore = false;
524  NShowerSeg = 0;
525  if ( IdentifyShowering( CSClist0B, usedCSClist0B, eta_temp, phi_temp, 0, NShowerSeg) ) showeringBefore = true ;
526  int NShowers = 0;
527  if ( showeringBefore ) {
528  // usedCSClist0B[index] = true;
529  NShowers++ ;
530  }
531 
532  SegmentContainer protoTrack;
533  protoTrack.push_back(*it);
534 
535  std::vector<int> layers;
536  layers.push_back(0);
537 
538  // Try adding segment from other station
539  if (foundMatchingSegment(1, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
540  if ( showeringBefore ) NShowers++ ;
541  if (foundMatchingSegment(1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
542  if ( showeringBefore ) NShowers++ ;
543  if (foundMatchingSegment(1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
544  if ( showeringBefore ) NShowers++ ;
545  if (foundMatchingSegment(1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
546  if ( showeringBefore ) NShowers++ ;
547 
548 
549  // adding showering information
550  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
551  for (size_t i=0; i< ShoweringSegments.size(); i++) {
552  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
553  protoTrack.push_back( ShoweringSegments[i] );
554  layers.push_back( ShoweringLayers[i] );
555  }
556  }
557  ShoweringSegments.clear() ;
558  ShoweringLayers.clear() ;
559 
560 
561  TrajectorySeed thisSeed;
562  if ( layers.size() < 2 ) {
563  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
564  }else {
565  if ( fabs( gp.eta() ) > 1.7 ) {
566  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg );
567  } else {
568  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
569  }
570  }
571 
572  // Add the seeds to master collection
573  rawSeeds.push_back(thisSeed);
574  etaOfSeed.push_back(eta_temp);
575  phiOfSeed.push_back(phi_temp);
576  nSegOnSeed.push_back( protoTrack.size() );
577 
578  // mark this segment as used
579  usedCSClist0B[index] = true;
580  }
581 
582 
583  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
584  index = -1;
585  for (SegmentContainer::iterator it = CSClist1B.begin(); it != CSClist1B.end(); ++it ){
586 
587  index++;
588 
589  if (usedCSClist1B[index] == true) continue;
590  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
591 
592  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
593  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
594 
595  // Global position of starting point for protoTrack
596  GlobalPoint gp = (*it)->globalPosition();
597  float eta_temp = gp.eta();
598  float phi_temp = gp.phi();
599  bool showeringBefore = false;
600  NShowerSeg = 0;
601  if ( IdentifyShowering( CSClist1B, usedCSClist1B, eta_temp, phi_temp, 1, NShowerSeg) ) showeringBefore = true ;
602  int NShowers = 0;
603  if ( showeringBefore ) {
604  // usedCSClist1B[index] = true;
605  NShowers++ ;
606  }
607 
608  SegmentContainer protoTrack;
609  protoTrack.push_back(*it);
610 
611  std::vector<int> layers;
612  layers.push_back(1);
613 
614  // Try adding segment from other stations
615  if (foundMatchingSegment(1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
616  if ( showeringBefore ) NShowers++ ;
617  if (foundMatchingSegment(1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
618  if ( showeringBefore ) NShowers++ ;
619  if (foundMatchingSegment(1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
620  if ( showeringBefore ) NShowers++ ;
621 
622  // adding showering information
623  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
624  for (size_t i=0; i< ShoweringSegments.size(); i++) {
625  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
626  protoTrack.push_back( ShoweringSegments[i] );
627  layers.push_back( ShoweringLayers[i] );
628  }
629  }
630  ShoweringSegments.clear() ;
631  ShoweringLayers.clear() ;
632 
633  TrajectorySeed thisSeed;
634  if ( layers.size() < 2 ) {
635  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
636  } else {
637  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
638  }
639  // Add the seeds to master collection
640  rawSeeds.push_back(thisSeed);
641  etaOfSeed.push_back(eta_temp);
642  phiOfSeed.push_back(phi_temp);
643  nSegOnSeed.push_back( protoTrack.size() );
644 
645  // mark this segment as used
646  usedCSClist1B[index] = true;
647 
648  }
649 
650 
651  // Loop over all possible ME2 segment to form seeds:
652  index = -1;
653  for (SegmentContainer::iterator it = CSClist2B.begin(); it != CSClist2B.end(); ++it ){
654 
655  index++;
656 
657  if (usedCSClist2B[index] == true) continue;
658  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
659 
660  double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
661  if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
662 
663  // Global position of starting point for protoTrack
664  GlobalPoint gp = (*it)->globalPosition();
665  float eta_temp = gp.eta();
666  float phi_temp = gp.phi();
667  bool showeringBefore = false;
668  NShowerSeg = 0;
669  if ( IdentifyShowering( CSClist2B, usedCSClist2B, eta_temp, phi_temp, 2, NShowerSeg) ) showeringBefore = true ;
670  int NShowers = 0;
671  if ( showeringBefore ) {
672  // usedCSClist2B[index] = true;
673  NShowers++ ;
674  }
675 
676  SegmentContainer protoTrack;
677  protoTrack.push_back(*it);
678 
679  std::vector<int> layers;
680  layers.push_back(2);
681 
682  // Try adding segment from other stations
683  if (foundMatchingSegment(1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
684  if ( showeringBefore ) NShowers++ ;
685  if (foundMatchingSegment(1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
686  if ( showeringBefore ) NShowers++ ;
687 
688  // adding showering information
689  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
690  for (size_t i=0; i< ShoweringSegments.size(); i++) {
691  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
692  protoTrack.push_back( ShoweringSegments[i] );
693  layers.push_back( ShoweringLayers[i] );
694  }
695  }
696  ShoweringSegments.clear() ;
697  ShoweringLayers.clear() ;
698 
699 
700  TrajectorySeed thisSeed;
701  if ( layers.size() < 2) {
702  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
703  } else {
704  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
705  }
706 
707  // Add the seeds to master collection
708  rawSeeds.push_back(thisSeed);
709  etaOfSeed.push_back(eta_temp);
710  phiOfSeed.push_back(phi_temp);
711  nSegOnSeed.push_back( protoTrack.size() );
712  // mark this segment as used
713  usedCSClist2B[index] = true;
714  }
715 
716  // Loop over all possible ME3 segment to form seeds:
717  index = -1;
718  for (SegmentContainer::iterator it = CSClist3B.begin(); it != CSClist3B.end(); ++it ){
719 
720  index++;
721 
722  if (usedCSClist3B[index] == true) continue;
723  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
724 
725  double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
726  if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
727 
728  // Global position of starting point for protoTrack
729  GlobalPoint gp = (*it)->globalPosition();
730  float eta_temp = gp.eta();
731  float phi_temp = gp.phi();
732  bool showeringBefore = false;
733  NShowerSeg = 0;
734  if ( IdentifyShowering( CSClist3B, usedCSClist3B, eta_temp, phi_temp, 3, NShowerSeg) ) showeringBefore = true ;
735  int NShowers = 0;
736  if ( showeringBefore ) {
737  // usedCSClist3B[index] = true;
738  NShowers++ ;
739  }
740 
741  SegmentContainer protoTrack;
742  protoTrack.push_back(*it);
743 
744  std::vector<int> layers;
745  layers.push_back(2);
746 
747  // Try adding segment from other stations
748  if (foundMatchingSegment(1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
749  if ( showeringBefore ) NShowers++ ;
750 
751 
752  // adding showering information
753  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
754  for (size_t i=0; i< ShoweringSegments.size(); i++) {
755  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
756  protoTrack.push_back( ShoweringSegments[i] );
757  layers.push_back( ShoweringLayers[i] );
758  }
759  }
760  ShoweringSegments.clear() ;
761  ShoweringLayers.clear() ;
762 
763  // mark this segment as used
764  usedCSClist3B[index] = true;
765 
766  if ( layers.size() < 2 ) continue;
767  TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
768 
769  // Add the seeds to master collection
770  rawSeeds.push_back(thisSeed);
771  etaOfSeed.push_back(eta_temp);
772  phiOfSeed.push_back(phi_temp);
773  nSegOnSeed.push_back( protoTrack.size() );
774  }
775 
776 
777  /* *****************************************************************************************************************
778  * Form seeds from forward endcap
779  *
780  * Proceed from inside -> out
781  * *****************************************************************************************************************/
782 
783  // Loop over all possible ME11 segment to form seeds:
784  index = -1;
785  for (SegmentContainer::iterator it = CSClist0F.begin(); it != CSClist0F.end(); ++it ){
786 
787  index++;
788 
789  if (usedCSClist0F[index] == true) continue;
790  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
791 
792  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
793  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
794 
795  // Global position of starting point for protoTrack
796  GlobalPoint gp = (*it)->globalPosition();
797  float eta_temp = gp.eta();
798  float phi_temp = gp.phi();
799  bool showeringBefore = false;
800  NShowerSeg = 0;
801  if ( IdentifyShowering( CSClist0F, usedCSClist0F, eta_temp, phi_temp, 0, NShowerSeg) ) showeringBefore = true ;
802  int NShowers = 0;
803  if ( showeringBefore ) {
804  // usedCSClist0F[index] = true;
805  NShowers++ ;
806  }
807 
808  SegmentContainer protoTrack;
809  protoTrack.push_back(*it);
810 
811  std::vector<int> layers;
812  layers.push_back(0);
813 
814  // Try adding segment from other station
815  if (foundMatchingSegment(1, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(1);
816  if ( showeringBefore ) NShowers++ ;
817  if (foundMatchingSegment(1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
818  if ( showeringBefore ) NShowers++ ;
819  if (foundMatchingSegment(1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
820  if ( showeringBefore ) NShowers++ ;
821  if (foundMatchingSegment(1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
822  if ( showeringBefore ) NShowers++ ;
823 
824 
825  // adding showering information
826  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
827  for (size_t i=0; i< ShoweringSegments.size(); i++) {
828  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
829  protoTrack.push_back( ShoweringSegments[i] );
830  layers.push_back( ShoweringLayers[i] );
831  }
832  }
833  ShoweringSegments.clear() ;
834  ShoweringLayers.clear() ;
835 
836 
837  TrajectorySeed thisSeed;
838  if ( layers.size() < 2 ) {
839  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
840  } else {
841  if ( fabs( gp.eta() ) > 1.7 ) {
842  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg );
843  } else {
844  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
845  }
846  }
847  // Add the seeds to master collection
848  rawSeeds.push_back(thisSeed);
849  etaOfSeed.push_back(eta_temp);
850  phiOfSeed.push_back(phi_temp);
851  nSegOnSeed.push_back( protoTrack.size() );
852 
853  // mark this segment as used
854  usedCSClist0F[index] = true;
855  }
856 
857 
858  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
859  index = -1;
860  for (SegmentContainer::iterator it = CSClist1F.begin(); it != CSClist1F.end(); ++it ){
861 
862  index++;
863 
864  if (usedCSClist1F[index] == true) continue;
865  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
866 
867  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
868  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
869 
870  // Global position of starting point for protoTrack
871  GlobalPoint gp = (*it)->globalPosition();
872  float eta_temp = gp.eta();
873  float phi_temp = gp.phi();
874  bool showeringBefore = false;
875  NShowerSeg = 0;
876  if ( IdentifyShowering( CSClist1F, usedCSClist1F, eta_temp, phi_temp, 1, NShowerSeg) ) showeringBefore = true ;
877  int NShowers = 0;
878  if ( showeringBefore ) {
879  // usedCSClist1F[index] = true;
880  NShowers++ ;
881  }
882 
883  SegmentContainer protoTrack;
884  protoTrack.push_back(*it);
885 
886  std::vector<int> layers;
887  layers.push_back(1);
888 
889  // Try adding segment from other stations
890  if (foundMatchingSegment(1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(2);
891  if ( showeringBefore ) NShowers++ ;
892  if (foundMatchingSegment(1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
893  if ( showeringBefore ) NShowers++ ;
894  if (foundMatchingSegment(1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
895  if ( showeringBefore ) NShowers++ ;
896 
897 
898  // adding showering information
899  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
900  for (size_t i=0; i< ShoweringSegments.size(); i++) {
901  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
902  protoTrack.push_back( ShoweringSegments[i] );
903  layers.push_back( ShoweringLayers[i] );
904  }
905  }
906  ShoweringSegments.clear() ;
907  ShoweringLayers.clear() ;
908 
909 
910  TrajectorySeed thisSeed;
911  if ( layers.size() < 2) {
912  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
913  } else {
914  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
915  }
916 
917  // Add the seeds to master collection
918  rawSeeds.push_back(thisSeed);
919  etaOfSeed.push_back(eta_temp);
920  phiOfSeed.push_back(phi_temp);
921  nSegOnSeed.push_back( protoTrack.size() );
922 
923  // mark this segment as used
924  usedCSClist1F[index] = true;
925  }
926 
927 
928  // Loop over all possible ME2 segment to form seeds:
929  index = -1;
930  for (SegmentContainer::iterator it = CSClist2F.begin(); it != CSClist2F.end(); ++it ){
931 
932  index++;
933 
934  if (usedCSClist2F[index] == true) continue;
935  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
936 
937  double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
938  if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
939 
940  // Global position of starting point for protoTrack
941  GlobalPoint gp = (*it)->globalPosition();
942  float eta_temp = gp.eta();
943  float phi_temp = gp.phi();
944  bool showeringBefore = false;
945  NShowerSeg = 0;
946  if ( IdentifyShowering( CSClist2F, usedCSClist2F, eta_temp, phi_temp, 2, NShowerSeg) ) showeringBefore = true ;
947  int NShowers = 0;
948  if ( showeringBefore ) {
949  // usedCSClist2F[index] = true;
950  NShowers++ ;
951  }
952 
953  SegmentContainer protoTrack;
954  protoTrack.push_back(*it);
955 
956  std::vector<int> layers;
957  layers.push_back(2);
958 
959  // Try adding segment from other stations
960  if (foundMatchingSegment(1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(3);
961  if ( showeringBefore ) NShowers++ ;
962  if (foundMatchingSegment(1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
963  if ( showeringBefore ) NShowers++ ;
964 
965 
966  // adding showering information
967  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
968  for (size_t i=0; i< ShoweringSegments.size(); i++) {
969  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
970  protoTrack.push_back( ShoweringSegments[i] );
971  layers.push_back( ShoweringLayers[i] );
972  }
973  }
974  ShoweringSegments.clear() ;
975  ShoweringLayers.clear() ;
976 
977  TrajectorySeed thisSeed;
978  if ( layers.size() < 2) {
979  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg );
980  } else {
981  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
982  }
983 
984  // Add the seeds to master collection
985  rawSeeds.push_back(thisSeed);
986  etaOfSeed.push_back(eta_temp);
987  phiOfSeed.push_back(phi_temp);
988  nSegOnSeed.push_back( protoTrack.size() );
989 
990  // mark this segment as used
991  usedCSClist2F[index] = true;
992  }
993 
994  // Loop over all possible ME3 segment to form seeds:
995  index = -1;
996  for (SegmentContainer::iterator it = CSClist3F.begin(); it != CSClist3F.end(); ++it ){
997 
998  index++;
999 
1000  if (usedCSClist3F[index] == true) continue;
1001  if ( int ((*it)->recHits().size()) < minCSCHitsPerSegment ) continue;
1002 
1003  double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
1004  if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
1005 
1006  // Global position of starting point for protoTrack
1007  GlobalPoint gp = (*it)->globalPosition();
1008  float eta_temp = gp.eta();
1009  float phi_temp = gp.phi();
1010  bool showeringBefore = false;
1011  NShowerSeg = 0;
1012  if ( IdentifyShowering( CSClist3F, usedCSClist3F, eta_temp, phi_temp, 3, NShowerSeg) ) showeringBefore = true ;
1013  int NShowers = 0;
1014  if ( showeringBefore ) {
1015  // usedCSClist3F[index] = true;
1016  NShowers++ ;
1017  }
1018 
1019  SegmentContainer protoTrack;
1020  protoTrack.push_back(*it);
1021 
1022  std::vector<int> layers;
1023  layers.push_back(2);
1024 
1025  // Try adding segment from other stations
1026  if (foundMatchingSegment(1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size()-1], showeringBefore )) layers.push_back(4);
1027  if ( showeringBefore ) NShowers++ ;
1028 
1029 
1030  // adding showering information
1031  if ( layers.size() < 2 && ShoweringSegments.size() > 0 ) {
1032  for (size_t i=0; i< ShoweringSegments.size(); i++) {
1033  if ( ShoweringLayers[i] <= layers[ layers.size()-1] ) continue;
1034  protoTrack.push_back( ShoweringSegments[i] );
1035  layers.push_back( ShoweringLayers[i] );
1036  }
1037  }
1038  ShoweringSegments.clear() ;
1039  ShoweringLayers.clear() ;
1040 
1041  // mark this segment as used
1042  usedCSClist3F[index] = true;
1043 
1044  if ( layers.size() < 2 ) continue;
1045 
1046  TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg );
1047 
1048  // Add the seeds to master collection
1049  rawSeeds.push_back(thisSeed);
1050  etaOfSeed.push_back(eta_temp);
1051  phiOfSeed.push_back(phi_temp);
1052  nSegOnSeed.push_back( protoTrack.size() );
1053 
1054  }
1055 
1056 
1057  /* *********************************************************************************************************************
1058  * Clean up raw seed collection and pass to master collection
1059  * *********************************************************************************************************************/
1060 
1061  if (debug) std::cout << "*** CLEAN UP " << std::endl;
1062  if (debug) std::cout << "Number of seeds BEFORE " << rawSeeds.size() << std::endl;
1063 
1064 
1065  int goodSeeds = 0;
1066 
1067  theSeeds = muonSeedClean_->seedCleaner(eventSetup,rawSeeds);
1068  goodSeeds = theSeeds.size();
1069 
1070  //std::cout << " === Before Cleaning : " << rawSeeds.size() <<std::endl;
1071  //std::cout << " => AFTER :" << goodSeeds << std::endl;
1072 
1073  if (debug) std::cout << "Number of seeds AFTER " << goodSeeds << std::endl;
1074 
1075 
1076  return goodSeeds;
1077 }
1078 
1079 
1080 
1081 /* *********************************************************************************************************************
1082  * Try to match segment from different station (layer)
1083  *
1084  * Note type = 1 --> endcap
1085  * = 2 --> overlap
1086  * = 3 --> barrel
1087  * *********************************************************************************************************************/
1088 
1090 bool MuonSeedBuilder::foundMatchingSegment( int type, SegmentContainer& protoTrack, SegmentContainer& segs,
1091  BoolContainer& usedSeg, float& eta_last, float& phi_last, int& lastLayer, bool& showeringBefore ) {
1092 
1093  bool ok = false;
1094  int scanlayer = (lastLayer < 0 ) ? (lastLayer-1) : (lastLayer+1) ;
1095 
1096  if ( IdentifyShowering( segs, usedSeg, eta_last, phi_last, scanlayer, NShowerSeg ) ) {
1097  showeringBefore = true;
1098  return ok ;
1099  }
1100 
1101  // Setup the searching cone-size
1102  // searching cone-size should changes with bending power
1103  double maxdEta;
1104  double maxdPhi;
1105  if ( type == 1 ) {
1106  // CSC
1107  maxdEta = maxDeltaEtaCSC;
1108  if ( lastLayer == 0 || lastLayer == 1 ) {
1109  if ( fabs(eta_last) < 2.1 ) { maxdPhi = maxDeltaPhiCSC; }
1110  else { maxdPhi = 0.06; }
1111  }
1112  else if (lastLayer== 2 ) { maxdPhi = 0.5*maxDeltaPhiCSC; }
1113  else { maxdPhi = 0.2*maxDeltaPhiCSC; }
1114 
1115  } else if ( type == 2 ) {
1116  // Overlap
1117  maxdEta = maxDeltaEtaOverlap;
1118  if ( lastLayer == -1 ) { maxdPhi = maxDeltaPhiDT; }
1119  else { maxdPhi = maxDeltaPhiOverlap; }
1120 
1121  } else {
1122  // DT
1123  maxdEta = maxDeltaEtaDT;
1124  if ( lastLayer == -1 ) { maxdPhi = maxDeltaPhiDT; }
1125  else if ( lastLayer == -2 ) { maxdPhi = 0.8*maxDeltaPhiDT; }
1126  else { maxdPhi = 0.4*maxDeltaPhiDT; }
1127 
1128  }
1129 
1130  // if previous layer showers, limite the maxdPhi < 0.06
1131  if ( showeringBefore && maxdPhi > 0.03 ) maxdPhi = 0.03;
1132  // reset the showering flag
1133  showeringBefore = false ;
1134 
1135  // global phi/eta from previous segment
1136  float eta_temp = eta_last;
1137  float phi_temp = phi_last;
1138 
1139  // Index counter to keep track of element used in segs
1140  int index = -1;
1141  int best_match = index;
1142  float best_R = sqrt( (maxdEta*maxdEta) + (maxdPhi*maxdPhi) );
1143  float best_chi2 = 200;
1144  int best_dimension = 2;
1145  int best_nhits = minDTHitsPerSegment;
1146  if( type == 1 ) best_nhits = minCSCHitsPerSegment;
1147  // Loop over segments in other station (layer) and find candidate match
1148  for (SegmentContainer::iterator it=segs.begin(); it!=segs.end(); ++it){
1149 
1150  index++;
1151 
1152  // Not to get confused: eta_last is from the previous layer.
1153  // This is only to find the best set of segments by comparing at the distance layer-by-layer
1154  GlobalPoint gp2 = (*it)->globalPosition();
1155  double dh = fabs( gp2.eta() - eta_temp );
1156  double df = fabs( gp2.phi() - phi_temp );
1157  double dR = sqrt( (dh*dh) + (df*df) );
1158 
1159  // dEta and dPhi should be within certain range
1160  bool case1 = ( dh < maxdEta && df < maxdPhi ) ? true:false ;
1161  // for DT station 4 ; CSCSegment is always 4D
1162  bool case2 = ( ((*it)->dimension()!= 4) && (dh< 0.5) && (df < maxdPhi) ) ? true:false ;
1163  if ( !case1 && !case2 ) continue;
1164 
1165  int NRechits = muonSeedClean_->NRecHitsFromSegment( &*(*it) ) ;
1166 
1167  if ( NRechits < best_nhits ) continue;
1168  best_nhits = NRechits ;
1169 
1170  // reject 2D segments if 4D segments are available
1171  if ( (*it)->dimension() < best_dimension ) continue;
1172  best_dimension = (*it)->dimension();
1173 
1174  // pick the segment with best chi2/dof within a fixed cone size
1175  if ( dR > best_R ) continue;
1176 
1177  // select smaller chi2/dof
1178  double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
1180  if ( (*it)->chi2()/dof < 0.001 && NRechits < 6 && type == 1) continue;
1181  if ( (*it)->chi2()/dof > best_chi2 ) continue;
1182  best_chi2 = (*it)->chi2()/dof ;
1183  best_match = index;
1184  // propagate the eta and phi to next layer
1185  if ((*it)->dimension() != 4 ) {
1186  phi_last = phi_last;
1187  eta_last = eta_last;
1188  } else {
1189  phi_last = gp2.phi();
1190  eta_last = gp2.eta();
1191  }
1192  }
1193 
1194  if (best_match < 0) return ok;
1195 
1196  // Add best matching segment to protoTrack:
1197  index = -1;
1198  for (SegmentContainer::iterator it=segs.begin(); it!=segs.end(); ++it){
1199  index++;
1200  if (index != best_match) continue;
1201  protoTrack.push_back(*it);
1202  usedSeg[best_match] = true;
1203  ok = true;
1204  }
1205  return ok;
1206 }
1207 
1208 
1209 bool MuonSeedBuilder::IdentifyShowering( SegmentContainer& segs, BoolContainer& usedSeg, float& eta_last, float& phi_last, int layer, int& NShoweringSegments ) {
1210 
1211  bool showering = false ;
1212 
1213  int nSeg = 0 ;
1214  int nRhits = 0 ;
1215  double nChi2 = 9999. ;
1216  int theOrigin = -1;
1217  std::vector<int> badtag;
1218  int index = -1;
1219  double aveEta = 0.0;
1220  for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it){
1221 
1222  index++;
1223  GlobalPoint gp = (*it)->globalPosition();
1224  double dh = gp.eta() - eta_last ;
1225  double df = gp.phi() - phi_last ;
1226  double dR = sqrt( (dh*dh) + (df*df) ) ;
1227 
1228  double dof = static_cast<double>( (*it)->degreesOfFreedom() );
1229  double nX2 = (*it)->chi2() / dof ;
1230 
1231  bool isDT = false ;
1232  DetId geoId = (*it)->geographicalId();
1233  if ( geoId.subdetId() == MuonSubdetId::DT ) isDT = true;
1234 
1235  if (dR < 0.3 ) {
1236  nSeg++ ;
1237  badtag.push_back( index ) ;
1238  aveEta += fabs( gp.eta() ) ;
1239  // pick up the best segment from showering chamber
1240  int rh = muonSeedClean_->NRecHitsFromSegment( &*(*it) );
1241  if (rh < 6 && !isDT) continue;
1242  if (rh < 12 && isDT) continue;
1243  if ( rh > nRhits ) {
1244  nRhits = rh ;
1245  if ( nX2 > nChi2 ) continue ;
1246  if (layer != 0 && layer != 1 && layer != -1 ) {
1247  theOrigin = index ;
1248  }
1249  }
1250  }
1251 
1252  }
1253  aveEta = aveEta/static_cast<double>(nSeg) ;
1254  bool isME11A = (aveEta >= 2.1 && layer == 0) ? true : false ;
1255  bool isME12 = (aveEta > 1.2 && aveEta <= 1.65 && layer == 1) ? true : false ;
1256  bool isME11 = (aveEta > 1.65 && aveEta <= 2.1 && layer == 0) ? true : false ;
1257  bool is1stLayer = (layer == -1 || layer == 0 || isME12 || isME11 || isME11A) ? true : false ;
1258 
1259  NShoweringSegments += nSeg;
1260 
1261  if ( nSeg > 3 && !isME11A ) showering = true ;
1262  if ( nSeg > 6 && isME11A ) showering = true ;
1263 
1264  // if showering, flag all segments in order to skip this layer for pt estimation except 1st layer
1265  //std::cout<<" from Showering "<<std::endl;
1266  if (showering && !is1stLayer ) {
1267  for (std::vector<int>::iterator it = badtag.begin(); it != badtag.end(); ++it ) {
1268  usedSeg[*it] = true;
1269  if ( (*it) != theOrigin ) continue;
1270  ShoweringSegments.push_back( segs[*it] );
1271  ShoweringLayers.push_back( layer );
1272  }
1273  }
1274  return showering ;
1275 
1276 }
1277 
1278 double MuonSeedBuilder::etaError(const GlobalPoint gp, double rErr) {
1279 
1280  double dHdTheta = 0.0;
1281  double dThetadR = 0.0;
1282  double etaErr = 1.0;
1283 
1284  if (gp.perp() != 0) {
1285 
1286  dHdTheta = ( gp.mag()+gp.z() )/gp.perp();
1287  dThetadR = gp.z() / gp.perp2() ;
1288  etaErr = 0.25 * (dHdTheta * dThetadR) * (dHdTheta * dThetadR) * rErr ;
1289  }
1290 
1291  return etaErr;
1292 }
1293 
MuonServiceProxy_cff.MuonServiceProxy
tuple MuonServiceProxy
Definition: MuonServiceProxy_cff.py:15
MuonSeedBuilder::debug
bool debug
group the seeds
Definition: MuonSeedBuilder.h:105
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
MuonSeedBuilder::etaError
double etaError(const GlobalPoint gp, double rErr)
calculate the eta error from global R error
Definition: MuonSeedBuilder.cc:1278
case2
#define case2(par1, par2)
Definition: vmac.h:92
LayerTriplets::layers
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
Definition: LayerTriplets.cc:4
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
MuonSeedBuilder::SegmentContainer
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
Definition: MuonSeedBuilder.h:38
MuonSeedBuilder::ShoweringSegments
SegmentContainer ShoweringSegments
Definition: MuonSeedBuilder.h:129
MuonSeedBuilder::IdentifyShowering
bool IdentifyShowering(SegmentContainer &segs, BoolContainer &usedSeg, float &eta_last, float &phi_last, int layer, int &NShoweringSegments)
identify the showering layer
Definition: MuonSeedBuilder.cc:1209
MuonSeedBuilder::build
int build(edm::Event &event, const edm::EventSetup &eventSetup, TrajectorySeedCollection &seeds)
Build seed collection.
Definition: MuonSeedBuilder.cc:121
MuonSeedBuilder::foundMatchingSegment
bool foundMatchingSegment(int type, SegmentContainer &protoTrack, SegmentContainer &segments, BoolContainer &usedSeg, float &eta_temp, float &phi_temp, int &lastLayer, bool &showeringBefore)
Find segment which matches protoTrack for endcap only.
Definition: MuonSeedBuilder.cc:1090
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
MuonSeedCreator::createSeed
TrajectorySeed createSeed(int type, const SegmentContainer &seg, const std::vector< int > &layers, int NShower, int NShowerSeg)
Create a seed from set of segments.
Definition: MuonSeedCreator.cc:153
PV3DBase::perp2
T perp2() const
Definition: PV3DBase.h:71
MuonSeedBuilder::muonMeasurements
MuonDetLayerMeasurements * muonMeasurements
Definition: MuonSeedBuilder.h:158
MuonSeedBuilder::muonSeedClean_
MuonSeedCleaner * muonSeedClean_
Definition: MuonSeedBuilder.h:139
MuonSeedBuilder::BField
const MagneticField * BField
Definition: MuonSeedBuilder.h:145
MuonSeedBuilder::theCSCSegmentLabel
edm::InputTag theCSCSegmentLabel
Name of the CSC segment collection.
Definition: MuonSeedBuilder.h:135
MuonSeedBuilder::BoolContainer
std::deque< bool > BoolContainer
Definition: MuonSeedBuilder.h:39
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
MuonSeedBuilder::theDTSegmentLabel
edm::InputTag theDTSegmentLabel
Name of the DT segment collection.
Definition: MuonSeedBuilder.h:132
MuonSeedBuilder::ShoweringLayers
std::vector< int > ShoweringLayers
Definition: MuonSeedBuilder.h:130
TrajectorySeedCollection
std::vector< TrajectorySeed > TrajectorySeedCollection
Definition: TrajectorySeedCollection.h:6
MuonDetLayerMeasurements::recHits
MuonRecHitContainer recHits(const DetLayer *layer, const edm::Event &iEvent)
returns the rechits which are on the layer
Definition: MuonDetLayerMeasurements.cc:345
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
MuonSeedCreator::setBField
void setBField(const MagneticField *theField)
Cache Magnetic Field for current event.
Definition: MuonSeedCreator.h:45
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
DetId
Definition: DetId.h:18
MuonSeedCleaner::NRecHitsFromSegment
int NRecHitsFromSegment(const TrackingRecHit &rhit)
Definition: MuonSeedCleaner.cc:481
MuonDetLayerGeometry::allDTLayers
const std::vector< const DetLayer * > & allDTLayers() const
return the DT DetLayers (barrel), inside-out
Definition: MuonDetLayerGeometry.cc:131
case1
#define case1(par)
Definition: vmac.h:81
MuonDetLayerGeometry::forwardCSCLayers
const std::vector< const DetLayer * > & forwardCSCLayers() const
return the forward (+Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:142
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
MuonSeedBuilder::MuonSeedBuilder
MuonSeedBuilder(const edm::ParameterSet &, edm::ConsumesCollector &)
Constructor.
Definition: MuonSeedBuilder.cc:52
MuonSeedBuilder::theService
MuonServiceProxy * theService
Definition: MuonSeedBuilder.h:148
MuonSeedBuilder::~MuonSeedBuilder
~MuonSeedBuilder()
Destructor.
Definition: MuonSeedBuilder.cc:104
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
MuonSubdetId::DT
static const int DT
Definition: MuonSubdetId.h:12
MuonSeedCleaner::seedCleaner
std::vector< TrajectorySeed > seedCleaner(const edm::EventSetup &eventSetup, std::vector< TrajectorySeed > &seeds)
Cache pointer to geometry.
Definition: MuonSeedCleaner.cc:74
MuonDetLayerGeometry::backwardCSCLayers
const std::vector< const DetLayer * > & backwardCSCLayers() const
return the backward (-Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:148
MuonSeedBuilder::muonSeedCreate_
MuonSeedCreator * muonSeedCreate_
Create seed according to region (CSC, DT, Overlap)
Definition: MuonSeedBuilder.h:138
MuonSeedBuilder::muonLayers
const MuonDetLayerGeometry * muonLayers
Definition: MuonSeedBuilder.h:142
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
cuy.dh
tuple dh
Definition: cuy.py:353
Generated for CMSSW Reference Manual by doxygen 1.8.5