CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
OutInConversionSeedFinder.cc
Go to the documentation of this file.
4 //
6 // Field
8 //
10 // Geometry
13 //
18 //
21 
22 //
23 #include "CLHEP/Units/GlobalPhysicalConstants.h"
24 #include "CLHEP/Geometry/Point3D.h"
25 #include "CLHEP/Geometry/Vector3D.h"
26  #include "CLHEP/Geometry/Transform3D.h"
27 #include <cfloat>
28 
29 
31 {
32 
33  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR " << "\n";
34 
35  //the2ndHitdphi_ = 0.01;
36  the2ndHitdphi_ = 0.03;
37  the2ndHitdzConst_ = 5.;
38  the2ndHitdznSigma_ = 2.;
40 
41 
42 }
43 
44 
45 
46 
48  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR " << "\n";
49 
50 }
51 
52 
53 
55 
56  // std::cout << " OutInConversionSeedFinder::makeSeeds() " << "\n";
57  theSeeds_.clear();
58 
59 
60  // debug
61  // std::cout << " OutInConversionSeedFinder::makeSeeds() SC position " << theSCPosition_.x() << " " << theSCPosition_.y() << " " << theSCPosition_.z() << "\n";
62  // std::cout << " SC eta " << theSCPosition_.eta() << " SC phi " << theSCPosition_.phi() << "\n";
63  // std::cout << " OutInConversionSeedFinder::makeSeeds() SC energy " << theSCenergy_ << "\n";
64  //
65 
66  findLayers();
67 
68 
69  // std::cout << " Check Calo cluster collection size " << allBC->size() << "\n";
70 
71  float theSCPhi=theSCPosition_.phi();
72  float theSCEta=theSCPosition_.eta();
73 
74 
75 
76  // Loop over the Calo Clusters in the event looking for seeds
77  reco::CaloClusterCollection::const_iterator bcItr;
78  LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() All BC in the event " << "\n";
79  for (unsigned i = 0; i < allBC->size(); ++i ) {
80 
81  //for(bcItr = allBC.begin(); bcItr != allBC.end(); bcItr++) {
82  nSeedsPerBC_=0;
83 
84  theBCEnergy_=allBC->ptrAt(i)->energy();
85  theBCPosition_ = GlobalPoint(allBC->ptrAt(i)->position().x(), allBC->ptrAt(i)->position().y(), allBC->ptrAt(i)->position().z() ) ;
86  float theBcEta= theBCPosition_.eta();
87  float theBcPhi= theBCPosition_.phi();
88  // float dPhi= theBcPhi-theSCPhi;
89 
90  // std::cout << " OutInConversionSeedFinder::makeSeeds() BC eta " << theBcEta << " phi " << theBcPhi << " BC energy " << theBCEnergy_ << " dPhi " << fabs(theBcPhi-theSCPhi) << " dEta " << fabs(theBcEta-theSCEta) << "\n";
91 
92  if ( theBCEnergy_ < 1.5 ) continue;
93 
94  LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() Passing the >=1.5 GeV cut BC eta " << theBcEta << " phi " << theBcPhi << " BC energy " << theBCEnergy_ << "\n";
95 
96  if ( fabs(theBcEta-theSCEta) < 0.015 && fabs(theBcPhi-theSCPhi) < 0.3 ) {
97  LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() in et and phi range passed to the analysis " << "\n";
98  fillClusterSeeds( allBC->ptrAt(i) );
99  }
100 
101 
102  }
103 
104 
105  // std::cout << "Built vector of seeds of size " << theSeeds_.size() << "\n" ;
106 
108  /*
109  int nSeed=0;
110  for ( std::vector<TrajectorySeed>::const_iterator iSeed= theSeeds_.begin(); iSeed != theSeeds_.end(); ++iSeed) {
111  nSeed++;
112  PTrajectoryStateOnDet ptsod=iSeed->startingState();
113  LogDebug("OutInConversionSeedFinder") << nSeed << ") Direction " << iSeed->direction() << " Num of hits " << iSeed->nHits() << " starting state position " << ptsod.parameters().position() << " R " << ptsod.parameters().position().perp() << " phi " << ptsod.parameters().position().phi() << " eta " << ptsod.parameters().position().eta() << "\n" ;
114 
115 
116  DetId tmpId = DetId( iSeed->startingState().detId());
117  const GeomDet* tmpDet = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
118  GlobalVector gv = tmpDet->surface().toGlobal( iSeed->startingState().parameters().momentum() );
119 
120  LogDebug("OutInConversionSeedFinder") << "seed perp,phi,eta : "
121  << gv.perp() << " , "
122  << gv.phi() << " , "
123  << gv.eta() << "\n" ; ;
124 
125 
126 
127 
128  TrajectorySeed::range hitRange = iSeed->recHits();
129  for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
130 
131  if ( ihit->isValid() ) {
132 
133  LogDebug("OutInConversionSeedFinder") << " Valid hit global position " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()) << " R " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).perp() << " phi " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).phi() << " eta " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).eta() << "\n" ;
134 
135  }
136  }
137  }
138 
139  */
140 
141 
142 
143 }
144 
145 
146 
148 
149  theSeeds_.clear();
150 
151  findLayers();
152 
153  float theSCPhi=theSCPosition_.phi();
154  float theSCEta=theSCPosition_.eta();
155 
156  nSeedsPerBC_=0;
157 
158  theBCEnergy_=aBC->energy();
159  theBCPosition_ = GlobalPoint(aBC->position().x(), aBC->position().y(), aBC->position().z() ) ;
160  float theBcEta= theBCPosition_.eta();
161  float theBcPhi= theBCPosition_.phi();
162  // float dPhi= theBcPhi-theSCPhi;
163 
164  if ( theBCEnergy_ < 1.5 ) return;
165 
166  if ( fabs(theBcEta-theSCEta) < 0.015 && fabs(theBcPhi-theSCPhi) < 0.25 ) {
167  fillClusterSeeds( aBC);
168  }
169 
170 }
171 
172 
173 
174 
175 
177 
178 
179  theFirstMeasurements_.clear();
180  FreeTrajectoryState fts;
181 
183  if ( makeTrackState(-1).second ) {
184  fts = makeTrackState(-1).first;
185  startSeed(fts);
186  }
188 
189  if ( makeTrackState(1).second ) {
190  fts = makeTrackState(1).first;
191  startSeed(fts);
192  }
193 
194 }
195 
196 
197 
198 std::pair<FreeTrajectoryState,bool> OutInConversionSeedFinder::makeTrackState(int charge) const {
199 
200  std::pair<FreeTrajectoryState,bool> result;
201  result.second=false;
202 
203 
204  //std::cout << " OutInConversionSeedFinder:makeTrackState " << "\n";
205 
206 
207  // Old GlobalPoint gpOrigine(theBCPosition_.x()*0.3, theBCPosition_.y()*0.3, theBCPosition_.z()*0.3) ;
208  // GlobalPoint gpOrigine(0.,0.,0.);
209 
211  GlobalVector gvBcRadius = theBCPosition_ - gpOrigine ;
212  HepGeom::Point3D<double> radiusBc(gvBcRadius.x(),gvBcRadius.y(),gvBcRadius.z()) ;
213  HepGeom::Point3D<double> momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_ ;
214 
215  // compute momentum direction at calo
216  double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp() ;
217  curvature /= 100. ; // in cm-1 !!
218 
219  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x() << " " << gpOrigine.y() << " " << gpOrigine.z() << " momentumWithoutCurvature " << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
220 
221  // define rotation angle
222  float R = theBCPosition_.perp();
223  float r = gpOrigine.perp();
224  float rho = 1./curvature;
225  // from the formula for the intersection of two circles
226  // turns out to be about 2/3 of the deflection of the old formula
227  float d = sqrt(r*r+rho*rho);
228  float u = rho + rho/d/d*(R*R-rho*rho) - r/d/d*sqrt((R*R-r*r+2*rho*R)*(R*R-r*r+2*rho*R));
229  //float u = rho + rho/d/d*(R*R-rho*rho) ;
230  if ( u <=R ) result.second=true;
231 
232  double sinAlpha = 0.5*u/R;
233  if ( sinAlpha>(1.-10*DBL_EPSILON) ) sinAlpha = 1.-10*DBL_EPSILON;
234  else if ( sinAlpha<-(1.-10*DBL_EPSILON) ) sinAlpha = -(1.-10*DBL_EPSILON);
235 
236  double newdphi = charge * asin( sinAlpha) ;
237 
238  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R << " u/R " << u/R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
239 
240  HepGeom::Transform3D rotation = HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double> (0., 0. ,1.));
241 
242 
243  HepGeom::Point3D<double> momentumInTracker = momentumWithoutCurvature.transform(rotation) ;
244  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState R " << R << " r " << r << " rho " << rho << " d " << d << " u " << u << " newdphi " << newdphi << " momentumInTracker " << momentumInTracker << "\n";
245 
246  HepGeom::Point3D<double> hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z()) ;
247 
248  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint " << hepStartingPoint << "\n";
249 
250  hepStartingPoint.transform(rotation);
251 
252  GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
253 
254  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint << " calo position " << theBCPosition_ << "\n";
255  GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
256  GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_) );
257  // error matrix
259  m(0,0) = 0.1; m(1,1) = 0.1 ; m(2,2) = 0.1 ;
260  m(3,3) = 0.1 ; m(4,4) = 0.1;
261 
262  //std::cout << "OutInConversionSeedFinder::makeTrackState " << FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) << "\n";
263 
264  result.first= FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) ;
265  return result;
266 
267 }
268 
269 
271 
272 
273  // std::cout << "OutInConversionSeedFinder::startSeed layer list " << this->layerList().size() << "\n";
274  //std::cout << "OutInConversionSeedFinder::startSeed fts " << fts << "\n";
275 
276  std::vector<const DetLayer*> myLayers=layerList();
277  if ( myLayers.size() > 3 ) {
278 
279  for(unsigned int ilayer = myLayers.size()-1; ilayer >= myLayers.size()-2; --ilayer) {
280  const DetLayer * layer = myLayers[ilayer];
281 
282 
283  // allow the z of the hit to be within a straight line from a vertex
284  // of +-15 cm to the cluster
285  // float dphi = 0.015;
286  float dphi = 0.030;
287 
288  MeasurementEstimator * newEstimator = makeEstimator(layer, dphi);
289 
290 
291  //std::cout << "OutInSeedFinder::startSeed propagationDirection " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";
292 
293  TSOS tsos(fts, layer->surface() );
294 
295  LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed after TSOS tsos(fts, layer->surface() ) " << "\n";
296 
297  LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
298  theFirstMeasurements_ = theLayerMeasurements_.measurements( *layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
299 
300  //std::cout << "OutInSeedFinder::startSeed after theFirstMeasurements_ " << theFirstMeasurements_.size() << "\n";
301 
302  if(theFirstMeasurements_.size() > 1) // always a dummy returned, too
303  LogDebug("OutInConversionSeedFinder") << " Found " << theFirstMeasurements_.size()-1 << " 1st hits in seed" << "\n";
304 
305  delete newEstimator;
306 
307  LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
308 
309  for(unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
311  if(m1.recHit()->isValid()) {
312 
313  // update the fts to start from this point. much better than starting from
314  // extrapolated point along the line
315  GlobalPoint hitPoint = m1.recHit()->globalPosition();
316  LogDebug("OutInConversionSeedFinder") << " Valid hit at R " << m1.recHit()->globalPosition().perp() << " Z " << m1.recHit()->globalPosition().z() << " eta " << m1.recHit()->globalPosition().eta() << " phi " << m1.recHit()->globalPosition().phi() << " xyz " << m1.recHit()->globalPosition() << "\n";
317 
318 
319  FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
320  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed newfts " << newfts << "\n";
321  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed propagationDirection after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
322  // std::cout << "OutInConversionSeedFinder::startSeed propagationDirection after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
323 
324 
325  completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-1);
326  // skip a layer, if you haven't already skipped the first layer
327  if(ilayer == myLayers.size()-1) {
328  completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-2);
329  }
330  }
331  }
332 
333  } // loop over layers
334  }
335 
336 
337 
338 
339 }
340 
341 
342 
344 
345  //std::cout << "OutInConversionSeedFinder::makeEstimator " << "\n";
346 
347  MeasurementEstimator * newEstimator=0;
348 
349  if (layer->location() == GeomDetEnumerators::barrel ) {
350 
351  const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
352  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Barrel r = " << barrelLayer->specificSurface().radius() << " " << "\n";
353  float r = barrelLayer->specificSurface().radius();
354  float zrange = 15.* (1.-r/theBCPosition_.perp());
355  newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
356  }
357 
358 
359 
360  if (layer->location() == GeomDetEnumerators::endcap ) {
361 
362  const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
363  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius() << " R " << forwardLayer->specificSurface().outerRadius() << " Z " << forwardLayer->specificSurface().position().z() << "\n";
364 
365  float zc = fabs(theBCPosition_.z());
366  float z = fabs(forwardLayer->surface().position().z());
367 
368  float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc*zc - 15.*zc);
369  newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
370  }
371 
372 
373 
374 
375  return newEstimator;
376 }
377 
378 
379 
380 
382  FreeTrajectoryState & fts,
383  const Propagator* propagator, int ilayer) const {
384 
385  //std::cout << "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
386 
387  MeasurementEstimator * newEstimator=0;
388  const DetLayer * layer = theLayerList_[ilayer];
389 
390 
391  if ( layer->location() == GeomDetEnumerators::barrel ) {
392  // z error for 2nd hit is 2 sigma quadded with 5 cm
393  LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
394  float dz = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1.recHit()->globalPositionError().czz()
396  newEstimator =
397  new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
398  }
399  else {
400  LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
401  // z error for 2nd hit is 2sigma quadded with 5 cm
402  //float m1dr = m1.recHit().globalPositionError().rerr(m1.recHit().globalPosition());
403  float m1dr = sqrt(m1.recHit()->localPositionError().yy());
404  float dr = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1dr*m1dr
406  // LogDebug("OutInConversionSeedFinder") << "second hit forward dr " << dr << " this hit " << m1dr << endl;
407  newEstimator =
408  new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
409  }
410 
411  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
412 
413  // Get the measurements consistent with the FTS and the Estimator
414  TSOS tsos(fts, layer->surface() );
415  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection " << int(propagator->propagationDirection() ) << "\n";
416  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
417 
418  LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
419  std::vector<TrajectoryMeasurement> measurements = theLayerMeasurements_.measurements( *layer, tsos, *propagator, *newEstimator);
420  //std::cout << "OutInConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
421  delete newEstimator;
422 
423  for(unsigned int i = 0; i < measurements.size(); ++i) {
424  if( measurements[i].recHit()->isValid() ) {
425  createSeed(m1, measurements[i]);
426  }
427  }
428 
429 
430 
431  //LogDebug("OutInConversionSeedFinder") << "COMPLETED " << theSeeds_.size() << " SEEDS " << "\n";
432 }
433 
434 
435 
437  const TrajectoryMeasurement & m2) const {
438 
439  //std::cout << "OutInConversionSeedFinder::createSeed from hit1 " << m1.recHit()->globalPosition() << " r1 " << m1.recHit()->globalPosition().perp() << " and hit2 " << m2.recHit()->globalPosition() << " r2 " << m2.recHit()->globalPosition().perp() << "\n";
440 
441 
442  FreeTrajectoryState fts = createSeedFTS(m1, m2);
443 
444 
445  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
446  LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts <<"\n";
447 
448 
449 
450  //std::cout << "OutInConversionSeedFinder::createSeed propagation dir " << int( thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
451  TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts, m1.recHit()->det()->surface());
452 
453  // LogDebug("OutInConversionSeedFinder") << "hit surface " << h1.det().surface().position() << endl;
454  // LogDebug("OutInConversionSeedFinder") << "prop to " << typeid(h1.det().surface()).name() << endl;
455  // LogDebug("OutInConversionSeedFinder") << "prop to first hit " << state1 << endl;
456  // LogDebug("OutInConversionSeedFinder") << "update to " << h1.globalPosition() << endl;
457 
458  if ( state1.isValid() ) {
459  TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit() );
460 
461  if ( updatedState1.isValid() ) {
462  TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(*updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
463 
464  if ( state2.isValid() ) {
465 
466  TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit() );
467  TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit() , m1.estimate(), m1.layer());
468  TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit() , m2.estimate(), m2.layer());
469 
471  myHits.push_back(meas1.recHit()->hit()->clone());
472  myHits.push_back(meas2.recHit()->hit()->clone());
473 
474  if ( nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_ ) return;
475 
476  TrajectoryStateTransform tsTransform;
477  PTrajectoryStateOnDet* ptsod= tsTransform.persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId() );
478 
479  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed from state " << state2.globalPosition() << "\n";
480  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed ptsod " << ptsod->parameters().position() << " R " << ptsod->parameters().position().perp() << " phi " << ptsod->parameters().position().phi() << " eta " << ptsod->parameters().position().eta() << "\n" ;
481 
482 
483 
484  theSeeds_.push_back(TrajectorySeed( *ptsod, myHits, oppositeToMomentum ));
485  nSeedsPerBC_++;
486 
487  delete ptsod;
488 
489 
490  }
491  }
492  }
493 
494 
495 
496 
497 
498 }
499 
500 
501 
502 
503 
505  const TrajectoryMeasurement & m2) const {
506 
507  //std::cout << "OutInConversionSeedFinder::createSeedFTS " << "\n";
508 
509  GlobalPoint xmeas = fixPointRadius(m1);
510  GlobalPoint xvert = fixPointRadius(m2);
511 
512 
513  float pt = theSCenergy_ * sin(theSCPosition_.theta());
514  float pz = theSCenergy_ * cos(theSCPosition_.theta());
515 
516 
517 
518  // doesn't work at all for endcap, where r is badly measured
519  //float dphidr = (p1.phi()-p2.phi())/(p1.perp()-p2.perp());
520  //int charge = (dphidr > 0.) ? -1 : 1;
521  int charge = m2.predictedState().charge();
522 
523  double BInTesla = theMF_->inTesla(xmeas).z();
524  GlobalVector xdiff = xmeas -xvert;
525 
526  double phi= xdiff.phi();
527  double pxOld = pt*cos(phi);
528  double pyOld = pt*sin(phi);
529  double RadCurv = 100*pt/(BInTesla*0.29979);
530  double alpha = asin(0.5*xdiff.perp()/RadCurv);
531  float ca = cos(charge*alpha);
532  float sa = sin(charge*alpha);
533  double pxNew = ca*pxOld + sa*pyOld;
534  double pyNew = -sa*pxOld + ca*pyOld;
535  GlobalVector pNew(pxNew, pyNew, pz);
536 
537  GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_) );
538 
540  m(0,0) = 0.05; m(1,1) = 0.02 ; m(2,2) = 0.007 ;
541  m(3,3) = 10. ; m(4,4) = 10. ;
543 
544 
545 }
546 
547 
548 
549 
551  GlobalPoint p1 = m1.recHit()->globalPosition();
552  GlobalPoint p2;
553  if(m1.layer()->location() == GeomDetEnumerators::barrel) {
554  p2 = p1;
555  } else {
556  float z = p1.z();
557  float phi = p1.phi();
558  float theta = theSCPosition_.theta();
559  float r = p1.z() * tan(theta);
560  p2 = GlobalPoint(r*cos(phi), r*sin(phi), z);
561  // LogDebug("OutInConversionSeedFinder") << "fixing point radius " << p2 << " from " << p1 << endl;
562  }
563  return p2;
564 }
565 
566 
567 
568 
569 
#define LogDebug(id)
const Propagator * thePropagatorAlongMomentum_
std::vector< TrajectoryMeasurement > measurements(const DetLayer &layer, const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
int i
Definition: DBlmapReader.cc:9
float alpha
Definition: AMPTWrapper.h:95
void fillClusterSeeds(const reco::CaloClusterPtr &bc) const
const MeasurementTracker * getMeasurementTracker() const
std::pair< FreeTrajectoryState, bool > makeTrackState(int charge) const
T perp() const
Definition: PV3DBase.h:66
virtual PropagationDirection propagationDirection() const
Definition: Propagator.h:143
virtual Location location() const =0
Which part of the detector (barrel, endcap)
LocalPoint position() const
Local x and y position coordinates.
void createSeed(const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2) const
FreeTrajectoryState * freeTrajectoryState(bool withErrors=true) const
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
virtual const BoundSurface & surface() const
The surface of the GeometricSearchDet.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Definition: DDAxes.h:10
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
MeasurementEstimator * makeEstimator(const DetLayer *, float dphi) const
GlobalPoint globalPosition() const
Geom::Theta< T > theta() const
OutInConversionSeedFinder(const edm::ParameterSet &config)
ConstRecHitPointer recHit() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
TrajectorySeedCollection theSeeds_
GlobalPoint fixPointRadius(const TrajectoryMeasurement &) const
const Propagator * thePropagatorOppositeToMomentum_
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
U second(std::pair< T, U > const &p)
void push_back(D *&d)
Definition: OwnVector.h:288
T curvature(T InversePt, const edm::EventSetup &iSetup)
Definition: DDAxes.h:10
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:55
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
tuple result
Definition: query.py:137
virtual const BoundDisk & specificSurface() const
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
const DetLayer * layer() const
std::vector< const DetLayer * > theLayerList_
virtual void makeSeeds(const edm::Handle< edm::View< reco::CaloCluster > > &allBc) const
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) const
TrajectoryStateOnSurface predictedState() const
double p2[4]
Definition: TauolaWrapper.h:90
std::vector< const DetLayer * > layerList() const
tuple conf
Definition: dbtoconf.py:185
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:9
edm::ESHandle< MagneticField > theMF_
virtual const BoundCylinder & specificSurface() const
Extension of the interface.
FreeTrajectoryState trackStateFromClusters(int aCharge, const GlobalPoint &gpOrigine, PropagationDirection dir, float scaleFactor) const
T eta() const
Definition: PV3DBase.h:70
double p1[4]
Definition: TauolaWrapper.h:89
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
Definition: KFUpdator.cc:10
FreeTrajectoryState createSeedFTS(const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2) const
float outerRadius() const
The outer radius of the disk.
Definition: BoundDisk.cc:10
float innerRadius() const
The inner radius of the disk.
Definition: BoundDisk.cc:6
const Point & position() const
position
Definition: BeamSpot.h:63
void completeSeed(const TrajectoryMeasurement &m1, FreeTrajectoryState &fts, const Propagator *, int layer) const
const PositionType & position() const
void startSeed(const FreeTrajectoryState &) const
std::vector< TrajectoryMeasurement > theFirstMeasurements_
const LocalTrajectoryParameters & parameters() const
Definition: DDAxes.h:10