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  // muon service
84  edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
85  theService = new MuonServiceProxy(serviceParameters);
86 
87  // Class for forming muon seeds
88  muonSeedCreate_ = new MuonSeedCreator(pset);
89  muonSeedClean_ = new MuonSeedCleaner(pset);
90 
91  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
92  muonMeasurements = new MuonDetLayerMeasurements(theDTSegmentLabel,
93  theCSCSegmentLabel,
94  edm::InputTag(),
95  edm::InputTag(),
96  edm::InputTag(),
97  iC,
98  enableDTMeasurement,
99  enableCSCMeasurement,
100  false,
101  false,
102  false);
103 }
104 
105 /*
106  * Destructor
107  */
108 MuonSeedBuilder::~MuonSeedBuilder() {
109  delete muonSeedCreate_;
110  delete muonSeedClean_;
111  if (theService)
112  delete theService;
113  if (muonMeasurements)
114  delete muonMeasurements;
115 }
116 
117 /*
118  * build
119  *
120  * Where the segments are sorted out to make a protoTrack (vector of matching segments in different
121  * stations/layers). The protoTrack is then passed to the seed creator to create CSC, DT or overlap seeds
122  *
123  */
124 //void MuonSeedBuilder::build( MuonDetLayerMeasurements muonMeasurements, TrajectorySeedCollection& theSeeds ) {
125 int MuonSeedBuilder::build(edm::Event& event, const edm::EventSetup& eventSetup, TrajectorySeedCollection& theSeeds) {
126  // Pass the Magnetic Field to where the seed is create
127  muonSeedCreate_->setBField(BField);
128 
129  // Create temporary collection of seeds which will be cleaned up to remove combinatorics
130  std::vector<TrajectorySeed> rawSeeds;
131  std::vector<float> etaOfSeed;
132  std::vector<float> phiOfSeed;
133  std::vector<int> nSegOnSeed;
134 
135  // Instantiate the accessor (get the segments: DT + CSC but not RPC=false)
136  // MuonDetLayerMeasurements muonMeasurements(enableDTMeasurement,enableCSCMeasurement,false,
137  // theDTSegmentLabel.label(),theCSCSegmentLabel.label());
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()
157  << std::endl;
158  std::cout << "In MB1: " << DTlist1.size() << std::endl;
159  std::cout << "In MB2: " << DTlist2.size() << std::endl;
160  std::cout << "In MB3: " << DTlist3.size() << std::endl;
161  std::cout << "In MB4: " << DTlist4.size() << std::endl;
162  }
163 
164  // 1b. get the CSC segments by stations (layers):
165  // 1b.1 Global z < 0
166  std::vector<const DetLayer*> cscBackwardLayers = muonLayers->backwardCSCLayers();
167  SegmentContainer CSClist4B = muonMeasurements->recHits(cscBackwardLayers[4], event);
168  SegmentContainer CSClist3B = muonMeasurements->recHits(cscBackwardLayers[3], event);
169  SegmentContainer CSClist2B = muonMeasurements->recHits(cscBackwardLayers[2], event);
170  SegmentContainer CSClist1B = muonMeasurements->recHits(cscBackwardLayers[1], event); // ME1/2 and 1/3
171  SegmentContainer CSClist0B = muonMeasurements->recHits(cscBackwardLayers[0], event); // ME11
172 
173  BoolContainer usedCSClist4B(CSClist4B.size(), false);
174  BoolContainer usedCSClist3B(CSClist3B.size(), false);
175  BoolContainer usedCSClist2B(CSClist2B.size(), false);
176  BoolContainer usedCSClist1B(CSClist1B.size(), false);
177  BoolContainer usedCSClist0B(CSClist0B.size(), false);
178 
179  // 1b.2 Global z > 0
180  std::vector<const DetLayer*> cscForwardLayers = muonLayers->forwardCSCLayers();
181  SegmentContainer CSClist4F = muonMeasurements->recHits(cscForwardLayers[4], event);
182  SegmentContainer CSClist3F = muonMeasurements->recHits(cscForwardLayers[3], event);
183  SegmentContainer CSClist2F = muonMeasurements->recHits(cscForwardLayers[2], event);
184  SegmentContainer CSClist1F = muonMeasurements->recHits(cscForwardLayers[1], event);
185  SegmentContainer CSClist0F = muonMeasurements->recHits(cscForwardLayers[0], event);
186 
187  BoolContainer usedCSClist4F(CSClist4F.size(), false);
188  BoolContainer usedCSClist3F(CSClist3F.size(), false);
189  BoolContainer usedCSClist2F(CSClist2F.size(), false);
190  BoolContainer usedCSClist1F(CSClist1F.size(), false);
191  BoolContainer usedCSClist0F(CSClist0F.size(), false);
192 
193  if (debug) {
194  std::cout << "*** Number of CSC segments is "
195  << CSClist4F.size() + CSClist3F.size() + CSClist2F.size() + CSClist1F.size() + CSClist0F.size() +
196  CSClist4B.size() + CSClist3B.size() + CSClist2B.size() + CSClist1B.size() + CSClist0B.size()
197  << std::endl;
198  std::cout << "In ME11: " << CSClist0F.size() + CSClist0B.size() << std::endl;
199  std::cout << "In ME12: " << CSClist1F.size() + CSClist1B.size() << std::endl;
200  std::cout << "In ME2 : " << CSClist2F.size() + CSClist2B.size() << std::endl;
201  std::cout << "In ME3 : " << CSClist3F.size() + CSClist3B.size() << std::endl;
202  std::cout << "In ME4 : " << CSClist4F.size() + CSClist4B.size() << std::endl;
203  }
204 
205  /* ******************************************************************************************************************
206  * Form seeds in barrel region
207  *
208  * Proceed from inside -> out
209  * ******************************************************************************************************************/
210 
211  // Loop over all possible MB1 segment to form seeds:
212  int index = -1;
213  for (SegmentContainer::iterator it = DTlist1.begin(); it != DTlist1.end(); ++it) {
214  index++;
215 
216  if (usedDTlist1[index] == true)
217  continue;
218  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
219  continue;
220  if ((*it)->dimension() != 4)
221  continue;
222 
223  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
224  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
225 
226  // Global position of starting point for protoTrack
227  GlobalPoint gp = (*it)->globalPosition();
228  float eta_temp = gp.eta();
229  float phi_temp = gp.phi();
230  bool showeringBefore = false;
231  NShowerSeg = 0;
232  if (IdentifyShowering(DTlist1, usedDTlist1, eta_temp, phi_temp, -1, NShowerSeg))
233  showeringBefore = true;
234  int NShowers = 0;
235  if (showeringBefore) {
236  //usedDTlist1[index] = true;
237  NShowers++;
238  }
239 
240  SegmentContainer protoTrack;
241  protoTrack.push_back(*it);
242 
243  std::vector<int> layers;
244  layers.push_back(-1);
245 
246  // Try adding segment from other stations
247  if (foundMatchingSegment(
248  3, protoTrack, DTlist2, usedDTlist2, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
249  layers.push_back(-2);
250  if (showeringBefore)
251  NShowers++;
252  if (foundMatchingSegment(
253  3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
254  layers.push_back(-3);
255  if (showeringBefore)
256  NShowers++;
257  if (foundMatchingSegment(
258  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
259  layers.push_back(-4);
260  if (showeringBefore)
261  NShowers++;
262 
263  // Check if in overlap region
264  if (eta_temp > 0.8) {
265  if (foundMatchingSegment(
266  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
267  layers.push_back(1);
268  if (showeringBefore)
269  NShowers++;
270  if (foundMatchingSegment(
271  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
272  layers.push_back(2);
273  if (showeringBefore)
274  NShowers++;
275  if (foundMatchingSegment(
276  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
277  layers.push_back(3);
278  if (showeringBefore)
279  NShowers++;
280  } else if (eta_temp < -0.8) {
281  if (foundMatchingSegment(
282  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
283  layers.push_back(1);
284  if (showeringBefore)
285  NShowers++;
286  if (foundMatchingSegment(
287  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
288  layers.push_back(2);
289  if (showeringBefore)
290  NShowers++;
291  if (foundMatchingSegment(
292  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
293  layers.push_back(3);
294  if (showeringBefore)
295  NShowers++;
296  }
297 
298  // adding showering information
299  if (layers.size() < 2 && !ShoweringSegments.empty()) {
300  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
301  if (ShoweringLayers[i] > 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  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
308  if (ShoweringLayers[i] >= layers[layers.size() - 1])
309  continue;
310  protoTrack.push_back(ShoweringSegments[i]);
311  layers.push_back(ShoweringLayers[i]);
312  }
313  }
314  }
315  ShoweringSegments.clear();
316  ShoweringLayers.clear();
317 
318  TrajectorySeed thisSeed;
319 
320  if (layers.size() < 2) {
321  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
322  } else {
323  if (layers[layers.size() - 1] > 0) {
324  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
325  } else {
326  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
327  }
328  }
329  // Add the seeds to master collection
330  rawSeeds.push_back(thisSeed);
331  etaOfSeed.push_back(eta_temp);
332  phiOfSeed.push_back(phi_temp);
333  nSegOnSeed.push_back(protoTrack.size());
334 
335  // Marked segment as used
336  usedDTlist1[index] = true;
337  }
338 
339  // Loop over all possible MB2 segment to form seeds:
340  index = -1;
341  for (SegmentContainer::iterator it = DTlist2.begin(); it != DTlist2.end(); ++it) {
342  index++;
343 
344  if (usedDTlist2[index] == true)
345  continue;
346  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
347  continue;
348  if ((*it)->dimension() != 4)
349  continue;
350 
351  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
352  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
353 
354  // Global position of starting point for protoTrack
355  GlobalPoint gp = (*it)->globalPosition();
356  float eta_temp = gp.eta();
357  float phi_temp = gp.phi();
358  bool showeringBefore = false;
359  NShowerSeg = 0;
360  if (IdentifyShowering(DTlist2, usedDTlist2, eta_temp, phi_temp, -2, NShowerSeg))
361  showeringBefore = true;
362  int NShowers = 0;
363  if (showeringBefore) {
364  // usedDTlist2[index] = true;
365  NShowers++;
366  }
367 
368  SegmentContainer protoTrack;
369  protoTrack.push_back(*it);
370 
371  std::vector<int> layers;
372  layers.push_back(-2);
373 
374  // Try adding segment from other stations
375  if (foundMatchingSegment(
376  3, protoTrack, DTlist3, usedDTlist3, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
377  layers.push_back(-3);
378  if (showeringBefore)
379  NShowers++;
380  if (foundMatchingSegment(
381  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
382  layers.push_back(-4);
383  if (showeringBefore)
384  NShowers++;
385 
386  // Check if in overlap region
387  if (eta_temp > 0.8) {
388  if (foundMatchingSegment(
389  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
390  layers.push_back(1);
391  if (showeringBefore)
392  NShowers++;
393  if (foundMatchingSegment(
394  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
395  layers.push_back(2);
396  if (showeringBefore)
397  NShowers++;
398  if (foundMatchingSegment(
399  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
400  layers.push_back(3);
401  if (showeringBefore)
402  NShowers++;
403  } else if (eta_temp < -0.8) {
404  if (foundMatchingSegment(
405  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
406  layers.push_back(1);
407  if (showeringBefore)
408  NShowers++;
409  if (foundMatchingSegment(
410  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
411  layers.push_back(2);
412  if (showeringBefore)
413  NShowers++;
414  if (foundMatchingSegment(
415  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
416  layers.push_back(3);
417  if (showeringBefore)
418  NShowers++;
419  }
420 
421  // adding showering information
422  if (layers.size() < 2 && !ShoweringSegments.empty()) {
423  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
424  if (ShoweringLayers[i] > 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  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
431  if (ShoweringLayers[i] >= layers[layers.size() - 1])
432  continue;
433  protoTrack.push_back(ShoweringSegments[i]);
434  layers.push_back(ShoweringLayers[i]);
435  }
436  }
437  }
438  ShoweringSegments.clear();
439  ShoweringLayers.clear();
440 
441  TrajectorySeed thisSeed;
442 
443  if (layers.size() < 2) {
444  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
445  } else {
446  if (layers[layers.size() - 1] > 0) {
447  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
448  } else {
449  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
450  }
451  }
452  // Add the seeds to master collection
453  rawSeeds.push_back(thisSeed);
454  etaOfSeed.push_back(eta_temp);
455  phiOfSeed.push_back(phi_temp);
456  nSegOnSeed.push_back(protoTrack.size());
457 
458  // Marked segment as used
459  usedDTlist2[index] = true;
460  }
461 
462  // Loop over all possible MB3 segment to form seeds:
463  index = -1;
464  for (SegmentContainer::iterator it = DTlist3.begin(); it != DTlist3.end(); ++it) {
465  index++;
466 
467  if (usedDTlist3[index] == true)
468  continue;
469  if (int((*it)->recHits().size()) < minDTHitsPerSegment)
470  continue;
471  if ((*it)->dimension() != 4)
472  continue;
473 
474  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
475  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
476 
477  // Global position of starting point for protoTrack
478  GlobalPoint gp = (*it)->globalPosition();
479  float eta_temp = gp.eta();
480  float phi_temp = gp.phi();
481  bool showeringBefore = false;
482  NShowerSeg = 0;
483  if (IdentifyShowering(DTlist3, usedDTlist3, eta_temp, phi_temp, -3, NShowerSeg))
484  showeringBefore = true;
485  int NShowers = 0;
486  if (showeringBefore) {
487  // usedDTlist3[index] = true;
488  NShowers++;
489  }
490 
491  SegmentContainer protoTrack;
492  protoTrack.push_back(*it);
493 
494  std::vector<int> layers;
495  layers.push_back(-3);
496 
497  // Try adding segment from other stations
498  if (foundMatchingSegment(
499  3, protoTrack, DTlist4, usedDTlist4, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
500  layers.push_back(-4);
501  if (showeringBefore)
502  NShowers++;
503 
504  // Check if in overlap region
505  if (eta_temp > 0.8) {
506  if (foundMatchingSegment(
507  2, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
508  layers.push_back(1);
509  if (showeringBefore)
510  NShowers++;
511  if (foundMatchingSegment(
512  2, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
513  layers.push_back(2);
514  if (showeringBefore)
515  NShowers++;
516  if (foundMatchingSegment(
517  2, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
518  layers.push_back(3);
519  if (showeringBefore)
520  NShowers++;
521  } else if (eta_temp < -0.8) {
522  if (foundMatchingSegment(
523  2, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
524  layers.push_back(1);
525  if (showeringBefore)
526  NShowers++;
527  if (foundMatchingSegment(
528  2, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
529  layers.push_back(2);
530  if (showeringBefore)
531  NShowers++;
532  if (foundMatchingSegment(
533  2, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
534  layers.push_back(3);
535  if (showeringBefore)
536  NShowers++;
537  }
538 
539  // adding showering information
540  if (layers.size() < 2 && !ShoweringSegments.empty()) {
541  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
542  if (ShoweringLayers[i] > 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  if (ShoweringLayers[i] < 0 && layers[layers.size() - 1] < 0) {
549  if (ShoweringLayers[i] >= layers[layers.size() - 1])
550  continue;
551  protoTrack.push_back(ShoweringSegments[i]);
552  layers.push_back(ShoweringLayers[i]);
553  }
554  }
555  }
556  ShoweringSegments.clear();
557  ShoweringLayers.clear();
558 
559  TrajectorySeed thisSeed;
560  if (layers.size() < 2) {
561  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
562  } else {
563  if (layers[layers.size() - 1] > 0) {
564  thisSeed = muonSeedCreate_->createSeed(2, protoTrack, layers, NShowers, NShowerSeg);
565  } else {
566  thisSeed = muonSeedCreate_->createSeed(3, protoTrack, layers, NShowers, NShowerSeg);
567  }
568  }
569 
570  // Add the seeds to master collection
571  rawSeeds.push_back(thisSeed);
572  etaOfSeed.push_back(eta_temp);
573  phiOfSeed.push_back(phi_temp);
574  nSegOnSeed.push_back(protoTrack.size());
575 
576  // Marked segment as used
577  usedDTlist3[index] = true;
578  }
579 
580  /* *********************************************************************************************************************
581  * Form seeds from backward endcap
582  *
583  * Proceed from inside -> out
584  * *********************************************************************************************************************/
585 
586  // Loop over all possible ME11 segment to form seeds:
587  index = -1;
588 
589  for (SegmentContainer::iterator it = CSClist0B.begin(); it != CSClist0B.end(); ++it) {
590  index++;
591 
592  if (usedCSClist0B[index] == true)
593  continue;
594  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
595  continue;
596 
597  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
598  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
599 
600  // Global position of starting point for protoTrack
601  GlobalPoint gp = (*it)->globalPosition();
602  float eta_temp = gp.eta();
603  float phi_temp = gp.phi();
604  bool showeringBefore = false;
605  NShowerSeg = 0;
606  if (IdentifyShowering(CSClist0B, usedCSClist0B, eta_temp, phi_temp, 0, NShowerSeg))
607  showeringBefore = true;
608  int NShowers = 0;
609  if (showeringBefore) {
610  // usedCSClist0B[index] = true;
611  NShowers++;
612  }
613 
614  SegmentContainer protoTrack;
615  protoTrack.push_back(*it);
616 
617  std::vector<int> layers;
618  layers.push_back(0);
619 
620  // Try adding segment from other station
621  if (foundMatchingSegment(
622  1, protoTrack, CSClist1B, usedCSClist1B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
623  layers.push_back(1);
624  if (showeringBefore)
625  NShowers++;
626  if (foundMatchingSegment(
627  1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
628  layers.push_back(2);
629  if (showeringBefore)
630  NShowers++;
631  if (foundMatchingSegment(
632  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
633  layers.push_back(3);
634  if (showeringBefore)
635  NShowers++;
636  if (foundMatchingSegment(
637  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
638  layers.push_back(4);
639  if (showeringBefore)
640  NShowers++;
641 
642  // adding showering information
643  if (layers.size() < 2 && !ShoweringSegments.empty()) {
644  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
645  if (ShoweringLayers[i] <= layers[layers.size() - 1])
646  continue;
647  protoTrack.push_back(ShoweringSegments[i]);
648  layers.push_back(ShoweringLayers[i]);
649  }
650  }
651  ShoweringSegments.clear();
652  ShoweringLayers.clear();
653 
654  TrajectorySeed thisSeed;
655  if (layers.size() < 2) {
656  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
657  } else {
658  if (fabs(gp.eta()) > 1.7) {
659  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
660  } else {
661  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
662  }
663  }
664 
665  // Add the seeds to master collection
666  rawSeeds.push_back(thisSeed);
667  etaOfSeed.push_back(eta_temp);
668  phiOfSeed.push_back(phi_temp);
669  nSegOnSeed.push_back(protoTrack.size());
670 
671  // mark this segment as used
672  usedCSClist0B[index] = true;
673  }
674 
675  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
676  index = -1;
677  for (SegmentContainer::iterator it = CSClist1B.begin(); it != CSClist1B.end(); ++it) {
678  index++;
679 
680  if (usedCSClist1B[index] == true)
681  continue;
682  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
683  continue;
684 
685  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
686  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
687 
688  // Global position of starting point for protoTrack
689  GlobalPoint gp = (*it)->globalPosition();
690  float eta_temp = gp.eta();
691  float phi_temp = gp.phi();
692  bool showeringBefore = false;
693  NShowerSeg = 0;
694  if (IdentifyShowering(CSClist1B, usedCSClist1B, eta_temp, phi_temp, 1, NShowerSeg))
695  showeringBefore = true;
696  int NShowers = 0;
697  if (showeringBefore) {
698  // usedCSClist1B[index] = true;
699  NShowers++;
700  }
701 
702  SegmentContainer protoTrack;
703  protoTrack.push_back(*it);
704 
705  std::vector<int> layers;
706  layers.push_back(1);
707 
708  // Try adding segment from other stations
709  if (foundMatchingSegment(
710  1, protoTrack, CSClist2B, usedCSClist2B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
711  layers.push_back(2);
712  if (showeringBefore)
713  NShowers++;
714  if (foundMatchingSegment(
715  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
716  layers.push_back(3);
717  if (showeringBefore)
718  NShowers++;
719  if (foundMatchingSegment(
720  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
721  layers.push_back(4);
722  if (showeringBefore)
723  NShowers++;
724 
725  // adding showering information
726  if (layers.size() < 2 && !ShoweringSegments.empty()) {
727  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
728  if (ShoweringLayers[i] <= layers[layers.size() - 1])
729  continue;
730  protoTrack.push_back(ShoweringSegments[i]);
731  layers.push_back(ShoweringLayers[i]);
732  }
733  }
734  ShoweringSegments.clear();
735  ShoweringLayers.clear();
736 
737  TrajectorySeed thisSeed;
738  if (layers.size() < 2) {
739  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
740  } else {
741  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
742  }
743  // Add the seeds to master collection
744  rawSeeds.push_back(thisSeed);
745  etaOfSeed.push_back(eta_temp);
746  phiOfSeed.push_back(phi_temp);
747  nSegOnSeed.push_back(protoTrack.size());
748 
749  // mark this segment as used
750  usedCSClist1B[index] = true;
751  }
752 
753  // Loop over all possible ME2 segment to form seeds:
754  index = -1;
755  for (SegmentContainer::iterator it = CSClist2B.begin(); it != CSClist2B.end(); ++it) {
756  index++;
757 
758  if (usedCSClist2B[index] == true)
759  continue;
760  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
761  continue;
762 
763  double dof = static_cast<double>((*it)->degreesOfFreedom());
764  if (((*it)->chi2() / dof) > 20000.0)
765  continue;
766 
767  // Global position of starting point for protoTrack
768  GlobalPoint gp = (*it)->globalPosition();
769  float eta_temp = gp.eta();
770  float phi_temp = gp.phi();
771  bool showeringBefore = false;
772  NShowerSeg = 0;
773  if (IdentifyShowering(CSClist2B, usedCSClist2B, eta_temp, phi_temp, 2, NShowerSeg))
774  showeringBefore = true;
775  int NShowers = 0;
776  if (showeringBefore) {
777  // usedCSClist2B[index] = true;
778  NShowers++;
779  }
780 
781  SegmentContainer protoTrack;
782  protoTrack.push_back(*it);
783 
784  std::vector<int> layers;
785  layers.push_back(2);
786 
787  // Try adding segment from other stations
788  if (foundMatchingSegment(
789  1, protoTrack, CSClist3B, usedCSClist3B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
790  layers.push_back(3);
791  if (showeringBefore)
792  NShowers++;
793  if (foundMatchingSegment(
794  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
795  layers.push_back(4);
796  if (showeringBefore)
797  NShowers++;
798 
799  // adding showering information
800  if (layers.size() < 2 && !ShoweringSegments.empty()) {
801  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
802  if (ShoweringLayers[i] <= layers[layers.size() - 1])
803  continue;
804  protoTrack.push_back(ShoweringSegments[i]);
805  layers.push_back(ShoweringLayers[i]);
806  }
807  }
808  ShoweringSegments.clear();
809  ShoweringLayers.clear();
810 
811  TrajectorySeed thisSeed;
812  if (layers.size() < 2) {
813  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
814  } else {
815  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
816  }
817 
818  // Add the seeds to master collection
819  rawSeeds.push_back(thisSeed);
820  etaOfSeed.push_back(eta_temp);
821  phiOfSeed.push_back(phi_temp);
822  nSegOnSeed.push_back(protoTrack.size());
823  // mark this segment as used
824  usedCSClist2B[index] = true;
825  }
826 
827  // Loop over all possible ME3 segment to form seeds:
828  index = -1;
829  for (SegmentContainer::iterator it = CSClist3B.begin(); it != CSClist3B.end(); ++it) {
830  index++;
831 
832  if (usedCSClist3B[index] == true)
833  continue;
834  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
835  continue;
836 
837  double dof = static_cast<double>((*it)->degreesOfFreedom());
838  if (((*it)->chi2() / dof) > 20000.0)
839  continue;
840 
841  // Global position of starting point for protoTrack
842  GlobalPoint gp = (*it)->globalPosition();
843  float eta_temp = gp.eta();
844  float phi_temp = gp.phi();
845  bool showeringBefore = false;
846  NShowerSeg = 0;
847  if (IdentifyShowering(CSClist3B, usedCSClist3B, eta_temp, phi_temp, 3, NShowerSeg))
848  showeringBefore = true;
849  int NShowers = 0;
850  if (showeringBefore) {
851  // usedCSClist3B[index] = true;
852  NShowers++;
853  }
854 
855  SegmentContainer protoTrack;
856  protoTrack.push_back(*it);
857 
858  std::vector<int> layers;
859  layers.push_back(2);
860 
861  // Try adding segment from other stations
862  if (foundMatchingSegment(
863  1, protoTrack, CSClist4B, usedCSClist4B, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
864  layers.push_back(4);
865  if (showeringBefore)
866  NShowers++;
867 
868  // adding showering information
869  if (layers.size() < 2 && !ShoweringSegments.empty()) {
870  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
871  if (ShoweringLayers[i] <= layers[layers.size() - 1])
872  continue;
873  protoTrack.push_back(ShoweringSegments[i]);
874  layers.push_back(ShoweringLayers[i]);
875  }
876  }
877  ShoweringSegments.clear();
878  ShoweringLayers.clear();
879 
880  // mark this segment as used
881  usedCSClist3B[index] = true;
882 
883  if (layers.size() < 2)
884  continue;
885  TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
886 
887  // Add the seeds to master collection
888  rawSeeds.push_back(thisSeed);
889  etaOfSeed.push_back(eta_temp);
890  phiOfSeed.push_back(phi_temp);
891  nSegOnSeed.push_back(protoTrack.size());
892  }
893 
894  /* *****************************************************************************************************************
895  * Form seeds from forward endcap
896  *
897  * Proceed from inside -> out
898  * *****************************************************************************************************************/
899 
900  // Loop over all possible ME11 segment to form seeds:
901  index = -1;
902  for (SegmentContainer::iterator it = CSClist0F.begin(); it != CSClist0F.end(); ++it) {
903  index++;
904 
905  if (usedCSClist0F[index] == true)
906  continue;
907  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
908  continue;
909 
910  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
911  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
912 
913  // Global position of starting point for protoTrack
914  GlobalPoint gp = (*it)->globalPosition();
915  float eta_temp = gp.eta();
916  float phi_temp = gp.phi();
917  bool showeringBefore = false;
918  NShowerSeg = 0;
919  if (IdentifyShowering(CSClist0F, usedCSClist0F, eta_temp, phi_temp, 0, NShowerSeg))
920  showeringBefore = true;
921  int NShowers = 0;
922  if (showeringBefore) {
923  // usedCSClist0F[index] = true;
924  NShowers++;
925  }
926 
927  SegmentContainer protoTrack;
928  protoTrack.push_back(*it);
929 
930  std::vector<int> layers;
931  layers.push_back(0);
932 
933  // Try adding segment from other station
934  if (foundMatchingSegment(
935  1, protoTrack, CSClist1F, usedCSClist1F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
936  layers.push_back(1);
937  if (showeringBefore)
938  NShowers++;
939  if (foundMatchingSegment(
940  1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
941  layers.push_back(2);
942  if (showeringBefore)
943  NShowers++;
944  if (foundMatchingSegment(
945  1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
946  layers.push_back(3);
947  if (showeringBefore)
948  NShowers++;
949  if (foundMatchingSegment(
950  1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
951  layers.push_back(4);
952  if (showeringBefore)
953  NShowers++;
954 
955  // adding showering information
956  if (layers.size() < 2 && !ShoweringSegments.empty()) {
957  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
958  if (ShoweringLayers[i] <= layers[layers.size() - 1])
959  continue;
960  protoTrack.push_back(ShoweringSegments[i]);
961  layers.push_back(ShoweringLayers[i]);
962  }
963  }
964  ShoweringSegments.clear();
965  ShoweringLayers.clear();
966 
967  TrajectorySeed thisSeed;
968  if (layers.size() < 2) {
969  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
970  } else {
971  if (fabs(gp.eta()) > 1.7) {
972  thisSeed = muonSeedCreate_->createSeed(5, protoTrack, layers, NShowers, NShowerSeg);
973  } else {
974  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
975  }
976  }
977  // Add the seeds to master collection
978  rawSeeds.push_back(thisSeed);
979  etaOfSeed.push_back(eta_temp);
980  phiOfSeed.push_back(phi_temp);
981  nSegOnSeed.push_back(protoTrack.size());
982 
983  // mark this segment as used
984  usedCSClist0F[index] = true;
985  }
986 
987  // Loop over all possible ME1/2 ME1/3 segment to form seeds:
988  index = -1;
989  for (SegmentContainer::iterator it = CSClist1F.begin(); it != CSClist1F.end(); ++it) {
990  index++;
991 
992  if (usedCSClist1F[index] == true)
993  continue;
994  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
995  continue;
996 
997  //double dof = static_cast<double>( (*it)->degreesOfFreedom() ) ;
998  //if ( ((*it)->chi2()/dof) > 20000.0 ) continue;
999 
1000  // Global position of starting point for protoTrack
1001  GlobalPoint gp = (*it)->globalPosition();
1002  float eta_temp = gp.eta();
1003  float phi_temp = gp.phi();
1004  bool showeringBefore = false;
1005  NShowerSeg = 0;
1006  if (IdentifyShowering(CSClist1F, usedCSClist1F, eta_temp, phi_temp, 1, NShowerSeg))
1007  showeringBefore = true;
1008  int NShowers = 0;
1009  if (showeringBefore) {
1010  // usedCSClist1F[index] = true;
1011  NShowers++;
1012  }
1013 
1014  SegmentContainer protoTrack;
1015  protoTrack.push_back(*it);
1016 
1017  std::vector<int> layers;
1018  layers.push_back(1);
1019 
1020  // Try adding segment from other stations
1021  if (foundMatchingSegment(
1022  1, protoTrack, CSClist2F, usedCSClist2F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1023  layers.push_back(2);
1024  if (showeringBefore)
1025  NShowers++;
1026  if (foundMatchingSegment(
1027  1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1028  layers.push_back(3);
1029  if (showeringBefore)
1030  NShowers++;
1031  if (foundMatchingSegment(
1032  1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1033  layers.push_back(4);
1034  if (showeringBefore)
1035  NShowers++;
1036 
1037  // adding showering information
1038  if (layers.size() < 2 && !ShoweringSegments.empty()) {
1039  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1040  if (ShoweringLayers[i] <= layers[layers.size() - 1])
1041  continue;
1042  protoTrack.push_back(ShoweringSegments[i]);
1043  layers.push_back(ShoweringLayers[i]);
1044  }
1045  }
1046  ShoweringSegments.clear();
1047  ShoweringLayers.clear();
1048 
1049  TrajectorySeed thisSeed;
1050  if (layers.size() < 2) {
1051  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
1052  } else {
1053  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1054  }
1055 
1056  // Add the seeds to master collection
1057  rawSeeds.push_back(thisSeed);
1058  etaOfSeed.push_back(eta_temp);
1059  phiOfSeed.push_back(phi_temp);
1060  nSegOnSeed.push_back(protoTrack.size());
1061 
1062  // mark this segment as used
1063  usedCSClist1F[index] = true;
1064  }
1065 
1066  // Loop over all possible ME2 segment to form seeds:
1067  index = -1;
1068  for (SegmentContainer::iterator it = CSClist2F.begin(); it != CSClist2F.end(); ++it) {
1069  index++;
1070 
1071  if (usedCSClist2F[index] == true)
1072  continue;
1073  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
1074  continue;
1075 
1076  double dof = static_cast<double>((*it)->degreesOfFreedom());
1077  if (((*it)->chi2() / dof) > 20000.0)
1078  continue;
1079 
1080  // Global position of starting point for protoTrack
1081  GlobalPoint gp = (*it)->globalPosition();
1082  float eta_temp = gp.eta();
1083  float phi_temp = gp.phi();
1084  bool showeringBefore = false;
1085  NShowerSeg = 0;
1086  if (IdentifyShowering(CSClist2F, usedCSClist2F, eta_temp, phi_temp, 2, NShowerSeg))
1087  showeringBefore = true;
1088  int NShowers = 0;
1089  if (showeringBefore) {
1090  // usedCSClist2F[index] = true;
1091  NShowers++;
1092  }
1093 
1094  SegmentContainer protoTrack;
1095  protoTrack.push_back(*it);
1096 
1097  std::vector<int> layers;
1098  layers.push_back(2);
1099 
1100  // Try adding segment from other stations
1101  if (foundMatchingSegment(
1102  1, protoTrack, CSClist3F, usedCSClist3F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1103  layers.push_back(3);
1104  if (showeringBefore)
1105  NShowers++;
1106  if (foundMatchingSegment(
1107  1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1108  layers.push_back(4);
1109  if (showeringBefore)
1110  NShowers++;
1111 
1112  // adding showering information
1113  if (layers.size() < 2 && !ShoweringSegments.empty()) {
1114  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1115  if (ShoweringLayers[i] <= layers[layers.size() - 1])
1116  continue;
1117  protoTrack.push_back(ShoweringSegments[i]);
1118  layers.push_back(ShoweringLayers[i]);
1119  }
1120  }
1121  ShoweringSegments.clear();
1122  ShoweringLayers.clear();
1123 
1124  TrajectorySeed thisSeed;
1125  if (layers.size() < 2) {
1126  thisSeed = muonSeedCreate_->createSeed(4, protoTrack, layers, NShowers, NShowerSeg);
1127  } else {
1128  thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1129  }
1130 
1131  // Add the seeds to master collection
1132  rawSeeds.push_back(thisSeed);
1133  etaOfSeed.push_back(eta_temp);
1134  phiOfSeed.push_back(phi_temp);
1135  nSegOnSeed.push_back(protoTrack.size());
1136 
1137  // mark this segment as used
1138  usedCSClist2F[index] = true;
1139  }
1140 
1141  // Loop over all possible ME3 segment to form seeds:
1142  index = -1;
1143  for (SegmentContainer::iterator it = CSClist3F.begin(); it != CSClist3F.end(); ++it) {
1144  index++;
1145 
1146  if (usedCSClist3F[index] == true)
1147  continue;
1148  if (int((*it)->recHits().size()) < minCSCHitsPerSegment)
1149  continue;
1150 
1151  double dof = static_cast<double>((*it)->degreesOfFreedom());
1152  if (((*it)->chi2() / dof) > 20000.0)
1153  continue;
1154 
1155  // Global position of starting point for protoTrack
1156  GlobalPoint gp = (*it)->globalPosition();
1157  float eta_temp = gp.eta();
1158  float phi_temp = gp.phi();
1159  bool showeringBefore = false;
1160  NShowerSeg = 0;
1161  if (IdentifyShowering(CSClist3F, usedCSClist3F, eta_temp, phi_temp, 3, NShowerSeg))
1162  showeringBefore = true;
1163  int NShowers = 0;
1164  if (showeringBefore) {
1165  // usedCSClist3F[index] = true;
1166  NShowers++;
1167  }
1168 
1169  SegmentContainer protoTrack;
1170  protoTrack.push_back(*it);
1171 
1172  std::vector<int> layers;
1173  layers.push_back(2);
1174 
1175  // Try adding segment from other stations
1176  if (foundMatchingSegment(
1177  1, protoTrack, CSClist4F, usedCSClist4F, eta_temp, phi_temp, layers[layers.size() - 1], showeringBefore))
1178  layers.push_back(4);
1179  if (showeringBefore)
1180  NShowers++;
1181 
1182  // adding showering information
1183  if (layers.size() < 2 && !ShoweringSegments.empty()) {
1184  for (size_t i = 0; i < ShoweringSegments.size(); i++) {
1185  if (ShoweringLayers[i] <= layers[layers.size() - 1])
1186  continue;
1187  protoTrack.push_back(ShoweringSegments[i]);
1188  layers.push_back(ShoweringLayers[i]);
1189  }
1190  }
1191  ShoweringSegments.clear();
1192  ShoweringLayers.clear();
1193 
1194  // mark this segment as used
1195  usedCSClist3F[index] = true;
1196 
1197  if (layers.size() < 2)
1198  continue;
1199 
1200  TrajectorySeed thisSeed = muonSeedCreate_->createSeed(1, protoTrack, layers, NShowers, NShowerSeg);
1201 
1202  // Add the seeds to master collection
1203  rawSeeds.push_back(thisSeed);
1204  etaOfSeed.push_back(eta_temp);
1205  phiOfSeed.push_back(phi_temp);
1206  nSegOnSeed.push_back(protoTrack.size());
1207  }
1208 
1209  /* *********************************************************************************************************************
1210  * Clean up raw seed collection and pass to master collection
1211  * *********************************************************************************************************************/
1212 
1213  if (debug)
1214  std::cout << "*** CLEAN UP " << std::endl;
1215  if (debug)
1216  std::cout << "Number of seeds BEFORE " << rawSeeds.size() << std::endl;
1217 
1218  int goodSeeds = 0;
1219 
1220  theSeeds = muonSeedClean_->seedCleaner(eventSetup, rawSeeds);
1221  goodSeeds = theSeeds.size();
1222 
1223  //std::cout << " === Before Cleaning : " << rawSeeds.size() <<std::endl;
1224  //std::cout << " => AFTER :" << goodSeeds << std::endl;
1225 
1226  if (debug)
1227  std::cout << "Number of seeds AFTER " << goodSeeds << std::endl;
1228 
1229  return goodSeeds;
1230 }
1231 
1232 /* *********************************************************************************************************************
1233  * Try to match segment from different station (layer)
1234  *
1235  * Note type = 1 --> endcap
1236  * = 2 --> overlap
1237  * = 3 --> barrel
1238  * *********************************************************************************************************************/
1239 
1241 bool MuonSeedBuilder::foundMatchingSegment(int type,
1242  SegmentContainer& protoTrack,
1243  SegmentContainer& segs,
1244  BoolContainer& usedSeg,
1245  float& eta_last,
1246  float& phi_last,
1247  int& lastLayer,
1248  bool& showeringBefore) {
1249  bool ok = false;
1250  int scanlayer = (lastLayer < 0) ? (lastLayer - 1) : (lastLayer + 1);
1251 
1252  if (IdentifyShowering(segs, usedSeg, eta_last, phi_last, scanlayer, NShowerSeg)) {
1253  showeringBefore = true;
1254  return ok;
1255  }
1256 
1257  // Setup the searching cone-size
1258  // searching cone-size should changes with bending power
1259  double maxdEta;
1260  double maxdPhi;
1261  if (type == 1) {
1262  // CSC
1263  maxdEta = maxDeltaEtaCSC;
1264  if (lastLayer == 0 || lastLayer == 1) {
1265  if (fabs(eta_last) < 2.1) {
1266  maxdPhi = maxDeltaPhiCSC;
1267  } else {
1268  maxdPhi = 0.06;
1269  }
1270  } else if (lastLayer == 2) {
1271  maxdPhi = 0.5 * maxDeltaPhiCSC;
1272  } else {
1273  maxdPhi = 0.2 * maxDeltaPhiCSC;
1274  }
1275 
1276  } else if (type == 2) {
1277  // Overlap
1278  maxdEta = maxDeltaEtaOverlap;
1279  if (lastLayer == -1) {
1280  maxdPhi = maxDeltaPhiDT;
1281  } else {
1282  maxdPhi = maxDeltaPhiOverlap;
1283  }
1284 
1285  } else {
1286  // DT
1287  maxdEta = maxDeltaEtaDT;
1288  if (lastLayer == -1) {
1289  maxdPhi = maxDeltaPhiDT;
1290  } else if (lastLayer == -2) {
1291  maxdPhi = 0.8 * maxDeltaPhiDT;
1292  } else {
1293  maxdPhi = 0.4 * maxDeltaPhiDT;
1294  }
1295  }
1296 
1297  // if previous layer showers, limite the maxdPhi < 0.06
1298  if (showeringBefore && maxdPhi > 0.03)
1299  maxdPhi = 0.03;
1300  // reset the showering flag
1301  showeringBefore = false;
1302 
1303  // global phi/eta from previous segment
1304  float eta_temp = eta_last;
1305  float phi_temp = phi_last;
1306 
1307  // Index counter to keep track of element used in segs
1308  int index = -1;
1309  int best_match = index;
1310  float best_R = sqrt((maxdEta * maxdEta) + (maxdPhi * maxdPhi));
1311  float best_chi2 = 200;
1312  int best_dimension = 2;
1313  int best_nhits = minDTHitsPerSegment;
1314  if (type == 1)
1315  best_nhits = minCSCHitsPerSegment;
1316  // Loop over segments in other station (layer) and find candidate match
1317  for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1318  index++;
1319 
1320  // Not to get confused: eta_last is from the previous layer.
1321  // This is only to find the best set of segments by comparing at the distance layer-by-layer
1322  GlobalPoint gp2 = (*it)->globalPosition();
1323  double dh = fabs(gp2.eta() - eta_temp);
1324  double df = fabs(gp2.phi() - phi_temp);
1325  double dR = sqrt((dh * dh) + (df * df));
1326 
1327  // dEta and dPhi should be within certain range
1328  bool case1 = (dh < maxdEta && df < maxdPhi) ? true : false;
1329  // for DT station 4 ; CSCSegment is always 4D
1330  bool case2 = (((*it)->dimension() != 4) && (dh < 0.5) && (df < maxdPhi)) ? true : false;
1331  if (!case1 && !case2)
1332  continue;
1333 
1334  int NRechits = muonSeedClean_->NRecHitsFromSegment(&*(*it));
1335 
1336  if (NRechits < best_nhits)
1337  continue;
1338  best_nhits = NRechits;
1339 
1340  // reject 2D segments if 4D segments are available
1341  if ((*it)->dimension() < best_dimension)
1342  continue;
1343  best_dimension = (*it)->dimension();
1344 
1345  // pick the segment with best chi2/dof within a fixed cone size
1346  if (dR > best_R)
1347  continue;
1348 
1349  // select smaller chi2/dof
1350  double dof = static_cast<double>((*it)->degreesOfFreedom());
1352  if ((*it)->chi2() / dof < 0.001 && NRechits < 6 && type == 1)
1353  continue;
1354  if ((*it)->chi2() / dof > best_chi2)
1355  continue;
1356  best_chi2 = (*it)->chi2() / dof;
1357  best_match = index;
1358  // propagate the eta and phi to next layer
1359  if ((*it)->dimension() != 4) {
1360  // Self assignment, is this a bug??
1361  // should this have been (phi/eta)_last = (phi/eta)_temp to make it reset?
1362  //phi_last = phi_last;
1363  //eta_last = eta_last;
1364  } else {
1365  phi_last = gp2.phi();
1366  eta_last = gp2.eta();
1367  }
1368  }
1369 
1370  if (best_match < 0)
1371  return ok;
1372 
1373  // Add best matching segment to protoTrack:
1374  index = -1;
1375  for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1376  index++;
1377  if (index != best_match)
1378  continue;
1379  protoTrack.push_back(*it);
1380  usedSeg[best_match] = true;
1381  ok = true;
1382  }
1383  return ok;
1384 }
1385 
1386 bool MuonSeedBuilder::IdentifyShowering(SegmentContainer& segs,
1387  BoolContainer& usedSeg,
1388  float& eta_last,
1389  float& phi_last,
1390  int layer,
1391  int& NShoweringSegments) {
1392  bool showering = false;
1393 
1394  int nSeg = 0;
1395  int nRhits = 0;
1396  double nChi2 = 9999.;
1397  int theOrigin = -1;
1398  std::vector<int> badtag;
1399  int index = -1;
1400  double aveEta = 0.0;
1401  for (SegmentContainer::iterator it = segs.begin(); it != segs.end(); ++it) {
1402  index++;
1403  GlobalPoint gp = (*it)->globalPosition();
1404  double dh = gp.eta() - eta_last;
1405  double df = gp.phi() - phi_last;
1406  double dR = sqrt((dh * dh) + (df * df));
1407 
1408  double dof = static_cast<double>((*it)->degreesOfFreedom());
1409  double nX2 = (*it)->chi2() / dof;
1410 
1411  bool isDT = false;
1412  DetId geoId = (*it)->geographicalId();
1413  if (geoId.subdetId() == MuonSubdetId::DT)
1414  isDT = true;
1415 
1416  if (dR < 0.3) {
1417  nSeg++;
1418  badtag.push_back(index);
1419  aveEta += fabs(gp.eta());
1420  // pick up the best segment from showering chamber
1421  int rh = muonSeedClean_->NRecHitsFromSegment(&*(*it));
1422  if (rh < 6 && !isDT)
1423  continue;
1424  if (rh < 12 && isDT)
1425  continue;
1426  if (rh > nRhits) {
1427  nRhits = rh;
1428  if (nX2 > nChi2)
1429  continue;
1430  if (layer != 0 && layer != 1 && layer != -1) {
1431  theOrigin = index;
1432  }
1433  }
1434  }
1435  }
1436  aveEta = aveEta / static_cast<double>(nSeg);
1437  bool isME11A = (aveEta >= 2.1 && layer == 0) ? true : false;
1438  bool isME12 = (aveEta > 1.2 && aveEta <= 1.65 && layer == 1) ? true : false;
1439  bool isME11 = (aveEta > 1.65 && aveEta <= 2.1 && layer == 0) ? true : false;
1440  bool is1stLayer = (layer == -1 || layer == 0 || isME12 || isME11 || isME11A) ? true : false;
1441 
1442  NShoweringSegments += nSeg;
1443 
1444  if (nSeg > 3 && !isME11A)
1445  showering = true;
1446  if (nSeg > 6 && isME11A)
1447  showering = true;
1448 
1449  // if showering, flag all segments in order to skip this layer for pt estimation except 1st layer
1450  //std::cout<<" from Showering "<<std::endl;
1451  if (showering && !is1stLayer) {
1452  for (std::vector<int>::iterator it = badtag.begin(); it != badtag.end(); ++it) {
1453  usedSeg[*it] = true;
1454  if ((*it) != theOrigin)
1455  continue;
1456  ShoweringSegments.push_back(segs[*it]);
1457  ShoweringLayers.push_back(layer);
1458  }
1459  }
1460  return showering;
1461 }
1462 
1463 double MuonSeedBuilder::etaError(const GlobalPoint gp, double rErr) {
1464  double dHdTheta = 0.0;
1465  double dThetadR = 0.0;
1466  double etaErr = 1.0;
1467 
1468  if (gp.perp() != 0) {
1469  dHdTheta = (gp.mag() + gp.z()) / gp.perp();
1470  dThetadR = gp.z() / gp.perp2();
1471  etaErr = 0.25 * (dHdTheta * dThetadR) * (dHdTheta * dThetadR) * rErr;
1472  }
1473 
1474  return etaErr;
1475 }
MuonSeedBuilder::debug
bool debug
group the seeds
Definition: MuonSeedBuilder.h:109
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
MuonSeedBuilder::etaError
double etaError(const GlobalPoint gp, double rErr)
calculate the eta error from global R error
Definition: MuonSeedBuilder.cc:1463
case2
#define case2(par1, par2)
Definition: vmac.h:94
mps_fire.i
i
Definition: mps_fire.py:341
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
MuonSeedBuilder::SegmentContainer
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
Definition: MuonSeedBuilder.h:35
MuonSeedBuilder::ShoweringSegments
SegmentContainer ShoweringSegments
Definition: MuonSeedBuilder.h:133
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:1386
MuonSeedBuilder::build
int build(edm::Event &event, const edm::EventSetup &eventSetup, TrajectorySeedCollection &seeds)
Build seed collection.
Definition: MuonSeedBuilder.cc:125
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:1241
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
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
PV3DBase::perp2
T perp2() const
Definition: PV3DBase.h:68
GeomDetEnumerators::isDT
bool isDT(GeomDetEnumerators::SubDetector m)
Definition: GeomDetEnumerators.cc:86
MuonSeedBuilder::muonMeasurements
MuonDetLayerMeasurements * muonMeasurements
Definition: MuonSeedBuilder.h:162
MuonSeedBuilder::muonSeedClean_
MuonSeedCleaner * muonSeedClean_
Definition: MuonSeedBuilder.h:143
MuonSeedBuilder::BField
const MagneticField * BField
Definition: MuonSeedBuilder.h:149
MuonSeedBuilder::theCSCSegmentLabel
edm::InputTag theCSCSegmentLabel
Name of the CSC segment collection.
Definition: MuonSeedBuilder.h:139
MuonSeedBuilder::BoolContainer
std::deque< bool > BoolContainer
Definition: MuonSeedBuilder.h:36
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
MuonSeedBuilder::theDTSegmentLabel
edm::InputTag theDTSegmentLabel
Name of the DT segment collection.
Definition: MuonSeedBuilder.h:136
MuonSeedBuilder::ShoweringLayers
std::vector< int > ShoweringLayers
Definition: MuonSeedBuilder.h:134
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:417
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
MuonSeedCreator::setBField
void setBField(const MagneticField *theField)
Cache Magnetic Field for current event.
Definition: MuonSeedCreator.h:43
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
DetId
Definition: DetId.h:17
MuonSeedCleaner::NRecHitsFromSegment
int NRecHitsFromSegment(const TrackingRecHit &rhit)
Definition: MuonSeedCleaner.cc:452
MuonDetLayerGeometry::allDTLayers
const std::vector< const DetLayer * > & allDTLayers() const
return the DT DetLayers (barrel), inside-out
Definition: MuonDetLayerGeometry.cc:150
MuonServiceProxy_cff.MuonServiceProxy
MuonServiceProxy
Definition: MuonServiceProxy_cff.py:15
case1
#define case1(par)
Definition: vmac.h:83
MuonDetLayerGeometry::forwardCSCLayers
const std::vector< const DetLayer * > & forwardCSCLayers() const
return the forward (+Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:154
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
MuonSeedBuilder::MuonSeedBuilder
MuonSeedBuilder(const edm::ParameterSet &, edm::ConsumesCollector &)
Constructor.
Definition: MuonSeedBuilder.cc:51
MuonSeedBuilder::theService
MuonServiceProxy * theService
Definition: MuonSeedBuilder.h:152
MuonSeedBuilder::~MuonSeedBuilder
~MuonSeedBuilder()
Destructor.
Definition: MuonSeedBuilder.cc:108
MuonSeedCleaner::seedCleaner
std::vector< TrajectorySeed > seedCleaner(const edm::EventSetup &eventSetup, std::vector< TrajectorySeed > &seeds)
Cache pointer to geometry.
Definition: MuonSeedCleaner.cc:72
gather_cfg.cout
cout
Definition: gather_cfg.py:144
MuonDetLayerGeometry::backwardCSCLayers
const std::vector< const DetLayer * > & backwardCSCLayers() const
return the backward (-Z) CSC DetLayers, inside-out
Definition: MuonDetLayerGeometry.cc:156
MuonSeedBuilder::muonSeedCreate_
MuonSeedCreator * muonSeedCreate_
Create seed according to region (CSC, DT, Overlap)
Definition: MuonSeedBuilder.h:142
MuonSeedBuilder::muonLayers
const MuonDetLayerGeometry * muonLayers
Definition: MuonSeedBuilder.h:146
MuonSubdetId::DT
static constexpr int DT
Definition: MuonSubdetId.h:11
cuy.dh
dh
Definition: cuy.py:355
hgcalTopologyTester_cfi.layers
layers
Definition: hgcalTopologyTester_cfi.py:8
event
Definition: event.py:1
Generated for CMSSW Reference Manual by doxygen 1.8.11