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