CMS 3D CMS Logo
MuonSeedBuilder.cc
Go to the documentation of this file.
1 
7 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedBuilder.h>
8 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCreator.h>
9 #include <RecoMuon/MuonSeedGenerator/src/MuonSeedCleaner.h>
10 
11 // Data Formats
12 #include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
13 #include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
14 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
15 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
16 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
17 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
18 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
19 #include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>
20 
21 // Geometry
22 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
23 #include <TrackingTools/DetLayers/interface/DetLayer.h>
24 #include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
25 #include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
26 #include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>
27 
28 // muon service
29 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
30 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
31 
32 // Framework
33 #include <FWCore/ParameterSet/interface/ParameterSet.h>
34 #include "FWCore/Framework/interface/Event.h"
35 #include <FWCore/Framework/interface/ESHandle.h>
36 #include <FWCore/MessageLogger/interface/MessageLogger.h>
37 #include <DataFormats/Common/interface/Handle.h>
38 
39 // C++
40 #include <vector>
41 #include <deque>
42 #include <utility>
43 
44 //typedef std::pair<double, TrajectorySeed> seedpr ;
45 //static bool ptDecreasing(const seedpr s1, const seedpr s2) { return ( s1.first > s2.first ); }
46 //static bool lengthSorting(const TrajectorySeed s1, const TrajectorySeed s2) { return ( s1.nHits() > s2.nHits() ); }
47 
48 /*
49  * Constructor
50  */
51 MuonSeedBuilder::MuonSeedBuilder(const edm::ParameterSet& pset, edm::ConsumesCollector& iC) {
52  // Local Debug flag
53  debug = pset.getParameter<bool>("DebugMuonSeed");
54 
55  // enable the DT chamber
56  enableDTMeasurement = pset.getParameter<bool>("EnableDTMeasurement");
57  theDTSegmentLabel = pset.getParameter<edm::InputTag>("DTSegmentLabel");
58 
59  // enable the CSC chamber
60  enableCSCMeasurement = pset.getParameter<bool>("EnableCSCMeasurement");
61  theCSCSegmentLabel = pset.getParameter<edm::InputTag>("CSCSegmentLabel");
62 
63  // Parameters for seed creation in endcap region
64  minCSCHitsPerSegment = pset.getParameter<int>("minCSCHitsPerSegment");
65  maxDeltaEtaCSC = pset.getParameter<double>("maxDeltaEtaCSC");
66  maxDeltaPhiCSC = pset.getParameter<double>("maxDeltaPhiCSC");
67 
68  // Parameters for seed creation in overlap region
69  maxDeltaEtaOverlap = pset.getParameter<double>("maxDeltaEtaOverlap");
70  maxDeltaPhiOverlap = pset.getParameter<double>("maxDeltaPhiOverlap");
71 
72  // Parameters for seed creation in barrel region
73  minDTHitsPerSegment = pset.getParameter<int>("minDTHitsPerSegment");
74  maxDeltaEtaDT = pset.getParameter<double>("maxDeltaEtaDT");
75  maxDeltaPhiDT = pset.getParameter<double>("maxDeltaPhiDT");
76 
77  // Parameters to suppress combinatorics (ghosts)
78  maxEtaResolutionDT = pset.getParameter<double>("maxEtaResolutionDT");
79  maxPhiResolutionDT = pset.getParameter<double>("maxPhiResolutionDT");
80  maxEtaResolutionCSC = pset.getParameter<double>("maxEtaResolutionCSC");
81  maxPhiResolutionCSC = pset.getParameter<double>("maxPhiResolutionCSC");
82 
83  // Class for forming muon seeds
84  muonSeedCreate_ = new MuonSeedCreator(pset);
85  muonSeedClean_ = new MuonSeedCleaner(pset, edm::ConsumesCollector(iC));
86 
87  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
88  muonMeasurements = new MuonDetLayerMeasurements(theDTSegmentLabel,
89  theCSCSegmentLabel,
90  edm::InputTag(),
91  edm::InputTag(),
92  edm::InputTag(),
93  iC,
94  enableDTMeasurement,
95  enableCSCMeasurement,
96  false,
97  false,
98  false);
99 }
100 
101 /*
102  * Destructor
103  */
104 MuonSeedBuilder::~MuonSeedBuilder() {
105  delete muonSeedCreate_;
106  delete muonSeedClean_;
107  if (muonMeasurements)
108  delete muonMeasurements;
109 }
110 
111 /*
112  * build
113  *
114  * Where the segments are sorted out to make a protoTrack (vector of matching segments in different
115  * stations/layers). The protoTrack is then passed to the seed creator to create CSC, DT or overlap seeds
116  *
117  */
118 //void MuonSeedBuilder::build( MuonDetLayerMeasurements muonMeasurements, TrajectorySeedCollection& theSeeds ) {
119 int MuonSeedBuilder::build(edm::Event& event, const edm::EventSetup& eventSetup, TrajectorySeedCollection& theSeeds) {
120  // Pass the Magnetic Field to where the seed is create
121  muonSeedCreate_->setBField(BField);
122 
123  // Create temporary collection of seeds which will be cleaned up to remove combinatorics
124  std::vector<TrajectorySeed> rawSeeds;
125  std::vector<float> etaOfSeed;
126  std::vector<float> phiOfSeed;
127  std::vector<int> nSegOnSeed;
128 
129  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
130  // MuonDetLayerMeasurements muonMeasurements(enableDTMeasurement,enableCSCMeasurement,false,
131  // theDTSegmentLabel.label(),theCSCSegmentLabel.label());
132 
133  // 1) Get the various stations and store segments in containers for each station (layers)
134 
135  // 1a. get the DT segments by stations (layers):
136  std::vector<const DetLayer*> dtLayers = muonLayers->allDTLayers();
137 
138  SegmentContainer DTlist4 = muonMeasurements->recHits(dtLayers[3], event);
139  SegmentContainer DTlist3 = muonMeasurements->recHits(dtLayers[2], event);
140  SegmentContainer DTlist2 = muonMeasurements->recHits(dtLayers[1], event);
141  SegmentContainer DTlist1 = muonMeasurements->recHits(dtLayers[0], event);
142 
143  // Initialize flags that a given segment has been allocated to a seed
144  BoolContainer usedDTlist4(DTlist4.size(), false);
145  BoolContainer usedDTlist3(DTlist3.size(), false);
146  BoolContainer usedDTlist2(DTlist2.size(), false);
147  BoolContainer usedDTlist1(DTlist1.size(), false);
148 
149  if (debug) {
150  std::cout << "*** Number of DT segments is: " << DTlist4.size() + DTlist3.size() + DTlist2.size() + DTlist1.size()
151  << std::endl;
152  std::cout << "In MB1: " << DTlist1.size() << std::endl;
153  std::cout << "In MB2: " << DTlist2.size() << std::endl;
154  std::cout << "In MB3: " << DTlist3.size() << std::endl;
155  std::cout << "In MB4: " << DTlist4.size() << std::endl;
156  }
157 
158  // 1b. get the CSC segments by stations (layers):
159  // 1b.1 Global z < 0
160  std::vector<const DetLayer*> cscBackwardLayers = muonLayers->backwardCSCLayers();
161  SegmentContainer CSClist4B = muonMeasurements->recHits(cscBackwardLayers[4], event);
162  SegmentContainer CSClist3B = muonMeasurements->recHits(cscBackwardLayers[3], event);
163  SegmentContainer CSClist2B = muonMeasurements->recHits(cscBackwardLayers[2], event);
164  SegmentContainer CSClist1B = muonMeasurements->recHits(cscBackwardLayers[1], event); // ME1/2 and 1/3
165  SegmentContainer CSClist0B = muonMeasurements->recHits(cscBackwardLayers[0], event); // ME11
166 
167  BoolContainer usedCSClist4B(CSClist4B.size(), false);
168  BoolContainer usedCSClist3B(CSClist3B.size(), false);
169  BoolContainer usedCSClist2B(CSClist2B.size(), false);
170  BoolContainer usedCSClist1B(CSClist1B.size(), false);
171  BoolContainer usedCSClist0B(CSClist0B.size(), false);
172 
173  // 1b.2 Global z > 0
174  std::vector<const DetLayer*> cscForwardLayers = muonLayers->forwardCSCLayers();
175  SegmentContainer CSClist4F = muonMeasurements->recHits(cscForwardLayers[4], event);
176  SegmentContainer CSClist3F = muonMeasurements->recHits(cscForwardLayers[3], event);
177  SegmentContainer CSClist2F = muonMeasurements->recHits(cscForwardLayers[2], event);
178  SegmentContainer CSClist1F = muonMeasurements->recHits(cscForwardLayers[1], event);
179  SegmentContainer CSClist0F = muonMeasurements->recHits(cscForwardLayers[0], event);
180 
181  BoolContainer usedCSClist4F(CSClist4F.size(), false);
182  BoolContainer usedCSClist3F(CSClist3F.size(), false);
183  BoolContainer usedCSClist2F(CSClist2F.size(), false);
184  BoolContainer usedCSClist1F(CSClist1F.size(), false);
185  BoolContainer usedCSClist0F(CSClist0F.size(), false);
186 
187  if (debug) {
188  std::cout << "*** Number of CSC segments is "
189  << CSClist4F.size() + CSClist3F.size() + CSClist2F.size() + CSClist1F.size() + CSClist0F.size() +
190  CSClist4B.size() + CSClist3B.size() + CSClist2B.size() + CSClist1B.size() + CSClist0B.size()
191  << std::endl;
192  std::cout << "In ME11: " << CSClist0F.size() + CSClist0B.size() << std::endl;
193  std::cout << "In ME12: " << CSClist1F.size() + CSClist1B.size() << std::endl;
194  std::cout << "In ME2 : " << CSClist2F.size() + CSClist2B.size() << std::endl;
195  std::cout << "In ME3 : " << CSClist3F.size() + CSClist3B.size() << std::endl;
196  std::cout << "In ME4 : " << CSClist4F.size() + CSClist4B.size() << std::endl;
197  }
198 
199  /* ******************************************************************************************************************
200  * Form seeds in barrel region
201  *
202  * Proceed from inside -> out
203  * ******************************************************************************************************************/
204 
205  // Loop over all possible MB1 segment to form seeds:
206  int index = -1;
207  for (SegmentContainer::iterator it = DTlist1.begin(); it != DTlist1.end(); ++it) {
208  index++;
209 
210  if (usedDTlist1[index] == true)
211  continue;
212  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
213  continue;
214  if ((*it)->dimension() != 4)
215  continue;
216 
217  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
218  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
219 
220  // Global position of starting point for protoTrack
221  GlobalPoint gp = (*it)->globalPosition();
222  float eta_temp = gp.eta();
223  float phi_temp = gp.phi();
224  bool showeringBefore = false;
225  NShowerSeg = 0;
226  if (IdentifyShowering(DTlist1, usedDTlist1, eta_temp, phi_temp, -1, NShowerSeg))
227  showeringBefore = true;
228  int NShowers = 0;
229  if (showeringBefore) {
230  //usedDTlist1[index] = true;
231  NShowers++;
232  }
233 
234  SegmentContainer protoTrack;
235  protoTrack.push_back(*it);
236 
237  std::vector<int> layers;
238  layers.push_back(-1);
239 
240  // Try adding segment from other stations
241  if (foundMatchingSegment(
242  3, protoTrack, DTlist2, usedDTlist2, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
243  layers.push_back(-2);
244  if (showeringBefore)
245  NShowers++;
246  if (foundMatchingSegment(
247  3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
248  layers.push_back(-3);
249  if (showeringBefore)
250  NShowers++;
251  if (foundMatchingSegment(
252  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
253  layers.push_back(-4);
254  if (showeringBefore)
255  NShowers++;
256 
257  // Check if in overlap region
258  if (eta_temp > 0.8) {
259  if (foundMatchingSegment(
260  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
261  layers.push_back(1);
262  if (showeringBefore)
263  NShowers++;
264  if (foundMatchingSegment(
265  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
266  layers.push_back(2);
267  if (showeringBefore)
268  NShowers++;
269  if (foundMatchingSegment(
270  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
271  layers.push_back(3);
272  if (showeringBefore)
273  NShowers++;
274  } else if (eta_temp < -0.8) {
275  if (foundMatchingSegment(
276  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
277  layers.push_back(1);
278  if (showeringBefore)
279  NShowers++;
280  if (foundMatchingSegment(
281  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
282  layers.push_back(2);
283  if (showeringBefore)
284  NShowers++;
285  if (foundMatchingSegment(
286  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
287  layers.push_back(3);
288  if (showeringBefore)
289  NShowers++;
290  }
291 
292  // adding showering information
293  if (layers.size() < 2 && !ShoweringSegments.empty()) {
294  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
295  if (ShoweringLayers[i] > 0) {
296  if (ShoweringLayers[i] <= layers[layers.size() - 1])
297  continue;
298  protoTrack.push_back(ShoweringSegments[i]);
299  layers.push_back(ShoweringLayers[i]);
300  }
301  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
302  if (ShoweringLayers[i] >= layers[layers.size() - 1])
303  continue;
304  protoTrack.push_back(ShoweringSegments[i]);
305  layers.push_back(ShoweringLayers[i]);
306  }
307  }
308  }
309  ShoweringSegments.clear();
310  ShoweringLayers.clear();
311 
312  TrajectorySeed thisSeed;
313 
314  if (layers.size() < 2) {
315  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
316  } else {
317  if (layers[layers.size() - 1] > 0) {
318  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
319  } else {
320  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
321  }
322  }
323  // Add the seeds to master collection
324  rawSeeds.push_back(thisSeed);
325  etaOfSeed.push_back(eta_temp);
326  phiOfSeed.push_back(phi_temp);
327  nSegOnSeed.push_back(protoTrack.size());
328 
329  // Marked segment as used
330  usedDTlist1[index] = true;
331  }
332 
333  // Loop over all possible MB2 segment to form seeds:
334  index = -1;
335  for (SegmentContainer::iterator it = DTlist2.begin(); it != DTlist2.end(); ++it) {
336  index++;
337 
338  if (usedDTlist2[index] == true)
339  continue;
340  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
341  continue;
342  if ((*it)->dimension() != 4)
343  continue;
344 
345  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
346  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
347 
348  // Global position of starting point for protoTrack
349  GlobalPoint gp = (*it)->globalPosition();
350  float eta_temp = gp.eta();
351  float phi_temp = gp.phi();
352  bool showeringBefore = false;
353  NShowerSeg = 0;
354  if (IdentifyShowering(DTlist2, usedDTlist2, eta_temp, phi_temp, -2, NShowerSeg))
355  showeringBefore = true;
356  int NShowers = 0;
357  if (showeringBefore) {
358  // usedDTlist2[index] = true;
359  NShowers++;
360  }
361 
362  SegmentContainer protoTrack;
363  protoTrack.push_back(*it);
364 
365  std::vector<int> layers;
366  layers.push_back(-2);
367 
368  // Try adding segment from other stations
369  if (foundMatchingSegment(
370  3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
371  layers.push_back(-3);
372  if (showeringBefore)
373  NShowers++;
374  if (foundMatchingSegment(
375  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
376  layers.push_back(-4);
377  if (showeringBefore)
378  NShowers++;
379 
380  // Check if in overlap region
381  if (eta_temp > 0.8) {
382  if (foundMatchingSegment(
383  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
384  layers.push_back(1);
385  if (showeringBefore)
386  NShowers++;
387  if (foundMatchingSegment(
388  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
389  layers.push_back(2);
390  if (showeringBefore)
391  NShowers++;
392  if (foundMatchingSegment(
393  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
394  layers.push_back(3);
395  if (showeringBefore)
396  NShowers++;
397  } else if (eta_temp < -0.8) {
398  if (foundMatchingSegment(
399  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
400  layers.push_back(1);
401  if (showeringBefore)
402  NShowers++;
403  if (foundMatchingSegment(
404  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
405  layers.push_back(2);
406  if (showeringBefore)
407  NShowers++;
408  if (foundMatchingSegment(
409  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
410  layers.push_back(3);
411  if (showeringBefore)
412  NShowers++;
413  }
414 
415  // adding showering information
416  if (layers.size() < 2 && !ShoweringSegments.empty()) {
417  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
418  if (ShoweringLayers[i] > 0) {
419  if (ShoweringLayers[i] <= layers[layers.size() - 1])
420  continue;
421  protoTrack.push_back(ShoweringSegments[i]);
422  layers.push_back(ShoweringLayers[i]);
423  }
424  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
425  if (ShoweringLayers[i] >= layers[layers.size() - 1])
426  continue;
427  protoTrack.push_back(ShoweringSegments[i]);
428  layers.push_back(ShoweringLayers[i]);
429  }
430  }
431  }
432  ShoweringSegments.clear();
433  ShoweringLayers.clear();
434 
435  TrajectorySeed thisSeed;
436 
437  if (layers.size() < 2) {
438  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
439  } else {
440  if (layers[layers.size() - 1] > 0) {
441  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
442  } else {
443  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
444  }
445  }
446  // Add the seeds to master collection
447  rawSeeds.push_back(thisSeed);
448  etaOfSeed.push_back(eta_temp);
449  phiOfSeed.push_back(phi_temp);
450  nSegOnSeed.push_back(protoTrack.size());
451 
452  // Marked segment as used
453  usedDTlist2[index] = true;
454  }
455 
456  // Loop over all possible MB3 segment to form seeds:
457  index = -1;
458  for (SegmentContainer::iterator it = DTlist3.begin(); it != DTlist3.end(); ++it) {
459  index++;
460 
461  if (usedDTlist3[index] == true)
462  continue;
463  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
464  continue;
465  if ((*it)->dimension() != 4)
466  continue;
467 
468  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
469  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
470 
471  // Global position of starting point for protoTrack
472  GlobalPoint gp = (*it)->globalPosition();
473  float eta_temp = gp.eta();
474  float phi_temp = gp.phi();
475  bool showeringBefore = false;
476  NShowerSeg = 0;
477  if (IdentifyShowering(DTlist3, usedDTlist3, eta_temp, phi_temp, -3, NShowerSeg))
478  showeringBefore = true;
479  int NShowers = 0;
480  if (showeringBefore) {
481  // usedDTlist3[index] = true;
482  NShowers++;
483  }
484 
485  SegmentContainer protoTrack;
486  protoTrack.push_back(*it);
487 
488  std::vector<int> layers;
489  layers.push_back(-3);
490 
491  // Try adding segment from other stations
492  if (foundMatchingSegment(
493  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
494  layers.push_back(-4);
495  if (showeringBefore)
496  NShowers++;
497 
498  // Check if in overlap region
499  if (eta_temp > 0.8) {
500  if (foundMatchingSegment(
501  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
502  layers.push_back(1);
503  if (showeringBefore)
504  NShowers++;
505  if (foundMatchingSegment(
506  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
507  layers.push_back(2);
508  if (showeringBefore)
509  NShowers++;
510  if (foundMatchingSegment(
511  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
512  layers.push_back(3);
513  if (showeringBefore)
514  NShowers++;
515  } else if (eta_temp < -0.8) {
516  if (foundMatchingSegment(
517  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
518  layers.push_back(1);
519  if (showeringBefore)
520  NShowers++;
521  if (foundMatchingSegment(
522  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
523  layers.push_back(2);
524  if (showeringBefore)
525  NShowers++;
526  if (foundMatchingSegment(
527  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
528  layers.push_back(3);
529  if (showeringBefore)
530  NShowers++;
531  }
532 
533  // adding showering information
534  if (layers.size() < 2 && !ShoweringSegments.empty()) {
535  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
536  if (ShoweringLayers[i] > 0) {
537  if (ShoweringLayers[i] <= layers[layers.size() - 1])
538  continue;
539  protoTrack.push_back(ShoweringSegments[i]);
540  layers.push_back(ShoweringLayers[i]);
541  }
542  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
543  if (ShoweringLayers[i] >= layers[layers.size() - 1])
544  continue;
545  protoTrack.push_back(ShoweringSegments[i]);
546  layers.push_back(ShoweringLayers[i]);
547  }
548  }
549  }
550  ShoweringSegments.clear();
551  ShoweringLayers.clear();
552 
553  TrajectorySeed thisSeed;
554  if (layers.size() < 2) {
555  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
556  } else {
557  if (layers[layers.size() - 1] > 0) {
558  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
559  } else {
560  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
561  }
562  }
563 
564  // Add the seeds to master collection
565  rawSeeds.push_back(thisSeed);
566  etaOfSeed.push_back(eta_temp);
567  phiOfSeed.push_back(phi_temp);
568  nSegOnSeed.push_back(protoTrack.size());
569 
570  // Marked segment as used
571  usedDTlist3[index] = true;
572  }
573 
574  /* *********************************************************************************************************************
575  * Form seeds from backward endcap
576  *
577  * Proceed from inside -> out
578  * *********************************************************************************************************************/
579 
580  // Loop over all possible ME11 segment to form seeds:
581  index = -1;
582 
583  for (SegmentContainer::iterator it = CSClist0B.begin(); it != CSClist0B.end(); ++it) {
584  index++;
585 
586  if (usedCSClist0B[index] == true)
587  continue;
588  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
589  continue;
590 
591  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
592  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
593 
594  // Global position of starting point for protoTrack
595  GlobalPoint gp = (*it)->globalPosition();
596  float eta_temp = gp.eta();
597  float phi_temp = gp.phi();
598  bool showeringBefore = false;
599  NShowerSeg = 0;
600  if (IdentifyShowering(CSClist0B, usedCSClist0B, eta_temp, phi_temp, 0, NShowerSeg))
601  showeringBefore = true;
602  int NShowers = 0;
603  if (showeringBefore) {
604  // usedCSClist0B[index] = true;
605  NShowers++;
606  }
607 
608  SegmentContainer protoTrack;
609  protoTrack.push_back(*it);
610 
611  std::vector<int> layers;
612  layers.push_back(0);
613 
614  // Try adding segment from other station
615  if (foundMatchingSegment(
616  1, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
617  layers.push_back(1);
618  if (showeringBefore)
619  NShowers++;
620  if (foundMatchingSegment(
621  1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
622  layers.push_back(2);
623  if (showeringBefore)
624  NShowers++;
625  if (foundMatchingSegment(
626  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
627  layers.push_back(3);
628  if (showeringBefore)
629  NShowers++;
630  if (foundMatchingSegment(
631  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
632  layers.push_back(4);
633  if (showeringBefore)
634  NShowers++;
635 
636  // adding showering information
637  if (layers.size() < 2 && !ShoweringSegments.empty()) {
638  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
639  if (ShoweringLayers[i] <= layers[layers.size() - 1])
640  continue;
641  protoTrack.push_back(ShoweringSegments[i]);
642  layers.push_back(ShoweringLayers[i]);
643  }
644  }
645  ShoweringSegments.clear();
646  ShoweringLayers.clear();
647 
648  TrajectorySeed thisSeed;
649  if (layers.size() < 2) {
650  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
651  } else {
652  if (fabs(gp.eta()) > 1.7) {
653  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
654  } else {
655  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
656  }
657  }
658 
659  // Add the seeds to master collection
660  rawSeeds.push_back(thisSeed);
661  etaOfSeed.push_back(eta_temp);
662  phiOfSeed.push_back(phi_temp);
663  nSegOnSeed.push_back(protoTrack.size());
664 
665  // mark this segment as used
666  usedCSClist0B[index] = true;
667  }
668 
669  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
670  index = -1;
671  for (SegmentContainer::iterator it = CSClist1B.begin(); it != CSClist1B.end(); ++it) {
672  index++;
673 
674  if (usedCSClist1B[index] == true)
675  continue;
676  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
677  continue;
678 
679  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
680  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
681 
682  // Global position of starting point for protoTrack
683  GlobalPoint gp = (*it)->globalPosition();
684  float eta_temp = gp.eta();
685  float phi_temp = gp.phi();
686  bool showeringBefore = false;
687  NShowerSeg = 0;
688  if (IdentifyShowering(CSClist1B, usedCSClist1B, eta_temp, phi_temp, 1, NShowerSeg))
689  showeringBefore = true;
690  int NShowers = 0;
691  if (showeringBefore) {
692  // usedCSClist1B[index] = true;
693  NShowers++;
694  }
695 
696  SegmentContainer protoTrack;
697  protoTrack.push_back(*it);
698 
699  std::vector<int> layers;
700  layers.push_back(1);
701 
702  // Try adding segment from other stations
703  if (foundMatchingSegment(
704  1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
705  layers.push_back(2);
706  if (showeringBefore)
707  NShowers++;
708  if (foundMatchingSegment(
709  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
710  layers.push_back(3);
711  if (showeringBefore)
712  NShowers++;
713  if (foundMatchingSegment(
714  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
715  layers.push_back(4);
716  if (showeringBefore)
717  NShowers++;
718 
719  // adding showering information
720  if (layers.size() < 2 && !ShoweringSegments.empty()) {
721  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
722  if (ShoweringLayers[i] <= layers[layers.size() - 1])
723  continue;
724  protoTrack.push_back(ShoweringSegments[i]);
725  layers.push_back(ShoweringLayers[i]);
726  }
727  }
728  ShoweringSegments.clear();
729  ShoweringLayers.clear();
730 
731  TrajectorySeed thisSeed;
732  if (layers.size() < 2) {
733  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
734  } else {
735  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
736  }
737  // Add the seeds to master collection
738  rawSeeds.push_back(thisSeed);
739  etaOfSeed.push_back(eta_temp);
740  phiOfSeed.push_back(phi_temp);
741  nSegOnSeed.push_back(protoTrack.size());
742 
743  // mark this segment as used
744  usedCSClist1B[index] = true;
745  }
746 
747  // Loop over all possible ME2 segment to form seeds:
748  index = -1;
749  for (SegmentContainer::iterator it = CSClist2B.begin(); it != CSClist2B.end(); ++it) {
750  index++;
751 
752  if (usedCSClist2B[index] == true)
753  continue;
754  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
755  continue;
756 
757  double dof = static_cast<double>((*it)->degreesOfFreedom());
758  if (((*it)->chi2() / dof) > 20000.0)
759  continue;
760 
761  // Global position of starting point for protoTrack
762  GlobalPoint gp = (*it)->globalPosition();
763  float eta_temp = gp.eta();
764  float phi_temp = gp.phi();
765  bool showeringBefore = false;
766  NShowerSeg = 0;
767  if (IdentifyShowering(CSClist2B, usedCSClist2B, eta_temp, phi_temp, 2, NShowerSeg))
768  showeringBefore = true;
769  int NShowers = 0;
770  if (showeringBefore) {
771  // usedCSClist2B[index] = true;
772  NShowers++;
773  }
774 
775  SegmentContainer protoTrack;
776  protoTrack.push_back(*it);
777 
778  std::vector<int> layers;
779  layers.push_back(2);
780 
781  // Try adding segment from other stations
782  if (foundMatchingSegment(
783  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
784  layers.push_back(3);
785  if (showeringBefore)
786  NShowers++;
787  if (foundMatchingSegment(
788  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
789  layers.push_back(4);
790  if (showeringBefore)
791  NShowers++;
792 
793  // adding showering information
794  if (layers.size() < 2 && !ShoweringSegments.empty()) {
795  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
796  if (ShoweringLayers[i] <= layers[layers.size() - 1])
797  continue;
798  protoTrack.push_back(ShoweringSegments[i]);
799  layers.push_back(ShoweringLayers[i]);
800  }
801  }
802  ShoweringSegments.clear();
803  ShoweringLayers.clear();
804 
805  TrajectorySeed thisSeed;
806  if (layers.size() < 2) {
807  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
808  } else {
809  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
810  }
811 
812  // Add the seeds to master collection
813  rawSeeds.push_back(thisSeed);
814  etaOfSeed.push_back(eta_temp);
815  phiOfSeed.push_back(phi_temp);
816  nSegOnSeed.push_back(protoTrack.size());
817  // mark this segment as used
818  usedCSClist2B[index] = true;
819  }
820 
821  // Loop over all possible ME3 segment to form seeds:
822  index = -1;
823  for (SegmentContainer::iterator it = CSClist3B.begin(); it != CSClist3B.end(); ++it) {
824  index++;
825 
826  if (usedCSClist3B[index] == true)
827  continue;
828  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
829  continue;
830 
831  double dof = static_cast<double>((*it)->degreesOfFreedom());
832  if (((*it)->chi2() / dof) > 20000.0)
833  continue;
834 
835  // Global position of starting point for protoTrack
836  GlobalPoint gp = (*it)->globalPosition();
837  float eta_temp = gp.eta();
838  float phi_temp = gp.phi();
839  bool showeringBefore = false;
840  NShowerSeg = 0;
841  if (IdentifyShowering(CSClist3B, usedCSClist3B, eta_temp, phi_temp, 3, NShowerSeg))
842  showeringBefore = true;
843  int NShowers = 0;
844  if (showeringBefore) {
845  // usedCSClist3B[index] = true;
846  NShowers++;
847  }
848 
849  SegmentContainer protoTrack;
850  protoTrack.push_back(*it);
851 
852  std::vector<int> layers;
853  layers.push_back(2);
854 
855  // Try adding segment from other stations
856  if (foundMatchingSegment(
857  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
858  layers.push_back(4);
859  if (showeringBefore)
860  NShowers++;
861 
862  // adding showering information
863  if (layers.size() < 2 && !ShoweringSegments.empty()) {
864  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
865  if (ShoweringLayers[i] <= layers[layers.size() - 1])
866  continue;
867  protoTrack.push_back(ShoweringSegments[i]);
868  layers.push_back(ShoweringLayers[i]);
869  }
870  }
871  ShoweringSegments.clear();
872  ShoweringLayers.clear();
873 
874  // mark this segment as used
875  usedCSClist3B[index] = true;
876 
877  if (layers.size() < 2)
878  continue;
879  TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
880 
881  // Add the seeds to master collection
882  rawSeeds.push_back(thisSeed);
883  etaOfSeed.push_back(eta_temp);
884  phiOfSeed.push_back(phi_temp);
885  nSegOnSeed.push_back(protoTrack.size());
886  }
887 
888  /* *****************************************************************************************************************
889  * Form seeds from forward endcap
890  *
891  * Proceed from inside -> out
892  * *****************************************************************************************************************/
893 
894  // Loop over all possible ME11 segment to form seeds:
895  index = -1;
896  for (SegmentContainer::iterator it = CSClist0F.begin(); it != CSClist0F.end(); ++it) {
897  index++;
898 
899  if (usedCSClist0F[index] == true)
900  continue;
901  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
902  continue;
903 
904  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
905  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
906 
907  // Global position of starting point for protoTrack
908  GlobalPoint gp = (*it)->globalPosition();
909  float eta_temp = gp.eta();
910  float phi_temp = gp.phi();
911  bool showeringBefore = false;
912  NShowerSeg = 0;
913  if (IdentifyShowering(CSClist0F, usedCSClist0F, eta_temp, phi_temp, 0, NShowerSeg))
914  showeringBefore = true;
915  int NShowers = 0;
916  if (showeringBefore) {
917  // usedCSClist0F[index] = true;
918  NShowers++;
919  }
920 
921  SegmentContainer protoTrack;
922  protoTrack.push_back(*it);
923 
924  std::vector<int> layers;
925  layers.push_back(0);
926 
927  // Try adding segment from other station
928  if (foundMatchingSegment(
929  1, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
930  layers.push_back(1);
931  if (showeringBefore)
932  NShowers++;
933  if (foundMatchingSegment(
934  1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
935  layers.push_back(2);
936  if (showeringBefore)
937  NShowers++;
938  if (foundMatchingSegment(
939  1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
940  layers.push_back(3);
941  if (showeringBefore)
942  NShowers++;
943  if (foundMatchingSegment(
944  1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
945  layers.push_back(4);
946  if (showeringBefore)
947  NShowers++;
948 
949  // adding showering information
950  if (layers.size() < 2 && !ShoweringSegments.empty()) {
951  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
952  if (ShoweringLayers[i] <= layers[layers.size() - 1])
953  continue;
954  protoTrack.push_back(ShoweringSegments[i]);
955  layers.push_back(ShoweringLayers[i]);
956  }
957  }
958  ShoweringSegments.clear();
959  ShoweringLayers.clear();
960 
961  TrajectorySeed thisSeed;
962  if (layers.size() < 2) {
963  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
964  } else {
965  if (fabs(gp.eta()) > 1.7) {
966  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
967  } else {
968  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
969  }
970  }
971  // Add the seeds to master collection
972  rawSeeds.push_back(thisSeed);
973  etaOfSeed.push_back(eta_temp);
974  phiOfSeed.push_back(phi_temp);
975  nSegOnSeed.push_back(protoTrack.size());
976 
977  // mark this segment as used
978  usedCSClist0F[index] = true;
979  }
980 
981  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
982  index = -1;
983  for (SegmentContainer::iterator it = CSClist1F.begin(); it != CSClist1F.end(); ++it) {
984  index++;
985 
986  if (usedCSClist1F[index] == true)
987  continue;
988  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
989  continue;
990 
991  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
992  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
993 
994  // Global position of starting point for protoTrack
995  GlobalPoint gp = (*it)->globalPosition();
996  float eta_temp = gp.eta();
997  float phi_temp = gp.phi();
998  bool showeringBefore = false;
999  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 
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());
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 }
CSCRecHit2DCollection.h
MuonSeedBuilder::etaError
double etaError(const GlobalPoint gp, double rErr)
calculate the eta error from global R error
Definition: MuonSeedBuilder.cc:1457
Handle.h
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:1235
MuonSeedBuilder::ShoweringSegments
SegmentContainer ShoweringSegments
Definition: MuonSeedBuilder.h:130
DOFs::dof
dof
Definition: AlignPCLThresholdsWriter.cc:37
mps_fire.i
i
Definition: mps_fire.py:355
MuonSeedBuilder::NShowerSeg
int NShowerSeg
Definition: MuonSeedBuilder.h:129
MessageLogger.h
MuonSeedBuilder::maxDeltaEtaCSC
float maxDeltaEtaCSC
Definition: MuonSeedBuilder.h:121
MuonSeedBuilder::enableDTMeasurement
bool enableDTMeasurement
Definition: MuonSeedBuilder.h:109
TrajectorySeedCollection
std::vector< TrajectorySeed > TrajectorySeedCollection
Definition: TrajectorySeedCollection.h:6
MuonSeedBuilder::maxDeltaPhiOverlap
float maxDeltaPhiOverlap
Definition: MuonSeedBuilder.h:124
ESHandle.h
MuonSeedBuilder::ShoweringLayers
std::vector< int > ShoweringLayers
Definition: MuonSeedBuilder.h:131
MuonDetLayerMeasurements
Definition: MuonDetLayerMeasurements.h:41
gather_cfg.cout
cout
Definition: gather_cfg.py:144
MuonSeedCleaner.h
MuonDetLayerGeometry::allDTLayers
const std::vector< const DetLayer * > & allDTLayers() const
return the DT DetLayers (barrel), inside-out
Definition: MuonDetLayerGeometry.cc:150
MuonSeedBuilder::enableCSCMeasurement
bool enableCSCMeasurement
Definition: MuonSeedBuilder.h:112
MuonSeedBuilder::BoolContainer
std::deque< bool > BoolContainer
Definition: MuonSeedBuilder.h:33
MuonSeedBuilder::SegmentContainer
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
Definition: MuonSeedBuilder.h:32
MuonDetLayerGeometry.h
MuonSeedBuilder::muonSeedClean_
MuonSeedCleaner * muonSeedClean_
Definition: MuonSeedBuilder.h:140
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
MuonSeedBuilder::maxDeltaEtaDT
float maxDeltaEtaDT
Definition: MuonSeedBuilder.h:125
case1
#define case1(par)
Definition: vmac.h:83
MuonSeedBuilder::theCSCSegmentLabel
edm::InputTag theCSCSegmentLabel
Name of the CSC segment collection.
Definition: MuonSeedBuilder.h:136
MuonDetLayerGeometry::backwardCSCLayers
const std::vector< const DetLayer * > & backwardCSCLayers() const
return the backward (-Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:156
MuonSeedBuilder::maxDeltaPhiCSC
float maxDeltaPhiCSC
Definition: MuonSeedBuilder.h:122
MuonSeedBuilder::maxEtaResolutionDT
float maxEtaResolutionDT
Definition: MuonSeedBuilder.h:150
DetId
Definition: DetId.h:17
case2
#define case2(par1, par2)
Definition: vmac.h:94
MuonSeedBuilder::maxDeltaPhiDT
float maxDeltaPhiDT
Definition: MuonSeedBuilder.h:126
PTrajectoryStateOnDet.h
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
MuonSeedCleaner::NRecHitsFromSegment
int NRecHitsFromSegment(const TrackingRecHit &rhit)
Definition: MuonSeedCleaner.cc:451
MuonSeedBuilder::maxPhiResolutionCSC
float maxPhiResolutionCSC
Definition: MuonSeedBuilder.h:153
MuonSubdetId::DT
static constexpr int DT
Definition: MuonSubdetId.h:11
Point3DBase< float, GlobalTag >
MuonDetLayerGeometry::forwardCSCLayers
const std::vector< const DetLayer * > & forwardCSCLayers() const
return the forward (+Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:154
TrajectorySeed.h
MuonSeedBuilder::muonMeasurements
MuonDetLayerMeasurements * muonMeasurements
Definition: MuonSeedBuilder.h:156
MuonSeedBuilder::theDTSegmentLabel
edm::InputTag theDTSegmentLabel
Name of the DT segment collection.
Definition: MuonSeedBuilder.h:133
MuonSeedBuilder::~MuonSeedBuilder
~MuonSeedBuilder()
Destructor.
Definition: MuonSeedBuilder.cc:104
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum)
Definition: DetId.h:48
MuonSeedBuilder::debug
bool debug
group the seeds
Definition: MuonSeedBuilder.h:106
runTauDisplay.gp
gp
Definition: runTauDisplay.py:431
MuonSeedBuilder::muonLayers
const MuonDetLayerGeometry * muonLayers
Definition: MuonSeedBuilder.h:143
funct::true
true
Definition: Factorize.h:173
edm::ParameterSet
Definition: ParameterSet.h:36
Event.h
MuonSeedBuilder.h
MuonSeedBuilder::MuonSeedBuilder
MuonSeedBuilder(const edm::ParameterSet &, edm::ConsumesCollector &)
Constructor.
Definition: MuonSeedBuilder.cc:51
MuonSeedBuilder::minDTHitsPerSegment
int minDTHitsPerSegment
Definition: MuonSeedBuilder.h:118
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
MuonSeedCleaner
Definition: MuonSeedCleaner.h:28
MuonRecoGeometryRecord.h
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:1380
edm::EventSetup
Definition: EventSetup.h:57
DTRecSegment4D.h
MuonSeedBuilder::maxEtaResolutionCSC
float maxEtaResolutionCSC
Definition: MuonSeedBuilder.h:151
MuonSeedBuilder::minCSCHitsPerSegment
int minCSCHitsPerSegment
Definition: MuonSeedBuilder.h:115
MuonSeedBuilder::maxDeltaEtaOverlap
float maxDeltaEtaOverlap
Definition: MuonSeedBuilder.h:123
TrajectorySeedCollection.h
MuonSeedBuilder::BField
const MagneticField * BField
Definition: MuonSeedBuilder.h:146
MuonDetLayerMeasurements::recHits
MuonRecHitContainer recHits(const DetLayer *layer, const edm::Event &iEvent)
returns the rechits which are on the layer
Definition: MuonDetLayerMeasurements.cc:417
MuonSeedCreator.h
MuonSeedCreator::setBField
void setBField(const MagneticField *theField)
Cache Magnetic Field for current event.
Definition: MuonSeedCreator.h:43
GeomDet.h
type
type
Definition: HCALResponse.h:21
MuonSeedCreator
Definition: MuonSeedCreator.h:30
MuonSeedBuilder::maxPhiResolutionDT
float maxPhiResolutionDT
Definition: MuonSeedBuilder.h:152
MuonDetLayerMeasurements.h
MuonSeedCleaner::seedCleaner
std::vector< TrajectorySeed > seedCleaner(const edm::EventSetup &eventSetup, std::vector< TrajectorySeed > &seeds)
Cache pointer to geometry.
Definition: MuonSeedCleaner.cc:71
DetLayer.h
TrajectorySeed
Definition: TrajectorySeed.h:17
TrajectoryStateTransform.h
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
CSCSegment.h
CSCRecHit2D.h
MuonSeedBuilder::build
int build(edm::Event &event, const edm::EventSetup &eventSetup, TrajectorySeedCollection &seeds)
Build seed collection.
Definition: MuonSeedBuilder.cc:119
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
ParameterSet.h
event
Definition: event.py:1
edm::Event
Definition: Event.h:73
GeomDetEnumerators::isDT
bool isDT(GeomDetEnumerators::SubDetector m)
Definition: GeomDetEnumerators.cc:86
MuonSeedBuilder::muonSeedCreate_
MuonSeedCreator * muonSeedCreate_
Create seed according to region (CSC, DT, Overlap)
Definition: MuonSeedBuilder.h:139
edm::InputTag
Definition: InputTag.h:15
edm::ConsumesCollector
Definition: ConsumesCollector.h:39
hgcalTopologyTester_cfi.layers
layers
Definition: hgcalTopologyTester_cfi.py:8
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
cuy.dh
dh
Definition: cuy.py:355
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:145
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
DTRecSegment4DCollection.h
CSCSegmentCollection.h
Generated for CMSSW Reference Manual by doxygen 1.8.16