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 
45 //for debug only
46 //#define FAMOS_DEBUG
47 
49 {
50 
51  // The input tag for the beam spot
52  theBeamSpot = conf.getParameter<edm::InputTag>("beamSpot");
53 
54  // The name of the TrajectorySeed Collections
55  seedingAlgo = conf.getParameter<std::vector<std::string> >("seedingAlgo");
56  for ( unsigned i=0; i<seedingAlgo.size(); ++i )
57  produces<TrajectorySeedCollection>(seedingAlgo[i]);
58 
59  // The smallest true pT for a track candidate
60  pTMin = conf.getParameter<std::vector<double> >("pTMin");
61  if ( pTMin.size() != seedingAlgo.size() )
62  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
63  << " WARNING : pTMin does not have the proper size "
64  << std::endl;
65 
66  for ( unsigned i=0; i<pTMin.size(); ++i )
67  pTMin[i] *= pTMin[i]; // Cut is done of perp2() - CPU saver
68 
69  // The smallest number of Rec Hits for a track candidate
70  minRecHits = conf.getParameter<std::vector<unsigned int> >("minRecHits");
71  if ( minRecHits.size() != seedingAlgo.size() )
72  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
73  << " WARNING : minRecHits does not have the proper size "
74  << std::endl;
75  // Set the overall number hits to be checked
76  absMinRecHits = 0;
77  for ( unsigned ialgo=0; ialgo<minRecHits.size(); ++ialgo )
78  if ( minRecHits[ialgo] > absMinRecHits ) absMinRecHits = minRecHits[ialgo];
79 
80  // The smallest true impact parameters (d0 and z0) for a track candidate
81  maxD0 = conf.getParameter<std::vector<double> >("maxD0");
82  if ( maxD0.size() != seedingAlgo.size() )
83  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
84  << " WARNING : maxD0 does not have the proper size "
85  << std::endl;
86 
87  maxZ0 = conf.getParameter<std::vector<double> >("maxZ0");
88  if ( maxZ0.size() != seedingAlgo.size() )
89  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
90  << " WARNING : maxZ0 does not have the proper size "
91  << std::endl;
92 
93  // The name of the hit producer
94  hitProducer = conf.getParameter<edm::InputTag>("HitProducer");
95 
96  // The cuts for seed cleaning
97  seedCleaning = conf.getParameter<bool>("seedCleaning");
98 
99  // Number of hits needed for a seed
100  numberOfHits = conf.getParameter<std::vector<unsigned int> >("numberOfHits");
101  if ( numberOfHits.size() != seedingAlgo.size() )
102  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
103  << " WARNING : numberOfHits does not have the proper size "
104  << std::endl;
105 
106  // Seeding based on muons
107  selectMuons = conf.getParameter<bool>("selectMuons");
108 
109  // Layers
110  newSyntax = conf.getParameter<bool>("newSyntax");
111  if (newSyntax) {
112  layerList = conf.getParameter<std::vector<std::string> >("layerList");
113  // (AG) for (unsigned i=0; i<layerList.size();i++) std::cout << "------- Layers = " << layerList[i] << std::endl;
114  } else {
115  // TO BE DELETED (AG)
117  conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectorNumber");
118  if ( firstHitSubDetectorNumber.size() != seedingAlgo.size() )
119  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
120  << " WARNING : firstHitSubDetectorNumber does not have the proper size "
121  << std::endl;
122 
123  std::vector<unsigned int> firstSubDets =
124  conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectors");
125  unsigned isub1 = 0;
126  unsigned check1 = 0;
127  firstHitSubDetectors.resize(seedingAlgo.size());
128  for ( unsigned ialgo=0; ialgo<firstHitSubDetectorNumber.size(); ++ialgo ) {
129  check1 += firstHitSubDetectorNumber[ialgo];
130  for ( unsigned idet=0; idet<firstHitSubDetectorNumber[ialgo]; ++idet ) {
131  firstHitSubDetectors[ialgo].push_back(firstSubDets[isub1++]);
132  }
133  }
134  if ( firstSubDets.size() != check1 )
135  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
136  << " WARNING : firstHitSubDetectors does not have the proper size (should be " << check1 << ")"
137  << std::endl;
138 
139 
141  conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectorNumber");
142  if ( secondHitSubDetectorNumber.size() != seedingAlgo.size() )
143  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
144  << " WARNING : secondHitSubDetectorNumber does not have the proper size "
145  << std::endl;
146 
147  std::vector<unsigned int> secondSubDets =
148  conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectors");
149  unsigned isub2 = 0;
150  unsigned check2 = 0;
151  secondHitSubDetectors.resize(seedingAlgo.size());
152  for ( unsigned ialgo=0; ialgo<secondHitSubDetectorNumber.size(); ++ialgo ) {
153  check2 += secondHitSubDetectorNumber[ialgo];
154  for ( unsigned idet=0; idet<secondHitSubDetectorNumber[ialgo]; ++idet ) {
155  secondHitSubDetectors[ialgo].push_back(secondSubDets[isub2++]);
156  }
157  }
158  if ( secondSubDets.size() != check2 )
159  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
160  << " WARNING : secondHitSubDetectors does not have the proper size (should be " << check2 << ")"
161  << std::endl;
162 
164  conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectorNumber");
165  if ( thirdHitSubDetectorNumber.size() != seedingAlgo.size() )
166  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
167  << " WARNING : thirdHitSubDetectorNumber does not have the proper size "
168  << std::endl;
169 
170  std::vector<unsigned int> thirdSubDets =
171  conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectors");
172  unsigned isub3 = 0;
173  unsigned check3 = 0;
174  thirdHitSubDetectors.resize(seedingAlgo.size());
175  for ( unsigned ialgo=0; ialgo<thirdHitSubDetectorNumber.size(); ++ialgo ) {
176  check3 += thirdHitSubDetectorNumber[ialgo];
177  for ( unsigned idet=0; idet<thirdHitSubDetectorNumber[ialgo]; ++idet ) {
178  thirdHitSubDetectors[ialgo].push_back(thirdSubDets[isub3++]);
179  }
180  }
181  if ( thirdSubDets.size() != check3 )
182  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
183  << " WARNING : thirdHitSubDetectors does not have the proper size (should be " << check3 << ")"
184  << std::endl;
185 
186  originRadius = conf.getParameter<std::vector<double> >("originRadius");
187  if ( originRadius.size() != seedingAlgo.size() )
188  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
189  << " WARNING : originRadius does not have the proper size "
190  << std::endl;
191 
192  originHalfLength = conf.getParameter<std::vector<double> >("originHalfLength");
193  if ( originHalfLength.size() != seedingAlgo.size() )
194  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
195  << " WARNING : originHalfLength does not have the proper size "
196  << std::endl;
197 
198  originpTMin = conf.getParameter<std::vector<double> >("originpTMin");
199  if ( originpTMin.size() != seedingAlgo.size() )
200  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
201  << " WARNING : originpTMin does not have the proper size "
202  << std::endl;
203 
204  primaryVertices = conf.getParameter<std::vector<edm::InputTag> >("primaryVertices");
205  if ( primaryVertices.size() != seedingAlgo.size() )
206  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
207  << " WARNING : primaryVertices does not have the proper size "
208  << std::endl;
209 
210  zVertexConstraint = conf.getParameter<std::vector<double> >("zVertexConstraint");
211  if ( zVertexConstraint.size() != seedingAlgo.size() )
212  throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
213  << " WARNING : zVertexConstraint does not have the proper size "
214  << std::endl;
215  }
216 
217 }
218 
219 
220 // Virtual destructor needed.
222 
223  if(thePropagator) delete thePropagator;
224 
225  // do nothing
226 #ifdef FAMOS_DEBUG
227  std::cout << "TrajectorySeedProducer destructed" << std::endl;
228 #endif
229 
230 }
231 
232 void
234 
235  //services
236  // es.get<TrackerRecoGeometryRecord>().get(theGeomSearchTracker);
237 
241 
242 
243  es.get<IdealMagneticFieldRecord>().get(magField);
244  es.get<TrackerDigiGeometryRecord>().get(geometry);
245  es.get<MagneticFieldMapRecord>().get(magFieldMap);
246 
247  theMagField = &(*magField);
248  theGeometry = &(*geometry);
249  theFieldMap = &(*magFieldMap);
250 
252 
253  const GlobalPoint g(0.,0.,0.);
254 
255 }
256 
257  // Functions that gets called by framework every event
258 void
260 
261 
262  // if( seedingAlgo[0] == "FourthPixelLessPairs") std::cout << "Seed producer in 4th iteration " << std::endl;
263 
264 #ifdef FAMOS_DEBUG
265  std::cout << "################################################################" << std::endl;
266  std::cout << " TrajectorySeedProducer produce init " << std::endl;
267 #endif
268 
269  //Retrieve tracker topology from geometry
271  es.get<IdealGeometryRecord>().get(tTopoHand);
272  const TrackerTopology *tTopo=tTopoHand.product();
273 
274 
275  unsigned nSimTracks = 0;
276  unsigned nTracksWithHits = 0;
277  unsigned nTracksWithPT = 0;
278  unsigned nTracksWithD0Z0 = 0;
279  // unsigned nTrackCandidates = 0;
280  PTrajectoryStateOnDet initialState;
281 
282  // Output
283  std::vector<TrajectorySeedCollection*>
284  output(seedingAlgo.size(),static_cast<TrajectorySeedCollection*>(0));
285  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
286  // std::auto_ptr<TrajectorySeedCollection> p(new TrajectorySeedCollection );
287  output[ialgo] = new TrajectorySeedCollection;
288  }
289 
290  // Beam spot
291  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
292  e.getByLabel(theBeamSpot,recoBeamSpotHandle);
293  math::XYZPoint BSPosition_ = recoBeamSpotHandle->position();
294 
295  //not used anymore. take the value from the py
296 
297  //double sigmaZ=recoBeamSpotHandle->sigmaZ();
298  //double sigmaZ0Error=recoBeamSpotHandle->sigmaZ0Error();
299  //double sigmaz0=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
300  x0 = BSPosition_.X();
301  y0 = BSPosition_.Y();
302  z0 = BSPosition_.Z();
303 
304  // SimTracks and SimVertices
306  e.getByLabel("famosSimHits",theSimTracks);
307 
309  e.getByLabel("famosSimHits",theSimVtx);
310 
311 #ifdef FAMOS_DEBUG
312  std::cout << " Step A: SimTracks found " << theSimTracks->size() << std::endl;
313 #endif
314 
315  // edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
317  e.getByLabel(hitProducer, theGSRecHits);
318 
319  // No tracking attempted if no hits (but put an empty collection in the event)!
320 #ifdef FAMOS_DEBUG
321  std::cout << " Step B: Full GS RecHits found " << theGSRecHits->size() << std::endl;
322 #endif
323  if(theGSRecHits->size() == 0) {
324  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
325  std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
326  e.put(p,seedingAlgo[ialgo]);
327  }
328  return;
329  }
330 
331  // Primary vertices
332  vertices = std::vector<const reco::VertexCollection*>
333  (seedingAlgo.size(),static_cast<const reco::VertexCollection*>(0));
334  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
335  //PAT Attempt!!!!
336 
337  //originHalfLength[ialgo] = 3.*sigmaz0; // Overrides the configuration
339  bool isVertexCollection = e.getByLabel(primaryVertices[ialgo],aHandle);
340  if (!isVertexCollection ) continue;
341  vertices[ialgo] = &(*aHandle);
342  }
343 
344 #ifdef FAMOS_DEBUG
345  std::cout << " Step C: Loop over the RecHits, track by track " << std::endl;
346 #endif
347 
348  // The vector of simTrack Id's carrying GSRecHits
349  const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
350 
351  // loop over SimTrack Id's
352  for ( unsigned tkId=0; tkId != theSimTrackIds.size(); ++tkId ) {
353 
354 #ifdef FAMOS_DEBUG
355  std::cout << "Track number " << tkId << std::endl;
356 #endif
357 
358  ++nSimTracks;
359  unsigned simTrackId = theSimTrackIds[tkId];
360  const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
361 #ifdef FAMOS_DEBUG
362  std::cout << "Pt = " << std::sqrt(theSimTrack.momentum().Perp2())
363  << " eta " << theSimTrack.momentum().Eta()
364  << " pdg ID " << theSimTrack.type()
365  << std::endl;
366 #endif
367 
368  // Select only muons, if requested
369  if (selectMuons && abs(theSimTrack.type()) != 13) continue;
370 
371  // Check that the sim track comes from the main vertex (loose cut)
372  int vertexIndex = theSimTrack.vertIndex();
373  const SimVertex& theSimVertex = (*theSimVtx)[vertexIndex];
374 #ifdef FAMOS_DEBUG
375  std::cout << " o SimTrack " << theSimTrack << std::endl;
376  std::cout << " o SimVertex " << theSimVertex << std::endl;
377 #endif
378 
379  BaseParticlePropagator theParticle =
381  RawParticle(XYZTLorentzVector(theSimTrack.momentum().px(),
382  theSimTrack.momentum().py(),
383  theSimTrack.momentum().pz(),
384  theSimTrack.momentum().e()),
385  XYZTLorentzVector(theSimVertex.position().x(),
386  theSimVertex.position().y(),
387  theSimVertex.position().z(),
388  theSimVertex.position().t())),
389  0.,0.,4.);
390  theParticle.setCharge((*theSimTracks)[simTrackId].charge());
391 
392  SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
393  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
394  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd = theRecHitRange.second;
399 
400  // Check the number of layers crossed
401  unsigned numberOfRecHits = 0;
402  TrackerRecHit previousHit, currentHit;
403  for ( iterRecHit = theRecHitRangeIteratorBegin;
404  iterRecHit != theRecHitRangeIteratorEnd;
405  ++iterRecHit) {
406  previousHit = currentHit;
407  currentHit = TrackerRecHit(&(*iterRecHit),theGeometry,tTopo);
408  if ( currentHit.isOnTheSameLayer(previousHit) ) continue;
409  ++numberOfRecHits;
410  if ( numberOfRecHits == absMinRecHits ) break;
411  }
412 
413  // Loop on the successive seedings
414  for ( unsigned int ialgo = 0; ialgo < seedingAlgo.size(); ++ialgo ) {
415 
416 #ifdef FAMOS_DEBUG
417  std::cout << "Algo " << seedingAlgo[ialgo] << std::endl;
418 #endif
419 
420  // Request a minimum number of RecHits for the track to give a seed.
421 #ifdef FAMOS_DEBUG
422  std::cout << "The number of RecHits = " << numberOfRecHits << std::endl;
423 #endif
424  if ( numberOfRecHits < minRecHits[ialgo] ) continue;
425  ++nTracksWithHits;
426 
427  // Request a minimum pT for the sim track
428  if ( theSimTrack.momentum().Perp2() < pTMin[ialgo] ) continue;
429  ++nTracksWithPT;
430 
431  // Cut on sim track impact parameters
432  if ( theParticle.xyImpactParameter(x0,y0) > maxD0[ialgo] ) continue;
433  if ( fabs( theParticle.zImpactParameter(x0,y0) - z0 ) > maxZ0[ialgo] ) continue;
434  ++nTracksWithD0Z0;
435 
436  std::vector<TrackerRecHit >
437  theSeedHits(numberOfHits[ialgo],
438  static_cast<TrackerRecHit >(TrackerRecHit()));
439  TrackerRecHit& theSeedHits0 = theSeedHits[0];
440  TrackerRecHit& theSeedHits1 = theSeedHits[1];
441  TrackerRecHit& theSeedHits2 = theSeedHits[2];
442  bool compatible = false;
443  for ( iterRecHit1 = theRecHitRangeIteratorBegin; iterRecHit1 != theRecHitRangeIteratorEnd; ++iterRecHit1) {
444  theSeedHits[0] = TrackerRecHit(&(*iterRecHit1),theGeometry,tTopo);
445 #ifdef FAMOS_DEBUG
446  std::cout << "The first hit position = " << theSeedHits0.globalPosition() << std::endl;
447  std::cout << "The first hit subDetId = " << theSeedHits0.subDetId() << std::endl;
448  std::cout << "The first hit layer = " << theSeedHits0.layerNumber() << std::endl;
449 #endif
450 
451  // Check if inside the requested detectors
452  bool isInside = true;
453  if (!selectMuons) {
454  if (newSyntax)
455  isInside = true; // AG placeholder
456  else
457  isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
458  // bool isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
459  if ( isInside ) continue;
460  }
461 
462  // Check if on requested detectors
463  // bool isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo]);
464  bool isOndet = true;
465  if (!selectMuons) {
466  if (newSyntax)
467  isOndet = theSeedHits0.isOnRequestedDet(layerList);
468  else
469  isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
470  // bool isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
471  // if ( !isOndet ) break;
472  if ( !isOndet ) continue;
473  }
474 
475 #ifdef FAMOS_DEBUG
476  std::cout << "Apparently the first hit is on the requested detector! " << std::endl;
477 #endif
478  for ( iterRecHit2 = iterRecHit1+1; iterRecHit2 != theRecHitRangeIteratorEnd; ++iterRecHit2) {
479  theSeedHits[1] = TrackerRecHit(&(*iterRecHit2),theGeometry,tTopo);
480 #ifdef FAMOS_DEBUG
481  std::cout << "The second hit position = " << theSeedHits1.globalPosition() << std::endl;
482  std::cout << "The second hit subDetId = " << theSeedHits1.subDetId() << std::endl;
483  std::cout << "The second hit layer = " << theSeedHits1.layerNumber() << std::endl;
484 #endif
485 
486  if (!selectMuons) {
487  // Check if inside the requested detectors
488  if (newSyntax)
489  isInside = true; // AG placeholder
490  else
491  isInside = theSeedHits1.subDetId() < secondHitSubDetectors[ialgo][0];
492  if ( isInside ) continue;
493 
494  // Check if on requested detectors
495  if (newSyntax)
496  isOndet = theSeedHits1.isOnRequestedDet(layerList);
497  else
498  isOndet = theSeedHits1.isOnRequestedDet(secondHitSubDetectors[ialgo], seedingAlgo[ialgo]);
499  if ( !isOndet ) break;
500  }
501 
502  // Check if on the same layer as previous hit
503  if ( theSeedHits1.isOnTheSameLayer(theSeedHits0) ) continue;
504 
505 #ifdef FAMOS_DEBUG
506  std::cout << "Apparently the second hit is on the requested detector! " << std::endl;
507 #endif
508  GlobalPoint gpos1 = theSeedHits0.globalPosition();
509  GlobalPoint gpos2 = theSeedHits1.globalPosition();
510  bool forward = theSeedHits0.isForward();
511  double error = std::sqrt(theSeedHits0.largerError()+theSeedHits1.largerError());
512  // compatible = compatibleWithVertex(gpos1,gpos2,ialgo);
513  //added out of desperation
514  if(seedingAlgo[0] == "PixelLess" || seedingAlgo[0] == "TobTecLayerPairs"){
515  compatible = true;
516  //std::cout << "Algo " << seedingAlgo[0] << "Det/layer = " << theSeedHits0.subDetId() << "/" << theSeedHits0.layerNumber() << std::endl;
517  } else {
518  compatible = compatibleWithBeamAxis(gpos1,gpos2,error,forward,ialgo);
519  }
520 
521 #ifdef FAMOS_DEBUG
522  std::cout << "Algo" << seedingAlgo[0] << "\t Are the two hits compatible with the PV? " << compatible << std::endl;
523 #endif
524 
525  if (!selectMuons) {
526  // Check if the pair is on the requested dets
527  if ( numberOfHits[ialgo] == 2 ) {
528 
529  if ( seedingAlgo[0] == "ThirdMixedPairs" ){
530  compatible = compatible && theSeedHits[0].makesAPairWith3rd(theSeedHits[1]);
531  } else {
532  compatible = compatible && theSeedHits[0].makesAPairWith(theSeedHits[1]);
533  //check
534  /*
535  if((seedingAlgo[0] == "PixelLess" || seedingAlgo[0] == "TobTecLayerPairs") && !compatible)
536  std::cout << "NOT Compatible " << seedingAlgo[0]
537  << "Hit 1 Det/layer/ring = " << theSeedHits0.subDetId() << "/" << theSeedHits0.layerNumber() << "/" << theSeedHits0.ringNumber()
538  << "\tHit 2 Det/layer/ring = " << theSeedHits1.subDetId() << "/" << theSeedHits1.layerNumber() << "/" << theSeedHits1.ringNumber() << std::endl;
539  */
540  }
541  }
542  }
543 
544  // Reject non suited pairs
545  if ( !compatible ) continue;
546 
547 #ifdef FAMOS_DEBUG
548  std::cout << "Pair kept! " << std::endl;
549 #endif
550 
551  // Leave here if only two hits are required.
552  if ( numberOfHits[ialgo] == 2 ) break;
553 
554  compatible = false;
555  // Check if there is a third satisfying hit otherwise
556  for ( iterRecHit3 = iterRecHit2+1; iterRecHit3 != theRecHitRangeIteratorEnd; ++iterRecHit3) {
557  theSeedHits[2] = TrackerRecHit(&(*iterRecHit3),theGeometry,tTopo);
558 #ifdef FAMOS_DEBUG
559  std::cout << "The third hit position = " << theSeedHits2.globalPosition() << std::endl;
560  std::cout << "The third hit subDetId = " << theSeedHits2.subDetId() << std::endl;
561  std::cout << "The third hit layer = " << theSeedHits2.layerNumber() << std::endl;
562 #endif
563 
564  if (!selectMuons) {
565  // Check if inside the requested detectors
566  if (newSyntax)
567  isInside = true; // AG placeholder
568  else
569  isInside = theSeedHits2.subDetId() < thirdHitSubDetectors[ialgo][0];
570  if ( isInside ) continue;
571 
572  // Check if on requested detectors
573  if (newSyntax)
574  isOndet = theSeedHits2.isOnRequestedDet(layerList);
575  else
576  isOndet = theSeedHits2.isOnRequestedDet(thirdHitSubDetectors[ialgo], seedingAlgo[ialgo]);
577  // if ( !isOndet ) break;
578  if ( !isOndet ) continue;
579  }
580 
581  // Check if on the same layer as previous hit
582  compatible = !(theSeedHits2.isOnTheSameLayer(theSeedHits1));
583 
584  // Check if the triplet is on the requested det combination
585  if (!selectMuons) compatible = compatible && theSeedHits[0].makesATripletWith(theSeedHits[1],theSeedHits[2]);
586 
587 #ifdef FAMOS_DEBUG
588  if ( compatible )
589  std::cout << "Apparently the third hit is on the requested detector! " << std::endl;
590 #endif
591 
592  if ( compatible ) break;
593 
594  }
595 
596  if ( compatible ) break;
597 
598  }
599 
600  if ( compatible ) break;
601 
602  }
603 
604  // There is no compatible seed for this track with this seeding algorithm
605  // Go to next algo
606  if ( !compatible ) continue;
607 #ifdef FAMOS_DEBUG
608  if ( compatible )
609  std::cout << "@@@ There is at least a compatible seed" << std::endl;
610  else
611  std::cout << "@@@ There is no compatible seed" << std::endl;
612 #endif
613 
614 
615 #ifdef FAMOS_DEBUG
616  std::cout << "Preparing to create the TrajectorySeed" << std::endl;
617 #endif
618  // The seed is validated -> include in the collection
619  // 1) Create the vector of RecHits
621  for ( unsigned ih=0; ih<theSeedHits.size(); ++ih ) {
622  TrackingRecHit* aTrackingRecHit = theSeedHits[ih].hit()->clone();
623  recHits.push_back(aTrackingRecHit);
624  }
625 #ifdef FAMOS_DEBUG
626  std::cout << "with " << recHits.size() << " hits." << std::endl;
627 #endif
628 
629  // 2) Create the initial state
630  // a) origin vertex
631  GlobalPoint position((*theSimVtx)[vertexIndex].position().x(),
632  (*theSimVtx)[vertexIndex].position().y(),
633  (*theSimVtx)[vertexIndex].position().z());
634 
635  // b) initial momentum
636  GlobalVector momentum( (*theSimTracks)[simTrackId].momentum().x() ,
637  (*theSimTracks)[simTrackId].momentum().y() ,
638  (*theSimTracks)[simTrackId].momentum().z() );
639  // c) electric charge
640  float charge = (*theSimTracks)[simTrackId].charge();
641  // -> inital parameters
642  GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
643  // -> large initial errors
644  AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();
645  // errorMatrix = errorMatrix * 10;
646 
647  //this line help the fit succeed in the case of pixelless tracks (4th and 5th iteration)
648  //for the future: probably the best thing is to use the mini-kalmanFilter
649  if(theSeedHits0.subDetId() !=1 || theSeedHits0.subDetId() !=2) errorMatrix = errorMatrix * 0.0000001;
650 
651 
652 
653 #ifdef FAMOS_DEBUG
654  std::cout << "TrajectorySeedProducer: SimTrack parameters " << std::endl;
655  std::cout << "\t\t pT = " << (*theSimTracks)[simTrackId].momentum().Pt() << std::endl;
656  std::cout << "\t\t eta = " << (*theSimTracks)[simTrackId].momentum().Eta() << std::endl;
657  std::cout << "\t\t phi = " << (*theSimTracks)[simTrackId].momentum().Phi() << std::endl;
658  std::cout << "TrajectorySeedProducer: AlgebraicSymMatrix " << errorMatrix << std::endl;
659 #endif
660  CurvilinearTrajectoryError initialError(errorMatrix);
661  // -> initial state
662  FreeTrajectoryState initialFTS(initialParams, initialError);
663 #ifdef FAMOS_DEBUG
664  std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
665 #endif
666  // const GeomDetUnit* initialLayer = theGeometry->idToDetUnit( recHits.front().geographicalId() );
667  const GeomDet* initialLayer = theGeometry->idToDet( recHits.front().geographicalId() );
668 
669  //this is wrong because the FTS is defined at vertex, and it need to be properly propagated.
670  // const TrajectoryStateOnSurface initialTSOS(initialFTS, initialLayer->surface());
671 
672  const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
673  if (!initialTSOS.isValid()) continue;
674 
675 #ifdef FAMOS_DEBUG
676  std::cout << "TrajectorySeedProducer: TSOS global momentum " << initialTSOS.globalMomentum() << std::endl;
677  std::cout << "\t\t\tpT = " << initialTSOS.globalMomentum().perp() << std::endl;
678  std::cout << "\t\t\teta = " << initialTSOS.globalMomentum().eta() << std::endl;
679  std::cout << "\t\t\tphi = " << initialTSOS.globalMomentum().phi() << std::endl;
680  std::cout << "TrajectorySeedProducer: TSOS local momentum " << initialTSOS.localMomentum() << std::endl;
681  std::cout << "TrajectorySeedProducer: TSOS local error " << initialTSOS.localError().positionError() << std::endl;
682  std::cout << "TrajectorySeedProducer: TSOS local error matrix " << initialTSOS.localError().matrix() << std::endl;
683  std::cout << "TrajectorySeedProducer: TSOS surface side " << initialTSOS.surfaceSide() << std::endl;
684 #endif
685  stateOnDet(initialTSOS,
686  recHits.front().geographicalId().rawId(),
687  initialState);
688  // Create a new Trajectory Seed
689  output[ialgo]->push_back(TrajectorySeed(initialState, recHits, alongMomentum));
690 #ifdef FAMOS_DEBUG
691  std::cout << "Trajectory seed created ! " << std::endl;
692 #endif
693  break;
694  // End of the loop over seeding algorithms
695  }
696  // End on the loop over simtracks
697  }
698 
699  for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
700  std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
701  e.put(p,seedingAlgo[ialgo]);
702  }
703 
704 }
705 
706 // This is a copy of a method in
707 // TrackingTools/TrajectoryState/src/TrajectoryStateTransform.cc
708 // but it does not return a pointer (thus avoiding a memory leak)
709 // In addition, it's also CPU more efficient, because
710 // ts.localError().matrix() is not copied
711 void
713  unsigned int detid,
714  PTrajectoryStateOnDet& pts) const
715 {
716 
717  const AlgebraicSymMatrix55& m = ts.localError().matrix();
718 
719  int dim = 5;
720 
721  float localErrors[15];
722  int k = 0;
723  for (int i=0; i<dim; ++i) {
724  for (int j=0; j<=i; ++j) {
725  localErrors[k++] = m(i,j);
726  }
727  }
728  int surfaceSide = static_cast<int>(ts.surfaceSide());
729 
731  localErrors, detid,
732  surfaceSide);
733 }
734 
735 bool
737  GlobalPoint& gpos2,
738  double error,
739  bool forward,
740  unsigned algo) const {
741 
742  if ( !seedCleaning ) return true;
743 
744  // The hits 1 and 2 positions, in HepLorentzVector's
745  XYZTLorentzVector thePos1(gpos1.x(),gpos1.y(),gpos1.z(),0.);
746  XYZTLorentzVector thePos2(gpos2.x(),gpos2.y(),gpos2.z(),0.);
747 #ifdef FAMOS_DEBUG
748  std::cout << "ThePos1 = " << thePos1 << std::endl;
749  std::cout << "ThePos2 = " << thePos2 << std::endl;
750 #endif
751 
752 
753  // Create new particles that pass through the second hit with pT = ptMin
754  // and charge = +/-1
755 
756  // The momentum direction is by default joining the two hits
757  XYZTLorentzVector theMom2 = (thePos2-thePos1);
758 
759  // The corresponding RawParticle, with an (irrelevant) electric charge
760  // (The charge is determined in the next step)
761  ParticlePropagator myPart(theMom2,thePos2,1.,theFieldMap);
762 
766  bool intersect = myPart.propagateToBeamCylinder(thePos1,originRadius[algo]*1.);
767  if ( !intersect ) return false;
768 
769 #ifdef FAMOS_DEBUG
770  std::cout << "MyPart R = " << myPart.R() << "\t Z = " << myPart.Z()
771  << "\t pT = " << myPart.Pt() << std::endl;
772 #endif
773 
774  // Check if the constraints are satisfied
775  // 1. pT at cylinder with radius originRadius
776  if ( myPart.Pt() < originpTMin[algo] ) return false;
777 
778  // 2. Z compatible with beam spot size
779  if ( fabs(myPart.Z()-z0) > originHalfLength[algo] ) return false;
780 
781  // 3. Z compatible with one of the primary vertices (always the case if no primary vertex)
782  const reco::VertexCollection* theVertices = vertices[algo];
783  if (!theVertices) return true;
784  unsigned nVertices = theVertices->size();
785  if ( !nVertices || zVertexConstraint[algo] < 0. ) return true;
786  // Radii of the two hits with respect to the beam spot position
787  double R1 = std::sqrt ( (thePos1.X()-x0)*(thePos1.X()-x0)
788  + (thePos1.Y()-y0)*(thePos1.Y()-y0) );
789  double R2 = std::sqrt ( (thePos2.X()-x0)*(thePos2.X()-x0)
790  + (thePos2.Y()-y0)*(thePos2.Y()-y0) );
791  // Loop on primary vertices
792  for ( unsigned iv=0; iv<nVertices; ++iv ) {
793  // Z position of the primary vertex
794  double zV = (*theVertices)[iv].z();
795  // Constraints on the inner hit
796  double checkRZ1 = forward ?
797  (thePos1.Z()-zV+zVertexConstraint[algo]) / (thePos2.Z()-zV+zVertexConstraint[algo]) * R2 :
798  -zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV+zVertexConstraint[algo]);
799  double checkRZ2 = forward ?
800  (thePos1.Z()-zV-zVertexConstraint[algo])/(thePos2.Z()-zV-zVertexConstraint[algo]) * R2 :
801  +zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV-zVertexConstraint[algo]);
802  double checkRZmin = std::min(checkRZ1,checkRZ2)-3.*error;
803  double checkRZmax = std::max(checkRZ1,checkRZ2)+3.*error;
804  // Check if the innerhit is within bounds
805  bool compat = forward ?
806  checkRZmin < R1 && R1 < checkRZmax :
807  checkRZmin < thePos1.Z()-zV && thePos1.Z()-zV < checkRZmax;
808  // If it is, just return ok
809  if ( compat ) return compat;
810  }
811  // Otherwise, return not ok
812  return false;
813 
814 }
815 
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:52
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
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.)
#define abs(x)
Definition: mlp_lapack.h:159
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
#define min(a, b)
Definition: mlp_lapack.h:161
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< std::string > layerList
std::vector< double > originpTMin
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:97
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:45
LocalVector localMomentum() const
PropagatorWithMaterial * thePropagator
double R() const
vertex radius
Definition: RawParticle.h:277
static const double pts[33]
Definition: Constants.h:31
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:45
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:94
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
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:361
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:91
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
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:82
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:79
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:36
double xyImpactParameter(double x0=0., double y0=0.) const
Transverse impact parameter.
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:15
const MagneticField * theMagField