CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrajectorySeedProducer.cc
Go to the documentation of this file.
1 #include <memory>
2 
7 
15 
18 
20 
23 
26 
30 
33 
36 
37 //Propagator withMaterial
39 //analyticalpropagator
40 //#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
41 
42 //
43 
44 //for debug only
45 //#define FAMOS_DEBUG
46 
48 {
49 
50  // The input tag for the beam spot
51  theBeamSpot = conf.getParameter<edm::InputTag>("beamSpot");
52 
53  // The name of the TrajectorySeed Collections
54  seedingAlgo = conf.getParameter<std::vector<std::string> >("seedingAlgo");
55  for ( unsigned i=0; i<seedingAlgo.size(); ++i )
56  produces<TrajectorySeedCollection>(seedingAlgo[i]);
57 
58  // The smallest true pT for a track candidate
59  pTMin = conf.getParameter<std::vector<double> >("pTMin");
60  if ( pTMin.size() != seedingAlgo.size() )
61  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
62  << " WARNING : pTMin does not have the proper size "
63  << std::endl;
64 
65  for ( unsigned i=0; i<pTMin.size(); ++i )
66  pTMin[i] *= pTMin[i]; // Cut is done of perp2() - CPU saver
67 
68  // The smallest number of Rec Hits for a track candidate
69  minRecHits = conf.getParameter<std::vector<unsigned int> >("minRecHits");
70  if ( minRecHits.size() != seedingAlgo.size() )
71  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
72  << " WARNING : minRecHits does not have the proper size "
73  << std::endl;
74  // Set the overall number hits to be checked
75  absMinRecHits = 0;
76  for ( unsigned ialgo=0; ialgo<minRecHits.size(); ++ialgo )
77  if ( minRecHits[ialgo] > absMinRecHits ) absMinRecHits = minRecHits[ialgo];
78 
79  // The smallest true impact parameters (d0 and z0) for a track candidate
80  maxD0 = conf.getParameter<std::vector<double> >("maxD0");
81  if ( maxD0.size() != seedingAlgo.size() )
82  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
83  << " WARNING : maxD0 does not have the proper size "
84  << std::endl;
85 
86  maxZ0 = conf.getParameter<std::vector<double> >("maxZ0");
87  if ( maxZ0.size() != seedingAlgo.size() )
88  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
89  << " WARNING : maxZ0 does not have the proper size "
90  << std::endl;
91 
92  // The name of the hit producer
93  hitProducer = conf.getParameter<edm::InputTag>("HitProducer");
94 
95  // The cuts for seed cleaning
96  seedCleaning = conf.getParameter<bool>("seedCleaning");
97 
98  // Number of hits needed for a seed
99  numberOfHits = conf.getParameter<std::vector<unsigned int> >("numberOfHits");
100  if ( numberOfHits.size() != seedingAlgo.size() )
101  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
102  << " WARNING : numberOfHits does not have the proper size "
103  << std::endl;
104 
105  // Seeding based on muons
106  selectMuons = conf.getParameter<bool>("selectMuons");
107 
108  // Layers
109  newSyntax = conf.getParameter<bool>("newSyntax");
110  if (newSyntax) {
111  std::vector<std::string> layerList = conf.getParameter<std::vector<std::string> >("layerList");
112  //for (unsigned i=0; i<layerList.size();i++) std::cout << "------- Layers = " << layerList[i] << std::endl;
113 
114  for(std::vector<std::string>::const_iterator it=layerList.begin(); it < layerList.end(); ++it) {
115  std::vector<LayerSpec> tempResult;
116  std::string line = *it;
118  while (pos != std::string::npos) {
119  pos=line.find("+");
120  std::string layer = line.substr(0, pos);
121  //
122  LayerSpec layerSpec;
123  layerSpec.name = line.substr(0, pos);;
124  // Possible names: BPix(1-3) || FPix(1-2)(pos, neg) || TIB(1-4) || TID(1-3)(pos, neg) || TOB(1-6) || TEC(1-9)(pos, neg)
125  // BPix(1-3)
126  if (layerSpec.name.substr(0,4)=="BPix" ) {
127  layerSpec.subDet=PXB;
128  layerSpec.side=BARREL;
129  layerSpec.idLayer = std::atoi(layerSpec.name.substr(4,1).c_str());
130  if (layerSpec.idLayer>3 || layerSpec.idLayer==0) {
131  throw cms::Exception("FastSimulation/Tracking/python")
132  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
133  << "Layers: " << layerSpec.name << ", number needs to be in range 1-3" << std::endl;
134  }
135  }
136  // FPix(1-2)(pos, neg)
137  else if (layerSpec.name.substr(0,4)=="FPix" ) {
138  layerSpec.subDet=PXD;
139  if(layerSpec.name.substr(layerSpec.name.size()-3)=="pos")
140  layerSpec.side = POS_ENDCAP;
141  else //no validation if it's not neg
142  layerSpec.side = NEG_ENDCAP;
143  layerSpec.idLayer = std::atoi(layerSpec.name.substr(4,1).c_str());
144  if (layerSpec.idLayer>2 || layerSpec.idLayer==0) {
145  throw cms::Exception("FastSimulation/Tracking/python")
146  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
147  << "Layers: " << layerSpec.name << ", number needs to be in range 1-2" << std::endl;
148  }
149  }
150  // TIB(1-4)
151  else if (layerSpec.name.substr(0,3)=="TIB" ) {
152  layerSpec.subDet=TIB;
153  layerSpec.side=BARREL;
154  layerSpec.idLayer = std::atoi(layerSpec.name.substr(3,1).c_str());
155  if (layerSpec.idLayer>4 || layerSpec.idLayer==0) {
156  throw cms::Exception("FastSimulation/Tracking/python")
157  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
158  << "Layers: " << layerSpec.name << ", number needs to be in range 1-4" << std::endl;
159  }
160  }
161  // TID(1-3)(pos, neg)
162  else if (layerSpec.name.substr(0,3)=="TID" ) {
163  layerSpec.subDet=TID;
164  if(layerSpec.name.substr(layerSpec.name.size()-3)=="pos")
165  layerSpec.side = POS_ENDCAP;
166  else
167  layerSpec.side = NEG_ENDCAP;
168  layerSpec.idLayer = std::atoi(layerSpec.name.substr(3,1).c_str());
169  if (layerSpec.idLayer>3 || layerSpec.idLayer==0) {
170  throw cms::Exception("FastSimulation/Tracking/python")
171  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
172  << "Layers: " << layerSpec.name << ", number needs to be in range 1-3" << std::endl;
173  }
174  }
175  // TOB(1-6)
176  else if (layerSpec.name.substr(0,3)=="TOB" ) {
177  layerSpec.subDet=TOB;
178  layerSpec.side=BARREL;
179  layerSpec.idLayer = std::atoi(layerSpec.name.substr(3,1).c_str());
180  if (layerSpec.idLayer>6 || layerSpec.idLayer==0) {
181  throw cms::Exception("FastSimulation/Tracking/python")
182  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
183  << "Layers: " << layerSpec.name << ", number needs to be in range 1-6" << std::endl;
184  }
185  }
186  // TEC(1-9)(pos, neg)
187  else if (layerSpec.name.substr(0,3)=="TEC" ) {
188  layerSpec.subDet=TEC;
189  if(layerSpec.name.substr(layerSpec.name.size()-3)=="pos")
190  layerSpec.side = POS_ENDCAP;
191  else
192  layerSpec.side = NEG_ENDCAP;
193  layerSpec.idLayer = std::atoi(layerSpec.name.substr(3,1).c_str());
194  if (layerSpec.idLayer>9 || layerSpec.idLayer==0) {
195  throw cms::Exception("FastSimulation/Tracking/python")
196  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
197  << "Layers: " << layerSpec.name << ", number needs to be in range 1-9" << std::endl;
198  }
199  }
200  else {
201  throw cms::Exception("FastSimulation/Tracking/python")
202  << "Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
203  << "Layer: " << layerSpec.name << ", shouldn't exist" << std::endl
204  << "Case sensitive names: BPix FPix TIB TID TOB TEC" << std::endl;
205  }
206  //
207  tempResult.push_back(layerSpec);
208  line=line.substr(pos+1,std::string::npos);
209  }
210  theLayersInSets.push_back(tempResult);
211  }
212 
213  //prints theLayersInSets
214  /* for (std::vector<std::vector<LayerSpec> >::const_iterator it = theLayersInSets.begin(); it != theLayersInSets.end(); ++it) {
215  for (std::vector<LayerSpec>::const_iterator is = it->begin(); is != it->end(); ++is) {
216  std::cout << is->name << " | " << is->subDet << " | " << is->idLayer << " | " << is->side << std::endl;
217  }
218  std::cout << "---" << std::endl;
219  } */
220  } else {
221  // TO BE DELETED (AG)
223  conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectorNumber");
224  if ( firstHitSubDetectorNumber.size() != seedingAlgo.size() )
225  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
226  << " WARNING : firstHitSubDetectorNumber does not have the proper size "
227  << std::endl;
228 
229  std::vector<unsigned int> firstSubDets =
230  conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectors");
231  unsigned isub1 = 0;
232  unsigned check1 = 0;
233  firstHitSubDetectors.resize(seedingAlgo.size());
234  for ( unsigned ialgo=0; ialgo<firstHitSubDetectorNumber.size(); ++ialgo ) {
235  check1 += firstHitSubDetectorNumber[ialgo];
236  for ( unsigned idet=0; idet<firstHitSubDetectorNumber[ialgo]; ++idet ) {
237  firstHitSubDetectors[ialgo].push_back(firstSubDets[isub1++]);
238  }
239  }
240  if ( firstSubDets.size() != check1 )
241  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
242  << " WARNING : firstHitSubDetectors does not have the proper size (should be " << check1 << ")"
243  << std::endl;
244 
245 
247  conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectorNumber");
248  if ( secondHitSubDetectorNumber.size() != seedingAlgo.size() )
249  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
250  << " WARNING : secondHitSubDetectorNumber does not have the proper size "
251  << std::endl;
252 
253  std::vector<unsigned int> secondSubDets =
254  conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectors");
255  unsigned isub2 = 0;
256  unsigned check2 = 0;
257  secondHitSubDetectors.resize(seedingAlgo.size());
258  for ( unsigned ialgo=0; ialgo<secondHitSubDetectorNumber.size(); ++ialgo ) {
259  check2 += secondHitSubDetectorNumber[ialgo];
260  for ( unsigned idet=0; idet<secondHitSubDetectorNumber[ialgo]; ++idet ) {
261  secondHitSubDetectors[ialgo].push_back(secondSubDets[isub2++]);
262  }
263  }
264  if ( secondSubDets.size() != check2 )
265  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
266  << " WARNING : secondHitSubDetectors does not have the proper size (should be " << check2 << ")"
267  << std::endl;
268 
270  conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectorNumber");
271  if ( thirdHitSubDetectorNumber.size() != seedingAlgo.size() )
272  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
273  << " WARNING : thirdHitSubDetectorNumber does not have the proper size "
274  << std::endl;
275 
276  std::vector<unsigned int> thirdSubDets =
277  conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectors");
278  unsigned isub3 = 0;
279  unsigned check3 = 0;
280  thirdHitSubDetectors.resize(seedingAlgo.size());
281  for ( unsigned ialgo=0; ialgo<thirdHitSubDetectorNumber.size(); ++ialgo ) {
282  check3 += thirdHitSubDetectorNumber[ialgo];
283  for ( unsigned idet=0; idet<thirdHitSubDetectorNumber[ialgo]; ++idet ) {
284  thirdHitSubDetectors[ialgo].push_back(thirdSubDets[isub3++]);
285  }
286  }
287  if ( thirdSubDets.size() != check3 )
288  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
289  << " WARNING : thirdHitSubDetectors does not have the proper size (should be " << check3 << ")"
290  << std::endl;
291  }
292 
293  originRadius = conf.getParameter<std::vector<double> >("originRadius");
294  if ( originRadius.size() != seedingAlgo.size() )
295  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
296  << " WARNING : originRadius does not have the proper size "
297  << std::endl;
298 
299  originHalfLength = conf.getParameter<std::vector<double> >("originHalfLength");
300  if ( originHalfLength.size() != seedingAlgo.size() )
301  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
302  << " WARNING : originHalfLength does not have the proper size "
303  << std::endl;
304 
305  originpTMin = conf.getParameter<std::vector<double> >("originpTMin");
306  if ( originpTMin.size() != seedingAlgo.size() )
307  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
308  << " WARNING : originpTMin does not have the proper size "
309  << std::endl;
310 
311  primaryVertices = conf.getParameter<std::vector<edm::InputTag> >("primaryVertices");
312  if ( primaryVertices.size() != seedingAlgo.size() )
313  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
314  << " WARNING : primaryVertices does not have the proper size "
315  << std::endl;
316 
317  zVertexConstraint = conf.getParameter<std::vector<double> >("zVertexConstraint");
318  if ( zVertexConstraint.size() != seedingAlgo.size() )
319  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
320  << " WARNING : zVertexConstraint does not have the proper size "
321  << std::endl;
322  //removed - }
323 }
324 
325 // Virtual destructor needed.
327 
328  if(thePropagator) delete thePropagator;
329 
330  // do nothing
331 #ifdef FAMOS_DEBUG
332  std::cout << "TrajectorySeedProducer destructed" << std::endl;
333 #endif
334 
335 }
336 
337 void
339 
340  //services
341  // es.get<TrackerRecoGeometryRecord>().get(theGeomSearchTracker);
342 
346 
347 
348  es.get<IdealMagneticFieldRecord>().get(magField);
349  es.get<TrackerDigiGeometryRecord>().get(geometry);
350  es.get<MagneticFieldMapRecord>().get(magFieldMap);
351 
352  theMagField = &(*magField);
353  theGeometry = &(*geometry);
354  theFieldMap = &(*magFieldMap);
355 
357 
358  const GlobalPoint g(0.,0.,0.);
359 
360 }
361 
362  // Functions that gets called by framework every event
363 void
365 
366 
367  // if( seedingAlgo[0] == "FourthPixelLessPairs") std::cout << "Seed producer in 4th iteration " << std::endl;
368 
369 #ifdef FAMOS_DEBUG
370  std::cout << "################################################################" << std::endl;
371  std::cout << " TrajectorySeedProducer produce init " << std::endl;
372 #endif
373 
374  //Retrieve tracker topology from geometry
376  es.get<IdealGeometryRecord>().get(tTopoHand);
377  const TrackerTopology *tTopo=tTopoHand.product();
378 
379 
380  unsigned nSimTracks = 0;
381  unsigned nTracksWithHits = 0;
382  unsigned nTracksWithPT = 0;
383  unsigned nTracksWithD0Z0 = 0;
384  // unsigned nTrackCandidates = 0;
385  PTrajectoryStateOnDet initialState;
386 
387  // Output
388  std::vector<TrajectorySeedCollection*>
389  output(seedingAlgo.size(),static_cast<TrajectorySeedCollection*>(0));
390  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
391  // std::auto_ptr<TrajectorySeedCollection> p(new TrajectorySeedCollection );
392  output[ialgo] = new TrajectorySeedCollection;
393  }
394 
395  // Beam spot
396  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
397  e.getByLabel(theBeamSpot,recoBeamSpotHandle);
398  math::XYZPoint BSPosition_ = recoBeamSpotHandle->position();
399 
400  //not used anymore. take the value from the py
401 
402  //double sigmaZ=recoBeamSpotHandle->sigmaZ();
403  //double sigmaZ0Error=recoBeamSpotHandle->sigmaZ0Error();
404  //double sigmaz0=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
405  x0 = BSPosition_.X();
406  y0 = BSPosition_.Y();
407  z0 = BSPosition_.Z();
408 
409  // SimTracks and SimVertices
411  e.getByLabel("famosSimHits",theSimTracks);
412 
414  e.getByLabel("famosSimHits",theSimVtx);
415 
416 #ifdef FAMOS_DEBUG
417  std::cout << " Step A: SimTracks found " << theSimTracks->size() << std::endl;
418 #endif
419 
420  // edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
422  e.getByLabel(hitProducer, theGSRecHits);
423 
424  // No tracking attempted if no hits (but put an empty collection in the event)!
425 #ifdef FAMOS_DEBUG
426  std::cout << " Step B: Full GS RecHits found " << theGSRecHits->size() << std::endl;
427 #endif
428  if(theGSRecHits->size() == 0) {
429  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
430  std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
431  e.put(p,seedingAlgo[ialgo]);
432  }
433  return;
434  }
435 
436  // Primary vertices
437  vertices = std::vector<const reco::VertexCollection*>
438  (seedingAlgo.size(),static_cast<const reco::VertexCollection*>(0));
439  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
440  //PAT Attempt!!!!
441  //originHalfLength[ialgo] = 3.*sigmaz0; // Overrides the configuration
443  bool isVertexCollection = e.getByLabel(primaryVertices[ialgo],aHandle);
444  if (!isVertexCollection ) continue;
445  vertices[ialgo] = &(*aHandle);
446  }
447 
448 #ifdef FAMOS_DEBUG
449  std::cout << " Step C: Loop over the RecHits, track by track " << std::endl;
450 #endif
451 
452  // The vector of simTrack Id's carrying GSRecHits
453  const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
454 
455  // loop over SimTrack Id's
456  for ( unsigned tkId=0; tkId != theSimTrackIds.size(); ++tkId ) {
457 
458 #ifdef FAMOS_DEBUG
459  std::cout << "Track number " << tkId << std::endl;
460 #endif
461 
462  ++nSimTracks;
463  unsigned simTrackId = theSimTrackIds[tkId];
464  const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
465 #ifdef FAMOS_DEBUG
466  std::cout << "Pt = " << std::sqrt(theSimTrack.momentum().Perp2())
467  << " eta " << theSimTrack.momentum().Eta()
468  << " pdg ID " << theSimTrack.type()
469  << std::endl;
470 #endif
471 
472  // Select only muons, if requested
473  if (selectMuons && abs(theSimTrack.type()) != 13) continue;
474 
475  // Check that the sim track comes from the main vertex (loose cut)
476  int vertexIndex = theSimTrack.vertIndex();
477  const SimVertex& theSimVertex = (*theSimVtx)[vertexIndex];
478 #ifdef FAMOS_DEBUG
479  std::cout << " o SimTrack " << theSimTrack << std::endl;
480  std::cout << " o SimVertex " << theSimVertex << std::endl;
481 #endif
482 
483  BaseParticlePropagator theParticle =
485  RawParticle(XYZTLorentzVector(theSimTrack.momentum().px(),
486  theSimTrack.momentum().py(),
487  theSimTrack.momentum().pz(),
488  theSimTrack.momentum().e()),
489  XYZTLorentzVector(theSimVertex.position().x(),
490  theSimVertex.position().y(),
491  theSimVertex.position().z(),
492  theSimVertex.position().t())),
493  0.,0.,4.);
494  theParticle.setCharge((*theSimTracks)[simTrackId].charge());
495 
496  SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
497  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
498  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd = theRecHitRange.second;
503 
504  // Check the number of layers crossed
505  unsigned numberOfRecHits = 0;
506  TrackerRecHit previousHit, currentHit;
507  for ( iterRecHit = theRecHitRangeIteratorBegin;
508  iterRecHit != theRecHitRangeIteratorEnd;
509  ++iterRecHit) {
510  previousHit = currentHit;
511  currentHit = TrackerRecHit(&(*iterRecHit),theGeometry,tTopo);
512  if ( currentHit.isOnTheSameLayer(previousHit) ) continue;
513  ++numberOfRecHits;
514  if ( numberOfRecHits == absMinRecHits ) break;
515  }
516 
517  // Loop on the successive seedings
518  for ( unsigned int ialgo = 0; ialgo < seedingAlgo.size(); ++ialgo ) {
519 #ifdef FAMOS_DEBUG
520  std::cout << "Algo " << seedingAlgo[ialgo] << std::endl;
521 #endif
522 
523  // Request a minimum number of RecHits for the track to give a seed.
524 #ifdef FAMOS_DEBUG
525  std::cout << "The number of RecHits = " << numberOfRecHits << std::endl;
526 #endif
527  if ( numberOfRecHits < minRecHits[ialgo] ) continue;
528  ++nTracksWithHits;
529 
530  // Request a minimum pT for the sim track
531  if ( theSimTrack.momentum().Perp2() < pTMin[ialgo] ) continue;
532  ++nTracksWithPT;
533 
534  // Cut on sim track impact parameters
535  if ( theParticle.xyImpactParameter(x0,y0) > maxD0[ialgo] ) continue;
536  if ( fabs( theParticle.zImpactParameter(x0,y0) - z0 ) > maxZ0[ialgo] ) continue;
537  ++nTracksWithD0Z0;
538 
539  std::vector<TrackerRecHit >
540  theSeedHits(numberOfHits[ialgo],
541  static_cast<TrackerRecHit >(TrackerRecHit()));
542  TrackerRecHit& theSeedHits0 = theSeedHits[0];
543  TrackerRecHit& theSeedHits1 = theSeedHits[1];
544  TrackerRecHit& theSeedHits2 = theSeedHits[2];
545 
546  bool compatible = false;
547 
548  for ( iterRecHit1 = theRecHitRangeIteratorBegin; iterRecHit1 != theRecHitRangeIteratorEnd; ++iterRecHit1) {
549  //std::cout << (*iterRecHit1).localPosition().phi() << " | J - iterRecHit1"<< std::endl; //
550  theSeedHits[0] = TrackerRecHit(&(*iterRecHit1),theGeometry,tTopo);
551 #ifdef FAMOS_DEBUG
552  std::cout << "The first hit position = " << theSeedHits0.globalPosition() << std::endl;
553  std::cout << "The first hit subDetId = " << theSeedHits0.subDetId() << std::endl;
554  std::cout << "The first hit layer = " << theSeedHits0.layerNumber() << std::endl;
555 #endif
556  // Check if inside the requested detectors
557  bool isInside = true;
558  if (!selectMuons) {
559  //(newSyntax) ? std::cout << "TRUE " : std::cout << "FALSE "; //J-
560  if (newSyntax)
561  isInside = false; // AG placeholder true
562  else
563  isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
564  // bool isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
565  if ( isInside ) continue;
566  }
567 
568  // Check if on requested detectors
569  // bool isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo]);
570  bool isOndet = true;
571  if (!selectMuons) {
572  if (newSyntax)
573  isOndet = theSeedHits[0].isOnRequestedDet(theLayersInSets);
574  else
575  isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
576  //std::cout << firstHitSubDetectors[ialgo][0] << " | " << seedingAlgo[ialgo] << " " << std::endl; //seedingAlgo[iAlgo]: PixelTriplet, LowPtPixelTriplets, PixelPair, DetachedPixelTriplets, MixedTriplets, PixelLessPairs, TobTecLayerPairs......
577  // bool isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
578  // if ( !isOndet ) break;
579  if ( !isOndet ) continue;
580  }
581 
582 #ifdef FAMOS_DEBUG
583  std::cout << "Apparently the first hit is on the requested detector! " << std::endl;
584 #endif
585  for ( iterRecHit2 = iterRecHit1+1; iterRecHit2 != theRecHitRangeIteratorEnd; ++iterRecHit2) {
586  theSeedHits[1] = TrackerRecHit(&(*iterRecHit2),theGeometry,tTopo);
587 #ifdef FAMOS_DEBUG
588  std::cout << "The second hit position = " << theSeedHits1.globalPosition() << std::endl;
589  std::cout << "The second hit subDetId = " << theSeedHits1.subDetId() << std::endl;
590  std::cout << "The second hit layer = " << theSeedHits1.layerNumber() << std::endl;
591 #endif
592 
593  if (!selectMuons) {
594  // Check if inside the requested detectors
595  if (newSyntax)
596  isInside = false; // AG placeholder true
597  else
598  isInside = theSeedHits1.subDetId() < secondHitSubDetectors[ialgo][0];
599  if ( isInside ) continue;
600 
601  // Check if on requested detectors
602  if (newSyntax)
603  isOndet = theSeedHits[0].isOnRequestedDet(theLayersInSets, theSeedHits[1]);
604  else
605  isOndet = theSeedHits1.isOnRequestedDet(secondHitSubDetectors[ialgo], seedingAlgo[ialgo]);
606  if ( !isOndet ) break;
607  }
608  // Check if on the same layer as previous hit
609  if ( theSeedHits1.isOnTheSameLayer(theSeedHits0) ) continue;
610 
611 #ifdef FAMOS_DEBUG
612  std::cout << "Apparently the second hit is on the requested detector! " << std::endl;
613 #endif
614  GlobalPoint gpos1 = theSeedHits0.globalPosition();
615  GlobalPoint gpos2 = theSeedHits1.globalPosition();
616  bool forward = theSeedHits0.isForward();
617  double error = std::sqrt(theSeedHits0.largerError()+theSeedHits1.largerError());
618  // compatible = compatibleWithVertex(gpos1,gpos2,ialgo);
619  //added out of desperation
620  if(seedingAlgo[0] == "PixelLess" || seedingAlgo[0] == "TobTecLayerPairs"){
621  compatible = true;
622  //std::cout << "Algo " << seedingAlgo[0] << "Det/layer = " << theSeedHits0.subDetId() << "/" << theSeedHits0.layerNumber() << std::endl;
623  } else {
624  compatible = compatibleWithBeamAxis(gpos1,gpos2,error,forward,ialgo);
625  }
626 
627 #ifdef FAMOS_DEBUG
628  std::cout << "Algo" << seedingAlgo[0] << "\t Are the two hits compatible with the PV? " << compatible << std::endl;
629 #endif
630 
631  if (!selectMuons) {
632  // Check if the pair is on the requested dets
633  if ( numberOfHits[ialgo] == 2 ) {
634 
635  if ( seedingAlgo[0] == "ThirdMixedPairs" ){
636  compatible = compatible && theSeedHits[0].makesAPairWith3rd(theSeedHits[1]);
637  } else {
638  compatible = compatible && theSeedHits[0].makesAPairWith(theSeedHits[1]);
639  //check
640  /*
641  if((seedingAlgo[0] == "PixelLess" || seedingAlgo[0] == "TobTecLayerPairs") && !compatible)
642  std::cout << "NOT Compatible " << seedingAlgo[0]
643  << "Hit 1 Det/layer/ring = " << theSeedHits0.subDetId() << "/" << theSeedHits0.layerNumber() << "/" << theSeedHits0.ringNumber()
644  << "\tHit 2 Det/layer/ring = " << theSeedHits1.subDetId() << "/" << theSeedHits1.layerNumber() << "/" << theSeedHits1.ringNumber() << std::endl;
645  */
646  }
647  }
648  }
649 
650  // Reject non suited pairs
651  if ( !compatible ) continue;
652 
653 #ifdef FAMOS_DEBUG
654  std::cout << "Pair kept! " << std::endl;
655 #endif
656 
657  // Leave here if only two hits are required.
658  if ( numberOfHits[ialgo] == 2 ) break;
659 
660  compatible = false;
661  // Check if there is a third satisfying hit otherwise
662  for ( iterRecHit3 = iterRecHit2+1; iterRecHit3 != theRecHitRangeIteratorEnd; ++iterRecHit3) {
663 
664  theSeedHits[2] = TrackerRecHit(&(*iterRecHit3),theGeometry,tTopo);
665 #ifdef FAMOS_DEBUG
666  std::cout << "The third hit position = " << theSeedHits2.globalPosition() << std::endl;
667  std::cout << "The third hit subDetId = " << theSeedHits2.subDetId() << std::endl;
668  std::cout << "The third hit layer = " << theSeedHits2.layerNumber() << std::endl;
669 #endif
670 
671  if (!selectMuons) {
672  // Check if inside the requested detectors
673  if (newSyntax)
674  isInside = false; // AG placeholder
675  else
676  isInside = theSeedHits2.subDetId() < thirdHitSubDetectors[ialgo][0];
677  if ( isInside ) continue;
678 
679  // Check if on requested detectors
680  if (newSyntax)
681  isOndet = theSeedHits[0].isOnRequestedDet(theLayersInSets, theSeedHits[1], theSeedHits[2]);
682  else
683  isOndet = theSeedHits2.isOnRequestedDet(thirdHitSubDetectors[ialgo], seedingAlgo[ialgo]);
684  // if ( !isOndet ) break;
685  if ( !isOndet ) continue;
686  }
687 
688  // Check if on the same layer as previous hit
689  compatible = !(theSeedHits2.isOnTheSameLayer(theSeedHits1));
690 
691  // Check if the triplet is on the requested det combination
692  if (!selectMuons) compatible = compatible && theSeedHits[0].makesATripletWith(theSeedHits[1],theSeedHits[2]); //J- maybe it's not necessary, as newSyntax layerlist is already checking?
693 
694 #ifdef FAMOS_DEBUG
695  if ( compatible )
696  std::cout << "Apparently the third hit is on the requested detector! " << std::endl;
697 #endif
698 
699  if ( compatible ) break;
700 
701  }
702 
703  if ( compatible ) break;
704 
705  }
706 
707  if ( compatible ) break;
708 
709  }
710 
711  // There is no compatible seed for this track with this seeding algorithm
712  // Go to next algo
713  if ( !compatible ) continue;
714 #ifdef FAMOS_DEBUG
715  if ( compatible )
716  std::cout << "@@@ There is at least a compatible seed" << std::endl;
717  else
718  std::cout << "@@@ There is no compatible seed" << std::endl;
719 #endif
720 
721 
722 #ifdef FAMOS_DEBUG
723  std::cout << "Preparing to create the TrajectorySeed" << std::endl;
724 #endif
725  // The seed is validated -> include in the collection
726  // 1) Create the vector of RecHits
728  for ( unsigned ih=0; ih<theSeedHits.size(); ++ih ) {
729  TrackingRecHit* aTrackingRecHit = theSeedHits[ih].hit()->clone();
730  recHits.push_back(aTrackingRecHit);
731  }
732 #ifdef FAMOS_DEBUG
733  std::cout << "with " << recHits.size() << " hits." << std::endl;
734 #endif
735 
736  // 2) Create the initial state
737  // a) origin vertex
738  GlobalPoint position((*theSimVtx)[vertexIndex].position().x(),
739  (*theSimVtx)[vertexIndex].position().y(),
740  (*theSimVtx)[vertexIndex].position().z());
741 
742  // b) initial momentum
743  GlobalVector momentum( (*theSimTracks)[simTrackId].momentum().x() ,
744  (*theSimTracks)[simTrackId].momentum().y() ,
745  (*theSimTracks)[simTrackId].momentum().z() );
746  // c) electric charge
747  float charge = (*theSimTracks)[simTrackId].charge();
748  // -> inital parameters
749  GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
750  // -> large initial errors
751  AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();
752  // errorMatrix = errorMatrix * 10;
753 
754  //this line help the fit succeed in the case of pixelless tracks (4th and 5th iteration)
755  //for the future: probably the best thing is to use the mini-kalmanFilter
756  if(theSeedHits0.subDetId() !=1 || theSeedHits0.subDetId() !=2) errorMatrix = errorMatrix * 0.0000001;
757 
758 
759 
760 #ifdef FAMOS_DEBUG
761  std::cout << "TrajectorySeedProducer: SimTrack parameters " << std::endl;
762  std::cout << "\t\t pT = " << (*theSimTracks)[simTrackId].momentum().Pt() << std::endl;
763  std::cout << "\t\t eta = " << (*theSimTracks)[simTrackId].momentum().Eta() << std::endl;
764  std::cout << "\t\t phi = " << (*theSimTracks)[simTrackId].momentum().Phi() << std::endl;
765  std::cout << "TrajectorySeedProducer: AlgebraicSymMatrix " << errorMatrix << std::endl;
766 #endif
767  CurvilinearTrajectoryError initialError(errorMatrix);
768  // -> initial state
769  FreeTrajectoryState initialFTS(initialParams, initialError);
770 #ifdef FAMOS_DEBUG
771  std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
772 #endif
773  // const GeomDetUnit* initialLayer = theGeometry->idToDetUnit( recHits.front().geographicalId() );
774  const GeomDet* initialLayer = theGeometry->idToDet( recHits.front().geographicalId() );
775 
776  //this is wrong because the FTS is defined at vertex, and it need to be properly propagated.
777  // const TrajectoryStateOnSurface initialTSOS(initialFTS, initialLayer->surface());
778 
779  const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
780  if (!initialTSOS.isValid()) continue;
781 
782 #ifdef FAMOS_DEBUG
783  std::cout << "TrajectorySeedProducer: TSOS global momentum " << initialTSOS.globalMomentum() << std::endl;
784  std::cout << "\t\t\tpT = " << initialTSOS.globalMomentum().perp() << std::endl;
785  std::cout << "\t\t\teta = " << initialTSOS.globalMomentum().eta() << std::endl;
786  std::cout << "\t\t\tphi = " << initialTSOS.globalMomentum().phi() << std::endl;
787  std::cout << "TrajectorySeedProducer: TSOS local momentum " << initialTSOS.localMomentum() << std::endl;
788  std::cout << "TrajectorySeedProducer: TSOS local error " << initialTSOS.localError().positionError() << std::endl;
789  std::cout << "TrajectorySeedProducer: TSOS local error matrix " << initialTSOS.localError().matrix() << std::endl;
790  std::cout << "TrajectorySeedProducer: TSOS surface side " << initialTSOS.surfaceSide() << std::endl;
791 #endif
792  stateOnDet(initialTSOS,
793  recHits.front().geographicalId().rawId(),
794  initialState);
795  // Create a new Trajectory Seed
796  output[ialgo]->push_back(TrajectorySeed(initialState, recHits, alongMomentum));
797 #ifdef FAMOS_DEBUG
798  std::cout << "Trajectory seed created ! " << std::endl;
799 #endif
800  break;
801  // End of the loop over seeding algorithms
802  }
803  // End on the loop over simtracks
804  }
805 
806  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
807  std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
808  e.put(p,seedingAlgo[ialgo]);
809  }
810 
811 }
812 
813 
814 // This is a copy of a method in
815 // TrackingTools/TrajectoryState/src/TrajectoryStateTransform.cc
816 // but it does not return a pointer (thus avoiding a memory leak)
817 // In addition, it's also CPU more efficient, because
818 // ts.localError().matrix() is not copied
819 void
821  unsigned int detid,
822  PTrajectoryStateOnDet& pts) const
823 {
824 
825  const AlgebraicSymMatrix55& m = ts.localError().matrix();
826 
827  int dim = 5;
828 
829  float localErrors[15];
830  int k = 0;
831  for (int i=0; i<dim; ++i) {
832  for (int j=0; j<=i; ++j) {
833  localErrors[k++] = m(i,j);
834  }
835  }
836  int surfaceSide = static_cast<int>(ts.surfaceSide());
837 
839  localErrors, detid,
840  surfaceSide);
841 }
842 
843 bool
845  GlobalPoint& gpos2,
846  double error,
847  bool forward,
848  unsigned algo) const {
849 
850  if ( !seedCleaning ) return true;
851 
852  // The hits 1 and 2 positions, in HepLorentzVector's
853  XYZTLorentzVector thePos1(gpos1.x(),gpos1.y(),gpos1.z(),0.);
854  XYZTLorentzVector thePos2(gpos2.x(),gpos2.y(),gpos2.z(),0.);
855 #ifdef FAMOS_DEBUG
856  std::cout << "ThePos1 = " << thePos1 << std::endl;
857  std::cout << "ThePos2 = " << thePos2 << std::endl;
858 #endif
859 
860 
861  // Create new particles that pass through the second hit with pT = ptMin
862  // and charge = +/-1
863 
864  // The momentum direction is by default joining the two hits
865  XYZTLorentzVector theMom2 = (thePos2-thePos1);
866 
867  // The corresponding RawParticle, with an (irrelevant) electric charge
868  // (The charge is determined in the next step)
869  ParticlePropagator myPart(theMom2,thePos2,1.,theFieldMap);
870 
874  bool intersect = myPart.propagateToBeamCylinder(thePos1,originRadius[algo]*1.);
875  if ( !intersect ) return false;
876 
877 #ifdef FAMOS_DEBUG
878  std::cout << "MyPart R = " << myPart.R() << "\t Z = " << myPart.Z()
879  << "\t pT = " << myPart.Pt() << std::endl;
880 #endif
881 
882  // Check if the constraints are satisfied
883  // 1. pT at cylinder with radius originRadius
884  if ( myPart.Pt() < originpTMin[algo] ) return false;
885 
886  // 2. Z compatible with beam spot size
887  if ( fabs(myPart.Z()-z0) > originHalfLength[algo] ) return false;
888 
889  // 3. Z compatible with one of the primary vertices (always the case if no primary vertex)
890  const reco::VertexCollection* theVertices = vertices[algo];
891  if (!theVertices) return true;
892  unsigned nVertices = theVertices->size();
893  if ( !nVertices || zVertexConstraint[algo] < 0. ) return true;
894  // Radii of the two hits with respect to the beam spot position
895  double R1 = std::sqrt ( (thePos1.X()-x0)*(thePos1.X()-x0)
896  + (thePos1.Y()-y0)*(thePos1.Y()-y0) );
897  double R2 = std::sqrt ( (thePos2.X()-x0)*(thePos2.X()-x0)
898  + (thePos2.Y()-y0)*(thePos2.Y()-y0) );
899  // Loop on primary vertices
900  for ( unsigned iv=0; iv<nVertices; ++iv ) {
901  // Z position of the primary vertex
902  double zV = (*theVertices)[iv].z();
903  // Constraints on the inner hit
904  double checkRZ1 = forward ?
905  (thePos1.Z()-zV+zVertexConstraint[algo]) / (thePos2.Z()-zV+zVertexConstraint[algo]) * R2 :
906  -zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV+zVertexConstraint[algo]);
907  double checkRZ2 = forward ?
908  (thePos1.Z()-zV-zVertexConstraint[algo])/(thePos2.Z()-zV-zVertexConstraint[algo]) * R2 :
909  +zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV-zVertexConstraint[algo]);
910  double checkRZmin = std::min(checkRZ1,checkRZ2)-3.*error;
911  double checkRZmax = std::max(checkRZ1,checkRZ2)+3.*error;
912  // Check if the innerhit is within bounds
913  bool compat = forward ?
914  checkRZmin < R1 && R1 < checkRZmax :
915  checkRZmin < thePos1.Z()-zV && thePos1.Z()-zV < checkRZmax;
916  // If it is, just return ok
917  if ( compat ) return compat;
918  }
919  // Otherwise, return not ok
920  return false;
921 
922 }
923 
924 
bool compatibleWithBeamAxis(GlobalPoint &gpos1, GlobalPoint &gpos2, double error, bool forward, unsigned algo) const
void setCharge(float q)
set the MEASURED charge
Definition: RawParticle.cc:138
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
std::vector< unsigned > minRecHits
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
double zImpactParameter(double x0=0, double y0=0.) const
Longitudinal impact parameter.
std::vector< unsigned int > firstHitSubDetectorNumber
T perp() const
Definition: PV3DBase.h:72
virtual void beginRun(edm::Run const &run, const edm::EventSetup &es) override
std::vector< std::vector< LayerSpec > > theLayersInSets
TrajectorySeedProducer(const edm::ParameterSet &conf)
const LocalTrajectoryParameters & localParameters() const
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
std::vector< double > pTMin
size_type size() const
Definition: OwnVector.h:247
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
bool propagateToBeamCylinder(const XYZTLorentzVector &v, double radius=0.)
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
std::vector< double > originRadius
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const MagneticFieldMap * theFieldMap
void stateOnDet(const TrajectoryStateOnSurface &ts, unsigned int detid, PTrajectoryStateOnDet &pts) const
A mere copy (without memory leak) of an existing tracking method.
std::vector< double > originpTMin
uint16_t size_type
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
double charge(const std::vector< uint8_t > &Ampls)
GlobalPoint globalPosition() const
The global position.
Definition: TrackerRecHit.h:99
double largerError()
LocalError positionError() const
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
float float float z
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
bool isOnTheSameLayer(const TrackerRecHit &other) const
Check if two hits are on the same layer of the same subdetector.
std::vector< unsigned int > numberOfHits
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
LocalVector localMomentum() const
PropagatorWithMaterial * thePropagator
double R() const
vertex radius
Definition: RawParticle.h:277
static const double pts[33]
Definition: Constants.h:30
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:43
void push_back(D *&d)
Definition: OwnVector.h:273
std::vector< std::vector< unsigned int > > secondHitSubDetectors
virtual void produce(edm::Event &e, const edm::EventSetup &es) override
std::vector< TrajectorySeed > TrajectorySeedCollection
const T & max(const T &a, const T &b)
std::vector< std::vector< unsigned int > > thirdHitSubDetectors
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< double > zVertexConstraint
SurfaceSide surfaceSide() const
Position relative to material, defined relative to momentum vector.
T z() const
Definition: PV3DBase.h:64
double Z() const
z of vertex
Definition: RawParticle.h:275
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
const AlgebraicSymMatrix55 & matrix() const
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:26
const LocalTrajectoryError & localError() const
const TrackerGeometry * theGeometry
virtual const GeomDet * idToDet(DetId) const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
GlobalVector momentum() const
tuple conf
Definition: dbtoconf.py:185
std::vector< unsigned int > thirdHitSubDetectorNumber
int k[5][pyjets_maxn]
virtual TrackingRecHit * clone() const =0
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:29
bool isForward() const
Is it a forward hit ?
Definition: TrackerRecHit.h:93
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
const T & get() const
Definition: EventSetup.h:55
std::vector< unsigned int > secondHitSubDetectorNumber
T const * product() const
Definition: ESHandle.h:62
unsigned int layerNumber() const
The Layer Number.
Definition: TrackerRecHit.h:84
T eta() const
Definition: PV3DBase.h:76
ESHandle< TrackerGeometry > geometry
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:40
bool isOnRequestedDet(const std::vector< unsigned int > &whichDet, const std::string &seedingAlgo) const
Check if the hit is on one of the requested detector.
GlobalVector globalMomentum() const
std::vector< const reco::VertexCollection * > vertices
std::vector< double > originHalfLength
const math::XYZTLorentzVectorD & momentum() const
particle info...
Definition: CoreSimTrack.h:36
std::vector< double > maxZ0
unsigned int subDetId() const
The subdet Id.
Definition: TrackerRecHit.h:81
std::vector< std::string > seedingAlgo
tuple cout
Definition: gather_cfg.py:121
std::vector< edm::InputTag > primaryVertices
std::vector< std::vector< unsigned int > > firstHitSubDetectors
std::vector< double > maxD0
DetId geographicalId() const
Definition: DDAxes.h:10
T x() const
Definition: PV3DBase.h:62
reference front()
Definition: OwnVector.h:348
Definition: Run.h:41
double xyImpactParameter(double x0=0., double y0=0.) const
Transverse impact parameter.
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:15
const MagneticField * theMagField