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