CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SimpleCosmicBONSeeder.cc
Go to the documentation of this file.
1 //
2 // Package: RecoTracker/TkSeedGenerator
3 // Class: GlobalPixelLessSeedGenerator
4 //
5 // From CosmicSeedGenerator+SimpleCosmicBONSeeder, with changes by Giovanni
6 // to seed Cosmics with B != 0
7 
12 
13 #include <numeric>
14 
15 using namespace std;
17  conf_(conf),
18  theLsb(conf.getParameter<edm::ParameterSet>("TripletsPSet")),
19  writeTriplets_(conf.getParameter<bool>("writeTriplets")),
20  seedOnMiddle_(conf.existsAs<bool>("seedOnMiddle") ? conf.getParameter<bool>("seedOnMiddle") : false),
21  rescaleError_(conf.existsAs<double>("rescaleError") ? conf.getParameter<double>("rescaleError") : 1.0),
22  tripletsVerbosity_(conf.getParameter<edm::ParameterSet>("TripletsPSet").getUntrackedParameter<uint32_t>("debugLevel",0)),
23  seedVerbosity_(conf.getUntrackedParameter<uint32_t>("seedDebugLevel",0)),
24  helixVerbosity_(conf.getUntrackedParameter<uint32_t>("helixDebugLevel",0)),
25  check_(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")),
26  maxTriplets_(conf.getParameter<int32_t>("maxTriplets")),
27  maxSeeds_(conf.getParameter<int32_t>("maxSeeds"))
28 {
29  edm::ParameterSet regionConf = conf_.getParameter<edm::ParameterSet>("RegionPSet");
30  float ptmin = regionConf.getParameter<double>("ptMin");
31  float originradius = regionConf.getParameter<double>("originRadius");
32  float halflength = regionConf.getParameter<double>("originHalfLength");
33  float originz = regionConf.getParameter<double>("originZPosition");
34  region_ = GlobalTrackingRegion(ptmin, originradius, halflength, originz);
35  pMin_ = regionConf.getParameter<double>("pMin");
36 
37  builderName = conf_.getParameter<std::string>("TTRHBuilder");
38 
39  //***top-bottom
40  positiveYOnly=conf_.getParameter<bool>("PositiveYOnly");
41  negativeYOnly=conf_.getParameter<bool>("NegativeYOnly");
42  //***
43 
44  produces<TrajectorySeedCollection>();
45  if (writeTriplets_) produces<edm::OwnVector<TrackingRecHit> >("cosmicTriplets");
46 
47  layerTripletNames_ = conf_.getParameter<edm::ParameterSet>("TripletsPSet").getParameter<std::vector<std::string> >("layerList");
48 
49  if (conf.existsAs<edm::ParameterSet>("ClusterChargeCheck")) {
50  edm::ParameterSet cccc = conf.getParameter<edm::ParameterSet>("ClusterChargeCheck");
51  checkCharge_ = cccc.getParameter<bool>("checkCharge");
52  matchedRecHitUsesAnd_ = cccc.getParameter<bool>("matchedRecHitsUseAnd");
53  chargeThresholds_.resize(7,0);
54  edm::ParameterSet ccct = cccc.getParameter<edm::ParameterSet>("Thresholds");
56  chargeThresholds_[StripSubdetector::TID] = ccct.getParameter<int32_t>("TID");
57  chargeThresholds_[StripSubdetector::TOB] = ccct.getParameter<int32_t>("TOB");
58  chargeThresholds_[StripSubdetector::TEC] = ccct.getParameter<int32_t>("TEC");
59  } else {
60  checkCharge_ = false;
61  }
62  if (conf.existsAs<edm::ParameterSet>("HitsPerModuleCheck")) {
63  edm::ParameterSet hpmcc = conf.getParameter<edm::ParameterSet>("HitsPerModuleCheck");
64  checkMaxHitsPerModule_ = hpmcc.getParameter<bool>("checkHitsPerModule");
66  edm::ParameterSet hpmct = hpmcc.getParameter<edm::ParameterSet>("Thresholds");
67  maxHitsPerModule_[StripSubdetector::TIB] = hpmct.getParameter<int32_t>("TIB");
68  maxHitsPerModule_[StripSubdetector::TID] = hpmct.getParameter<int32_t>("TID");
69  maxHitsPerModule_[StripSubdetector::TOB] = hpmct.getParameter<int32_t>("TOB");
70  maxHitsPerModule_[StripSubdetector::TEC] = hpmct.getParameter<int32_t>("TEC");
71  } else {
72  checkMaxHitsPerModule_ = false;
73  }
75  goodHitsPerSeed_ = conf.getParameter<int32_t>("minimumGoodHitsInSeed");
76  } else {
77  goodHitsPerSeed_ = 0;
78  }
79 
80 }
81 
82 // Functions that gets called by framework every event
84 {
85  std::auto_ptr<TrajectorySeedCollection> output(new TrajectorySeedCollection());
86  std::auto_ptr<edm::OwnVector<TrackingRecHit> > outtriplets(new edm::OwnVector<TrackingRecHit>());
87 
89  if (magfield->inTesla(GlobalPoint(0,0,0)).mag() > 0.01) {
90  size_t clustsOrZero = check_.tooManyClusters(ev);
91  if (clustsOrZero) {
92  edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
93  } else {
94  init(es);
95  bool tripletsOk = triplets(ev,es);
96  if (tripletsOk) {
97 
98  bool seedsOk = seeds(*output,es);
99  if (!seedsOk) { }
100 
101  if (writeTriplets_) {
102  for (OrderedHitTriplets::const_iterator it = hitTriplets.begin(); it != hitTriplets.end(); ++it) {
103  const TrackingRecHit * hit1 = it->inner()->hit();
104  const TrackingRecHit * hit2 = it->middle()->hit();
105  const TrackingRecHit * hit3 = it->outer()->hit();
106  outtriplets->push_back(hit1->clone());
107  outtriplets->push_back(hit2->clone());
108  outtriplets->push_back(hit3->clone());
109  }
110  }
111  }
112  done();
113  }
114  }
115 
116  if (writeTriplets_) {
117  ev.put(outtriplets, "cosmicTriplets");
118  }
119  ev.put(output);
120 }
121 
122 void
124 {
125  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
127 
128  // FIXME: these should come from ES too!!
131  theUpdator = new KFUpdator();
132 
133 }
134 
136  bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const {
137  //FIXME: inner gives a SEGV
138 #if 0
139  //const TransientTrackingRecHit::ConstRecHitPointer &ihit1 = trip1.inner();
140  //const TransientTrackingRecHit::ConstRecHitPointer &ihit2 = trip2.inner();
143  const TransientTrackingRecHit::ConstRecHitPointer &ohit1 = trip1.outer();
144  const TransientTrackingRecHit::ConstRecHitPointer &ohit2 = trip2.outer();
145 #endif
150  float iy1 = ihit1->globalPosition().y();
151  float oy1 = ohit1->globalPosition().y();
152  float iy2 = ihit2->globalPosition().y();
153  float oy2 = ohit2->globalPosition().y();
154  if (oy1 - iy1 > 0) { // 1 Downgoing
155  if (oy2 - iy2 > 0) { // 2 Downgoing
156  // sort by inner, or by outer if inners are the same
157  return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
158  } else return true; // else prefer downgoing
159  } else if (oy2 - iy2 > 0) {
160  return false; // prefer downgoing
161  } else { // both upgoing
162  // sort by inner, or by outer
163  return (iy1 != iy2 ? (iy1 < iy2) : (oy1 < oy2));
164  }
165  }
166 };
167 
169  using namespace ctfseeding;
170 
171  hitTriplets.clear();
172  hitTriplets.reserve(0);
173  SeedingLayerSets lss = theLsb.layers(es);
174  SeedingLayerSets::const_iterator iLss;
175 
176  double minRho = region_.ptMin() / ( 0.003 * magfield->inTesla(GlobalPoint(0,0,0)).z() );
177 
178  for (iLss = lss.begin(); iLss != lss.end(); iLss++){
179  SeedingLayers ls = *iLss;
180  if (ls.size() != 3){
181  throw cms::Exception("CtfSpecialSeedGenerator") << "You are using " << ls.size() <<" layers in set instead of 3 ";
182  }
183 
185  std::vector<SeedingHit> innerHits = region_.hits(e, es, &ls[0]);
186  std::vector<SeedingHit> middleHits = region_.hits(e, es, &ls[1]);
187  std::vector<SeedingHit> outerHits = region_.hits(e, es, &ls[2]);
188  std::vector<SeedingHit>::const_iterator iOuterHit,iMiddleHit,iInnerHit;
189 
190  if (tripletsVerbosity_ > 0) {
191  std::cout << "GenericTripletGenerator iLss = " << layerTripletNames_[iLss - lss.begin()]
192  << " (" << (iLss - lss.begin()) << "): # = "
193  << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size() << std::endl;
194  }
195 
198  std::vector<TTRH> innerTTRHs, middleTTRHs, outerTTRHs;
199 
201  std::vector<bool> innerOk( innerHits.size(), true);
202  std::vector<bool> middleOk(middleHits.size(), true);
203  std::vector<bool> outerOk( outerHits.size(), true);
204 
205  size_t sizBefore = hitTriplets.size();
207  int idx = 0;
208  for (iOuterHit = outerHits.begin(), idx = 0; iOuterHit != outerHits.end(); ++idx, ++iOuterHit){
209  outerTTRHs.push_back(ls[2].hitBuilder()->build((**iOuterHit).hit()));
210  if (checkCharge_ && !checkCharge(outerTTRHs.back()->hit())) outerOk[idx] = false;
211  }
212  for (iMiddleHit = middleHits.begin(), idx = 0; iMiddleHit != middleHits.end(); ++idx, ++iMiddleHit){
213  middleTTRHs.push_back(ls[1].hitBuilder()->build((**iMiddleHit).hit()));
214  if (checkCharge_ && !checkCharge(middleTTRHs.back()->hit())) middleOk[idx] = false;
215  }
216  for (iInnerHit = innerHits.begin(), idx = 0; iInnerHit != innerHits.end(); ++idx, ++iInnerHit){
217  innerTTRHs.push_back(ls[0].hitBuilder()->build((**iInnerHit).hit()));
218  if (checkCharge_ && !checkCharge(innerTTRHs.back()->hit())) innerOk[idx] = false;
219  }
221  checkNoisyModules(innerTTRHs, innerOk);
222  checkNoisyModules(middleTTRHs, middleOk);
223  checkNoisyModules(outerTTRHs, outerOk);
224  }
225 
226  for (iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++){
227  idx = iOuterHit - outerHits.begin();
228  TTRH & outerTTRH = outerTTRHs[idx];
229  GlobalPoint outerpos = outerTTRH->globalPosition(); // this caches by itself
230  bool outerok = outerOk[idx];
231  if (outerok < goodHitsPerSeed_ - 2) {
232  if (tripletsVerbosity_ > 2)
233  std::cout << "Skipping at first hit: " << (outerok) << " < " << (goodHitsPerSeed_ - 2) << std::endl;
234  continue;
235  }
236 
237  for (iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++){
238  idx = iMiddleHit - middleHits.begin();
239  TTRH & middleTTRH = middleTTRHs[idx];
240  GlobalPoint middlepos = middleTTRH->globalPosition(); // this caches by itself
241  bool middleok = middleOk[idx];
242  if (outerok+middleok < goodHitsPerSeed_ - 1) {
243  if (tripletsVerbosity_ > 2)
244  std::cout << "Skipping at second hit: " << (outerok+middleok) << " < " << (goodHitsPerSeed_ - 1) << std::endl;
245  continue;
246  }
247 
248  for (iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++){
249  idx = iInnerHit - innerHits.begin();
250  TTRH & innerTTRH = innerTTRHs[idx];
251  GlobalPoint innerpos = innerTTRH->globalPosition(); // this caches by itself
252  bool innerok = innerOk[idx];
253  if (outerok+middleok+innerok < goodHitsPerSeed_) {
254  if (tripletsVerbosity_ > 2)
255  std::cout << "Skipping at third hit: " << (outerok+middleok+innerok) << " < " << (goodHitsPerSeed_) << std::endl;
256  continue;
257  }
258 
259  //***top-bottom
260  if (positiveYOnly && (innerpos.y()<0 || middlepos.y()<0 || outerpos.y()<0
261  || outerpos.y() < innerpos.y()
262  ) ) continue;
263  if (negativeYOnly && (innerpos.y()>0 || middlepos.y()>0 || outerpos.y()>0
264  || outerpos.y() > innerpos.y()
265  ) ) continue;
266  //***
267 
268  if (tripletsVerbosity_ > 2) std::cout << "Trying seed with: " << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
269  if (goodTriplet(innerpos,middlepos,outerpos,minRho)) {
270  OrderedHitTriplet oht(*iInnerHit,*iMiddleHit,*iOuterHit);
271  hitTriplets.push_back(oht);
272  if ((maxTriplets_ > 0) && (hitTriplets.size() > size_t(maxTriplets_))) {
273  hitTriplets.clear(); // clear
274  //OrderedHitTriplets().swap(hitTriplets); // really clear
275  edm::LogError("TooManyTriplets") << "Found too many triplets, bailing out.\n";
276  return false;
277  }
278  if (tripletsVerbosity_ > 3) {
279  std::cout << " accepted seed #" << (hitTriplets.size()-1) << " w/: "
280  << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
281  }
282  if (tripletsVerbosity_ == 2) {
283  std::cout << " good seed #" << (hitTriplets.size()-1) << " w/: "
284  << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
285  }
286  if (tripletsVerbosity_ > 3 && (helixVerbosity_ > 0)) { // debug the momentum here too
287  pqFromHelixFit(innerpos,middlepos,outerpos,es);
288  }
289  }
290  }
291  }
292  }
293  if ((tripletsVerbosity_ > 0) && (hitTriplets.size() > sizBefore)) {
294  std::cout << " iLss = " << layerTripletNames_[iLss - lss.begin()]
295  << " (" << (iLss - lss.begin()) << "): # = "
296  << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size()
297  << ": Found " << (hitTriplets.size() - sizBefore) << " seeds [running total: " << hitTriplets.size() << "]"
298  << std::endl ;
299  }
300 
301  }
303  return true;
304 }
306  DetId detid(hit->geographicalId());
307  if (detid.det() != DetId::Tracker) return false; // should not happen
308  int subdet = detid.subdetId();
309  if (subdet < 3) { // pixel
310  return true;
311  } else {
312  if (typeid(*hit) == typeid(SiStripMatchedRecHit2D)) {
313  const SiStripMatchedRecHit2D *mhit = static_cast<const SiStripMatchedRecHit2D *>(hit);
314  if (matchedRecHitUsesAnd_) {
315  return checkCharge(mhit->monoHit(), subdet) && checkCharge(mhit->stereoHit(), subdet);
316  } else {
317  return checkCharge(mhit->monoHit(), subdet) || checkCharge(mhit->stereoHit(), subdet);
318  }
319  } else if (typeid(*hit) == typeid(SiStripRecHit2D)) {
320  return checkCharge(static_cast<const SiStripRecHit2D &>(*hit), subdet);
321  } else {
322  return true;
323  }
324  }
325 }
326 
327 // to be fixed to use OmniCluster
328 bool SimpleCosmicBONSeeder::checkCharge(const SiStripRecHit2D &hit, int subdetid) const {
329  const SiStripCluster *clust = (hit.cluster().isNonnull() ? hit.cluster().get() : hit.cluster_regional().get());
330  int charge = std::accumulate(clust->amplitudes().begin(), clust->amplitudes().end(), int(0));
331  if (tripletsVerbosity_ > 1) {
332  std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
333  << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
334  } else if ((tripletsVerbosity_ == 1) && (charge < chargeThresholds_[subdetid])) {
335  std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
336  << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
337  }
338  return charge > chargeThresholds_[subdetid];
339 }
340 
341 void SimpleCosmicBONSeeder::checkNoisyModules(const std::vector<TransientTrackingRecHit::RecHitPointer> &hits, std::vector<bool> &oks) const {
343  std::vector<TTRH>::const_iterator it = hits.begin(), start = it, end = hits.end();
344  std::vector<bool>::iterator ok = oks.begin(), okStart = ok;
345  while (start < end) {
346  DetId lastid = (*start)->geographicalId();
347  for (it = start + 1; (it < end) && ((*it)->geographicalId() == lastid); ++it) {
348  ++ok;
349  }
350  if ( (it - start) > maxHitsPerModule_[lastid.subdetId()] ) {
351  if (tripletsVerbosity_ > 0) {
352  std::cerr << "SimpleCosmicBONSeeder: Marking noisy module " << lastid.rawId() << ", it has " << (it-start) << " rechits"
353  << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
354  }
355  std::fill(okStart,ok,false);
356  } else if (tripletsVerbosity_ > 0) {
357  if ( (it - start) > std::min(4,maxHitsPerModule_[lastid.subdetId()]/4) ) {
358  std::cerr << "SimpleCosmicBONSeeder: Not marking noisy module " << lastid.rawId() << ", it has " << (it-start) << " rechits"
359  << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
360  }
361  }
362  start = it; okStart = ok;
363  }
364 }
365 
366 bool SimpleCosmicBONSeeder::goodTriplet(const GlobalPoint &inner, const GlobalPoint & middle, const GlobalPoint & outer, const double &minRho) const {
367  float dyOM = outer.y() - middle.y(), dyIM = inner.y() - middle.y();
368  if ((dyOM * dyIM > 0) && (fabs(dyOM)>10) && (fabs(dyIM)>10)) {
369  if (tripletsVerbosity_ > 2) std::cout << " fail for non coherent dy" << std::endl;
370  return false;
371  }
372  float dzOM = outer.z() - middle.z(), dzIM = inner.z() - middle.z();
373  if ((dzOM * dzIM > 0) && (fabs(dzOM)>50) && (fabs(dzIM)>50)) {
374  if (tripletsVerbosity_ > 2) std::cout << " fail for non coherent dz" << std::endl;
375  return false;
376  }
377  if (minRho > 0) {
378  FastCircle theCircle(inner,middle,outer);
379  if (theCircle.rho() < minRho) {
380  if (tripletsVerbosity_ > 2) std::cout << " fail for pt cut" << std::endl;
381  return false;
382  }
383  }
384  return true;
385 }
386 
387 std::pair<GlobalVector,int>
389  const edm::EventSetup& iSetup) const {
390  if (helixVerbosity_ > 0) {
391  std::cout << "DEBUG PZ =====" << std::endl;
392  FastHelix helix(inner,middle,outer,iSetup);
393  GlobalVector gv=helix.stateAtVertex().parameters().momentum(); // status on inner hit
394  std::cout << "FastHelix P = " << gv << "\n";
395  std::cout << "FastHelix Q = " << helix.stateAtVertex().parameters().charge() << "\n";
396  }
397 
398  // My attempt (with different approx from FastHelix)
399  // 1) fit the circle
400  FastCircle theCircle(inner,middle,outer);
401  double rho = theCircle.rho();
402 
403  // 2) Get the PT
404  GlobalVector tesla = magfield->inTesla(middle);
405  double pt = 0.01 * rho * (0.3*tesla.z());
406 
407  // 3) Get the PX,PY at OUTER hit (VERTEX)
408  double dx1 = outer.x()-theCircle.x0();
409  double dy1 = outer.y()-theCircle.y0();
410  double py = pt*dx1/rho, px = -pt*dy1/rho;
411  if(px*(middle.x() - outer.x()) + py*(middle.y() - outer.y()) < 0.) {
412  px *= -1.; py *= -1.;
413  }
414 
415  // 4) Get the PZ through pz = pT*(dz/d(R*phi)))
416  double dz = inner.z() - outer.z();
417  double sinphi = ( dx1*(inner.y()-theCircle.y0()) - dy1*(inner.x()-theCircle.x0())) / (rho * rho);
418  double dphi = std::abs(std::asin(sinphi));
419  double pz = pt * dz / (dphi * rho);
420 
421  int myq = ((theCircle.x0()*py - theCircle.y0()*px) / tesla.z()) > 0. ? +1 : -1;
422 
423  std::pair<GlobalVector,int> mypq(GlobalVector(px,py,pz),myq);
424 
425  if (helixVerbosity_ > 1) {
426  std::cout << "Gio: pt = " << pt << std::endl;
427  std::cout << "Gio: dz = " << dz << ", sinphi = " << sinphi << ", dphi = " << dphi << ", dz/drphi = " << (dz/dphi/rho) << std::endl;
428  }
429  if (helixVerbosity_ > 0) {
430  std::cout << "Gio's fit P = " << mypq.first << "\n";
431  std::cout << "Gio's fit Q = " << myq << "\n";
432  }
433 
434  return mypq;
435 }
436 
438 {
440 
441  for (size_t it=0;it<hitTriplets.size();it++){
442  const OrderedHitTriplet &trip = hitTriplets[it];
443 
444  GlobalPoint inner = tracker->idToDet((*(trip.inner())).geographicalId())->surface().
445  toGlobal((*(trip.inner())).localPosition());
446 
447  GlobalPoint middle = tracker->idToDet((*(trip.middle())).geographicalId())->surface().
448  toGlobal((*(trip.middle())).localPosition());
449 
450  GlobalPoint outer = tracker->idToDet((*(trip.outer())).geographicalId())->surface().
451  toGlobal((*(trip.outer())).localPosition());
452 
453  if (seedVerbosity_ > 1)
454  std::cout << "Processing triplet " << it << ": " << inner << " + " << middle << " + " << outer << std::endl;
455 
456  if ( (outer.y()-inner.y())*outer.y() < 0 ) {
457  std::swap(inner,outer);
458  std::swap(const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.inner()),
459  const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.outer()) );
460 
461 // std::swap(const_cast<ctfseeding::SeedingHit &>(trip.inner()),
462 // const_cast<ctfseeding::SeedingHit &>(trip.outer()) );
463  if (seedVerbosity_ > 1) {
464  std::cout << "The seed was going away from CMS! swapped in <-> out" << std::endl;
465  std::cout << "Processing swapped triplet " << it << ": " << inner << " + " << middle << " + " << outer << std::endl;
466  }
467  }
468 
469  // First use FastHelix out of the box
470  std::pair<GlobalVector,int> pq = pqFromHelixFit(inner,middle,outer,iSetup);
471  GlobalVector gv = pq.first;
472  float ch = pq.second;
473  float Mom = sqrt( gv.x()*gv.x() + gv.y()*gv.y() + gv.z()*gv.z() );
474 
475  if(Mom > 10000 || isnan(Mom)) {
476  if (seedVerbosity_ > 1)
477  std::cout << "Processing triplet " << it << ": fail for momentum." << std::endl;
478  continue;
479  }
480 
481  if (gv.perp() < region_.ptMin()) {
482  if (seedVerbosity_ > 1)
483  std::cout << "Processing triplet " << it << ": fail for pt = " << gv.perp() << " < ptMin = " << region_.ptMin() << std::endl;
484  continue;
485  }
486 
487  const Propagator * propagator = 0;
488  if((outer.y()-inner.y())>0){
489  if (seedVerbosity_ > 1)
490  std::cout << "Processing triplet " << it << ": downgoing." << std::endl;
491  propagator = thePropagatorAl;
492  } else {
493  gv = -1*gv; ch = -1.*ch;
494  propagator = thePropagatorOp;
495  if (seedVerbosity_ > 1)
496  std::cout << "Processing triplet " << it << ": upgoing." << std::endl;
497  }
498 
499  if (seedVerbosity_ > 1) {
500  if (( gv.z() * (outer.z()-inner.z()) > 0 ) && ( fabs(outer.z()-inner.z()) > 5) && (fabs(gv.z()) > .01)) {
501  std::cout << "ORRORE: outer.z()-inner.z() = " << (outer.z()-inner.z()) << ", gv.z() = " << gv.z() << std::endl;
502  }
503  }
504 
505  GlobalTrajectoryParameters Gtp(outer,
506  gv,int(ch),
507  &(*magfield));
508  FreeTrajectoryState CosmicSeed(Gtp,
510  CosmicSeed.rescaleError(100);
511  if (seedVerbosity_ > 2) {
512  std::cout << "Processing triplet " << it << ". start from " << std::endl;
513  std::cout << " X = " << outer << ", P = " << gv << std::endl;
514  std::cout << " Cartesian error (X,P) = \n" << CosmicSeed.cartesianError().matrix() << std::endl;
515  }
516 
518  OrderedHitTriplet seedHits(trip.outer(),trip.middle(),trip.inner());
519  TSOS propagated, updated;
520  bool fail = false;
521  for (size_t ih = 0; ih < 3; ++ih) {
522  if ((ih == 2) && seedOnMiddle_) {
523  if (seedVerbosity_ > 2)
524  std::cout << "Stopping at middle hit, as requested." << std::endl;
525  break;
526  }
527  if (seedVerbosity_ > 2)
528  std::cout << "Processing triplet " << it << ", hit " << ih << "." << std::endl;
529  if (ih == 0) {
530  propagated = propagator->propagate(CosmicSeed, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
531  } else {
532  propagated = propagator->propagate(updated, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
533  }
534  if (!propagated.isValid()) {
535  if (seedVerbosity_ > 1)
536  std::cout << "Processing triplet " << it << ", hit " << ih << ": failed propagation." << std::endl;
537  fail = true; break;
538  } else {
539  if (seedVerbosity_ > 2)
540  std::cout << "Processing triplet " << it << ", hit " << ih << ": propagated state = " << propagated;
541  }
542  const TransientTrackingRecHit::ConstRecHitPointer & tthp = seedHits[ih];
543  TransientTrackingRecHit::RecHitPointer newtth = tthp->clone(propagated);
544  hits.push_back(newtth->hit()->clone());
545  updated = theUpdator->update(propagated, *newtth);
546  if (!updated.isValid()) {
547  if (seedVerbosity_ > 1)
548  std::cout << "Processing triplet " << it << ", hit " << ih << ": failed update." << std::endl;
549  fail = true; break;
550  } else {
551  if (seedVerbosity_ > 2)
552  std::cout << "Processing triplet " << it << ", hit " << ih << ": updated state = " << updated;
553  }
554  }
555  if (!fail && updated.isValid() && (updated.globalMomentum().perp() < region_.ptMin())) {
556  if (seedVerbosity_ > 1)
557  std::cout << "Processing triplet " << it <<
558  ": failed for final pt " << updated.globalMomentum().perp() << " < " << region_.ptMin() << std::endl;
559  fail = true;
560  }
561  if (!fail && updated.isValid() && (updated.globalMomentum().mag() < pMin_)) {
562  if (seedVerbosity_ > 1)
563  std::cout << "Processing triplet " << it <<
564  ": failed for final p " << updated.globalMomentum().perp() << " < " << pMin_ << std::endl;
565  fail = true;
566  }
567  if (!fail) {
568  if (rescaleError_ != 1.0) {
569  if (seedVerbosity_ > 2) {
570  std::cout << "Processing triplet " << it << ", rescale error by " << rescaleError_ << ": state BEFORE rescaling " << updated;
571  std::cout << " Cartesian error (X,P) before rescaling= \n" << updated.cartesianError().matrix() << std::endl;
572  }
573  updated.rescaleError(rescaleError_);
574  }
575  if (seedVerbosity_ > 0) {
576  std::cout << "Processed triplet " << it << ": success (saved as #"<<output.size()<<") : "
577  << inner << " + " << middle << " + " << outer << std::endl;
578  std::cout << " pt = " << updated.globalMomentum().perp() <<
579  " eta = " << updated.globalMomentum().eta() <<
580  " phi = " << updated.globalMomentum().phi() <<
581  " ch = " << updated.charge() << std::endl;
582  if (seedVerbosity_ > 1) {
583  std::cout << " State:" << updated;
584  } else {
585  std::cout << " X = " << updated.globalPosition() << ", P = " << updated.globalMomentum() << std::endl;
586  }
587  std::cout << " Cartesian error (X,P) = \n" << updated.cartesianError().matrix() << std::endl;
588  }
589 
591  (*(seedOnMiddle_ ? trip.middle() : trip.inner())).geographicalId().rawId());
592  output.push_back(TrajectorySeed(PTraj,hits,
593  ( (outer.y()-inner.y()>0) ? alongMomentum : oppositeToMomentum) ));
594  if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
595  output.clear();
596  edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
597  return false;
598  }
599  }
600  }
601  return true;
602 }
603 
605  delete thePropagatorAl;
606  delete thePropagatorOp;
607  delete theUpdator;
608 }
609 
610 
T getParameter(std::string const &) const
GlobalTrackingRegion region_
FTS stateAtVertex() const
Definition: FastHelix.cc:7
const OuterRecHit & outer() const
bool triplets(const edm::Event &e, const edm::EventSetup &c)
void init(const edm::EventSetup &c)
string fill
Definition: lumiContext.py:319
CartesianTrajectoryError cartesianError() const
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:187
T perp() const
Definition: PV3DBase.h:71
const GlobalTrajectoryParameters & parameters() const
double x0() const
Definition: FastCircle.h:50
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
OrderedHitTriplets hitTriplets
Definition: DDAxes.h:10
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
T y() const
Definition: PV3DBase.h:62
const InnerRecHit & inner() const
virtual void produce(edm::Event &e, const edm::EventSetup &c)
const MiddleRecHit & middle() const
#define abs(x)
Definition: mlp_lapack.h:159
#define min(a, b)
Definition: mlp_lapack.h:161
uint16_t firstStrip() const
void checkNoisyModules(const std::vector< TransientTrackingRecHit::RecHitPointer > &hits, std::vector< bool > &oks) const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
bool checkCharge(const TrackingRecHit *hit) const
double charge(const std::vector< uint8_t > &Ampls)
double rho() const
Definition: FastCircle.h:54
SimpleCosmicBONSeeder(const edm::ParameterSet &conf)
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
edm::ESHandle< TrackerGeometry > tracker
std::pair< GlobalVector, int > pqFromHelixFit(const GlobalPoint &inner, const GlobalPoint &middle, const GlobalPoint &outer, const edm::EventSetup &iSetup) const
virtual TrackingRegion::Hits hits(const edm::Event &ev, const edm::EventSetup &es, const ctfseeding::SeedingLayer *layer) const
get hits from layer compatible with region constraints
void push_back(D *&d)
Definition: OwnVector.h:273
SeedingLayerSetsBuilder theLsb
std::vector< TrajectorySeed > TrajectorySeedCollection
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
bool isnan(float x)
Definition: math.h:13
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
virtual float ptMin() const
minimal pt of interest
T sqrt(T t)
Definition: SSEVec.h:46
PropagatorWithMaterial * thePropagatorOp
T z() const
Definition: PV3DBase.h:63
edm::ESHandle< MagneticField > magfield
size_t tooManyClusters(const edm::Event &e) const
edm::ESHandle< TransientTrackingRecHitBuilder > TTTRHBuilder
#define end
Definition: vmac.h:38
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
tuple conf
Definition: dbtoconf.py:185
PropagatorWithMaterial * thePropagatorAl
const AlgebraicSymMatrix66 & matrix() const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:12
virtual TrackingRecHit * clone() const =0
Definition: DetId.h:20
TrajectoryStateOnSurface TSOS
Definition: TestHits.cc:23
void rescaleError(double factor)
const T & get() const
Definition: EventSetup.h:55
bool seeds(TrajectorySeedCollection &output, const edm::EventSetup &iSetup)
bool goodTriplet(const GlobalPoint &inner, const GlobalPoint &middle, const GlobalPoint &outer, const double &minRho) const
bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const
double ptmin
Definition: HydjetWrapper.h:86
double y0() const
Definition: FastCircle.h:52
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
Definition: KFUpdator.cc:10
ctfseeding::SeedingLayerSets layers(const edm::EventSetup &es) const
std::vector< std::string > layerTripletNames_
TransientTrackingRecHit::ConstRecHitPointer SeedingHit
std::vector< int32_t > chargeThresholds_
tuple cout
Definition: gather_cfg.py:121
std::vector< int32_t > maxHitsPerModule_
DetId geographicalId() const
virtual unsigned int size() const
T x() const
Definition: PV3DBase.h:61
const std::vector< uint8_t > & amplitudes() const
std::vector< SeedingLayer > SeedingLayers
Global3DVector GlobalVector
Definition: GlobalVector.h:10
std::vector< std::vector< SeedingLayer > > SeedingLayerSets