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