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  FreeTrajectoryState fts;
174 
176  if (makeTrackState(-1).second) {
177  fts = makeTrackState(-1).first;
178  startSeed(fts);
179  }
181 
182  if (makeTrackState(1).second) {
183  fts = makeTrackState(1).first;
184  startSeed(fts);
185  }
186  theFirstMeasurements_.clear();
187 }
188 
189 std::pair<FreeTrajectoryState, bool> OutInConversionSeedFinder::makeTrackState(int charge) const {
190  std::pair<FreeTrajectoryState, bool> result;
191  result.second = false;
192 
193  //std::cout << " OutInConversionSeedFinder:makeTrackState " << "\n";
194 
195  // Old GlobalPoint gpOrigine(theBCPosition_.x()*0.3, theBCPosition_.y()*0.3, theBCPosition_.z()*0.3) ;
196  // GlobalPoint gpOrigine(0.,0.,0.);
197 
198  GlobalPoint gpOrigine(theBeamSpot_.position().x(), theBeamSpot_.position().y(), theBeamSpot_.position().z());
199  GlobalVector gvBcRadius = theBCPosition_ - gpOrigine;
200  HepGeom::Point3D<double> radiusBc(gvBcRadius.x(), gvBcRadius.y(), gvBcRadius.z());
201  HepGeom::Point3D<double> momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_;
202 
203  // compute momentum direction at calo
204  double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp();
205  curvature /= 100.; // in cm-1 !!
206 
207  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x()
208  << " " << gpOrigine.y() << " " << gpOrigine.z() << " momentumWithoutCurvature "
209  << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
210 
211  // define rotation angle
212  float R = theBCPosition_.perp();
213  float r = gpOrigine.perp();
214  float rho = 1. / curvature;
215  // from the formula for the intersection of two circles
216  // turns out to be about 2/3 of the deflection of the old formula
217  float d = sqrt(r * r + rho * rho);
218  float u = rho + rho / d / d * (R * R - rho * rho) -
219  r / d / d * sqrt((R * R - r * r + 2 * rho * R) * (R * R - r * r + 2 * rho * R));
220  //float u = rho + rho/d/d*(R*R-rho*rho) ;
221  if (u <= R)
222  result.second = true;
223 
224  double sinAlpha = 0.5 * u / R;
225  if (sinAlpha > (1. - 10 * DBL_EPSILON))
226  sinAlpha = 1. - 10 * DBL_EPSILON;
227  else if (sinAlpha < -(1. - 10 * DBL_EPSILON))
228  sinAlpha = -(1. - 10 * DBL_EPSILON);
229 
230  double newdphi = charge * asin(sinAlpha);
231 
232  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R
233  << " u/R " << u / R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
234 
235  HepGeom::Transform3D rotation = HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double>(0., 0., 1.));
236 
237  HepGeom::Point3D<double> momentumInTracker = momentumWithoutCurvature.transform(rotation);
238  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState R " << R << " r " << r << " rho "
239  << rho << " d " << d << " u " << u << " newdphi " << newdphi
240  << " momentumInTracker " << momentumInTracker << "\n";
241 
242  HepGeom::Point3D<double> hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z());
243 
244  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint "
245  << hepStartingPoint << "\n";
246 
247  hepStartingPoint.transform(rotation);
248 
249  GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
250 
251  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint
252  << " calo position " << theBCPosition_ << "\n";
253  GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
254  GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_));
255  // error matrix
256  AlgebraicSymMatrix55 m = AlgebraicMatrixID();
257  m(0, 0) = 0.1;
258  m(1, 1) = 0.1;
259  m(2, 2) = 0.1;
260  m(3, 3) = 0.1;
261  m(4, 4) = 0.1;
262 
263  // std::cout << "OutInConversionSeedFinder::makeTrackState " << FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) << "\n";
264 
265  result.first = FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m));
266  return result;
267 }
268 
269 void OutInConversionSeedFinder::startSeed(const FreeTrajectoryState& fts) {
270  // std::cout << "OutInConversionSeedFinder::startSeed layer list " << this->layerList().size() << "\n";
271  //std::cout << "OutInConversionSeedFinder::startSeed fts " << fts << "\n";
272 
273  std::vector<const DetLayer*> myLayers = layerList();
274  if (myLayers.size() > 3) {
275  for (unsigned int ilayer = myLayers.size() - 1; ilayer >= myLayers.size() - 2; --ilayer) {
276  const DetLayer* layer = myLayers[ilayer];
277 
278  // allow the z of the hit to be within a straight line from a vertex
279  // of +-15 cm to the cluster
280  // float dphi = 0.015;
281  float dphi = 0.030;
282 
283  MeasurementEstimator* newEstimator = makeEstimator(layer, dphi);
284 
285  //std::cout << "OutInSeedFinder::startSeed propagationDirection " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";
286 
287  TSOS tsos(fts, layer->surface());
288 
289  LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed after TSOS tsos(fts, layer->surface() ) "
290  << "\n";
291 
292  LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
293  theFirstMeasurements_ =
294  theLayerMeasurements_.measurements(*layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
295 
296  //std::cout << "OutInSeedFinder::startSeed after theFirstMeasurements_ " << theFirstMeasurements_.size() << "\n";
297 
298  if (theFirstMeasurements_.size() > 1) // always a dummy returned, too
299  LogDebug("OutInConversionSeedFinder") << " Found " << theFirstMeasurements_.size() - 1 << " 1st hits in seed"
300  << "\n";
301 
302  delete newEstimator;
303 
304  LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer
305  << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
306 
307  for (unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
308  TrajectoryMeasurement m1 = theFirstMeasurements_[i];
309  if (m1.recHit()->isValid()) {
310  // update the fts to start from this point. much better than starting from
311  // extrapolated point along the line
312  GlobalPoint hitPoint = m1.recHit()->globalPosition();
313  LogDebug("OutInConversionSeedFinder")
314  << " Valid hit at R " << m1.recHit()->globalPosition().perp() << " Z "
315  << m1.recHit()->globalPosition().z() << " eta " << m1.recHit()->globalPosition().eta() << " phi "
316  << m1.recHit()->globalPosition().phi() << " xyz " << m1.recHit()->globalPosition() << "\n";
317 
318  FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
319  //std::cout << "OutInConversionSeedFinder::startSeed newfts " << newfts << "\n";
320  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed newfts " << newfts << "\n";
321  LogDebug("OutInConversionSeedFinder")
322  << "OutInConversionSeedFinder::startSeed propagationDirection after switching "
323  << int(thePropagatorOppositeToMomentum_->propagationDirection()) << "\n";
324  // std::cout << "OutInConversionSeedFinder::startSeed propagationDirection after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
325 
326  completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer - 1);
327  // skip a layer, if you haven't already skipped the first layer
328  if (ilayer == myLayers.size() - 1) {
329  completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer - 2);
330  }
331  }
332  }
333 
334  } // loop over layers
335  }
336 }
337 
338 MeasurementEstimator* OutInConversionSeedFinder::makeEstimator(const DetLayer* layer, float dphi) const {
339  //std::cout << "OutInConversionSeedFinder::makeEstimator " << "\n";
340 
341  MeasurementEstimator* newEstimator = nullptr;
342 
343  if (layer->location() == GeomDetEnumerators::barrel) {
344  const BarrelDetLayer* barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
345  LogDebug("OutInConversionSeedFinder")
346  << "OutInConversionSeedFinder::makeEstimator Barrel r = " << barrelLayer->specificSurface().radius() << " "
347  << "\n";
348  float r = barrelLayer->specificSurface().radius();
349  float zrange = 15. * (1. - r / theBCPosition_.perp());
350  newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
351  }
352 
353  if (layer->location() == GeomDetEnumerators::endcap) {
354  const ForwardDetLayer* forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
355  LogDebug("OutInConversionSeedFinder")
356  << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius()
357  << " R " << forwardLayer->specificSurface().outerRadius() << " Z "
358  << forwardLayer->specificSurface().position().z() << "\n";
359 
360  float zc = fabs(theBCPosition_.z());
361  float z = fabs(forwardLayer->surface().position().z());
362 
363  float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc * zc - 15. * zc);
364  newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
365  }
366 
367  return newEstimator;
368 }
369 
370 void OutInConversionSeedFinder::completeSeed(const TrajectoryMeasurement& m1,
371  const FreeTrajectoryState& fts,
372  const Propagator* propagator,
373  int ilayer) {
374  //std::cout << "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
375 
376  MeasurementEstimator* newEstimator = nullptr;
377  const DetLayer* layer = theLayerList_[ilayer];
378 
379  if (layer->location() == GeomDetEnumerators::barrel) {
380  // z error for 2nd hit is 2 sigma quadded with 5 cm
381  LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_
382  << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
383  float dz = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1.recHit()->globalPositionError().czz() +
384  the2ndHitdzConst_ * the2ndHitdzConst_);
385  newEstimator = new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
386  } else {
387  LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_
388  << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
389  // z error for 2nd hit is 2sigma quadded with 5 cm
390  //float m1dr = m1.recHit().globalPositionError().rerr(m1.recHit().globalPosition());
391  float m1dr = sqrt(m1.recHit()->localPositionError().yy());
392  float dr = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1dr * m1dr + the2ndHitdzConst_ * the2ndHitdznSigma_);
393  // LogDebug("OutInConversionSeedFinder") << "second hit forward dr " << dr << " this hit " << m1dr << endl;
394  newEstimator = new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
395  }
396 
397  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
398 
399  // Get the measurements consistent with the FTS and the Estimator
400  TSOS tsos(fts, layer->surface());
401  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection "
402  << int(propagator->propagationDirection()) << "\n";
403  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator "
404  << newEstimator << "\n";
405 
406  LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
407  std::vector<TrajectoryMeasurement> measurements =
408  theLayerMeasurements_.measurements(*layer, tsos, *propagator, *newEstimator);
409  //std::cout << "OutInConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
410  delete newEstimator;
411 
412  for (unsigned int i = 0; i < measurements.size(); ++i) {
413  if (measurements[i].recHit()->isValid()) {
414  createSeed(m1, measurements[i]);
415  }
416  }
417 
418  //LogDebug("OutInConversionSeedFinder") << "COMPLETED " << theSeeds_.size() << " SEEDS " << "\n";
419 }
420 
421 void OutInConversionSeedFinder::createSeed(const TrajectoryMeasurement& m1, const TrajectoryMeasurement& m2) {
422  //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";
423 
424  FreeTrajectoryState fts = createSeedFTS(m1, m2);
425 
426  //std::cout << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
427  // std::cout << "original cluster FTS " << fts <<"\n";
428  LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors "
429  << m1.recHit()->parametersError() << "\n";
430  LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts << "\n";
431 
432  //std::cout << "OutInConversionSeedFinder::createSeed propagation dir " << int( thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
433  TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts, m1.recHit()->det()->surface());
434 
435  // LogDebug("OutInConversionSeedFinder") << "hit surface " << h1.det().surface().position() << endl;
436  // LogDebug("OutInConversionSeedFinder") << "prop to " << typeid(h1.det().surface()).name() << endl;
437  // LogDebug("OutInConversionSeedFinder") << "prop to first hit " << state1 << endl;
438  // LogDebug("OutInConversionSeedFinder") << "update to " << h1.globalPosition() << endl;
439 
440  if (state1.isValid()) {
441  TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit());
442 
443  if (updatedState1.isValid()) {
444  TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(
445  *updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
446 
447  if (state2.isValid()) {
448  TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit());
449  TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit(), m1.estimate(), m1.layer());
450  TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit(), m2.estimate(), m2.layer());
451 
452  edm::OwnVector<TrackingRecHit> myHits;
453  myHits.push_back(meas1.recHit()->hit()->clone());
454  myHits.push_back(meas2.recHit()->hit()->clone());
455 
456  if (nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_)
457  return;
458 
459  PTrajectoryStateOnDet ptsod =
460  trajectoryStateTransform::persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId());
461 
462  LogDebug("OutInConversionSeedFinder")
463  << "OutInConversionSeedFinder::createSeed new seed from state " << state2.globalPosition() << "\n";
464  LogDebug("OutInConversionSeedFinder")
465  << "OutInConversionSeedFinder::createSeed new seed ptsod " << ptsod.parameters().position() << " R "
466  << ptsod.parameters().position().perp() << " phi " << ptsod.parameters().position().phi() << " eta "
467  << ptsod.parameters().position().eta() << "\n";
468 
469  theSeeds_.push_back(TrajectorySeed(ptsod, myHits, oppositeToMomentum));
470  nSeedsPerBC_++;
471  }
472  }
473  }
474 }
475 
476 FreeTrajectoryState OutInConversionSeedFinder::createSeedFTS(const TrajectoryMeasurement& m1,
477  const TrajectoryMeasurement& m2) const {
478  GlobalPoint xmeas = fixPointRadius(m1);
479  GlobalPoint xvert = fixPointRadius(m2);
480 
481  float pt = theSCenergy_ * sin(theSCPosition_.theta());
482  float pz = theSCenergy_ * cos(theSCPosition_.theta());
483 
484  // doesn't work at all for endcap, where r is badly measured
485  //float dphidr = (p1.phi()-p2.phi())/(p1.perp()-p2.perp());
486  //int charge = (dphidr > 0.) ? -1 : 1;
487  int charge = m2.predictedState().charge();
488 
489  double BInTesla = theMF_->inTesla(xmeas).z();
490  GlobalVector xdiff = xmeas - xvert;
491 
492  double phi = xdiff.phi();
493  double pxOld = pt * cos(phi);
494  double pyOld = pt * sin(phi);
495  double RadCurv = 100 * pt / (BInTesla * 0.29979);
496  double alpha = asin(0.5 * xdiff.perp() / RadCurv);
497  float ca = cos(charge * alpha);
498  float sa = sin(charge * alpha);
499  double pxNew = ca * pxOld + sa * pyOld;
500  double pyNew = -sa * pxOld + ca * pyOld;
501  GlobalVector pNew(pxNew, pyNew, pz);
502 
503  GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_));
504 
505  AlgebraicSymMatrix55 m = AlgebraicMatrixID();
506  m(0, 0) = 0.05;
507  m(1, 1) = 0.02;
508  m(2, 2) = 0.007;
509  m(3, 3) = 10.;
510  m(4, 4) = 10.;
511  //std::cout << "OutInConversionSeedFinder::createSeedFTS " << FreeTrajectoryState(gp, CurvilinearTrajectoryError(m)) << "\n";
512  return FreeTrajectoryState(gp, CurvilinearTrajectoryError(m));
513 }
514 
515 GlobalPoint OutInConversionSeedFinder::fixPointRadius(const TrajectoryMeasurement& m1) const {
516  GlobalPoint p1 = m1.recHit()->globalPosition();
517  GlobalPoint p2;
518  if (m1.layer()->location() == GeomDetEnumerators::barrel) {
519  p2 = p1;
520  } else {
521  float z = p1.z();
522  float phi = p1.phi();
523  float theta = theSCPosition_.theta();
524  float r = p1.z() * tan(theta);
525  p2 = GlobalPoint(r * cos(phi), r * sin(phi), z);
526  // LogDebug("OutInConversionSeedFinder") << "fixing point radius " << p2 << " from " << p1 << endl;
527  }
528  return p2;
529 }
ConversionSeedFinder::layerList
std::vector< const DetLayer * > const & layerList() const
Definition: ConversionSeedFinder.h:61
Vector3DBase
Definition: Vector3DBase.h:8
ConversionSeedFinder::theBeamSpot_
reco::BeamSpot theBeamSpot_
Definition: ConversionSeedFinder.h:112
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:122
mps_fire.i
i
Definition: mps_fire.py:355
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:189
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:515
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
reco::deltaPhi
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
zMuMuMuonUserData.alpha
alpha
zGenParticlesMatch = cms.InputTag(""),
Definition: zMuMuMuonUserData.py:9
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:113
OutInConversionSeedFinder::startSeed
void startSeed(const FreeTrajectoryState &)
Definition: OutInConversionSeedFinder.cc:269
TrajectoryStateOnSurface::charge
TrackCharge charge() const
Definition: TrajectoryStateOnSurface.h:68
oppositeToMomentum
Definition: PropagationDirection.h:4
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:215
OutInConversionSeedFinder::createSeed
void createSeed(const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2)
Definition: OutInConversionSeedFinder.cc:421
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:93
edm::Handle
Definition: AssociativeIterator.h:50
ConversionSeedFinder::thePropagatorOppositeToMomentum_
const Propagator * thePropagatorOppositeToMomentum_
Definition: ConversionSeedFinder.h:110
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:370
LayerMeasurements
Definition: LayerMeasurements.h:18
TrajectoryMeasurement::predictedState
TrajectoryStateOnSurface const & predictedState() const
Definition: TrajectoryMeasurement.h:174
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:72
OutInConversionSeedFinder::makeEstimator
MeasurementEstimator * makeEstimator(const DetLayer *, float dphi) const
Definition: OutInConversionSeedFinder.cc:338
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:99
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
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
ConversionSeedFinder::theBCEnergy_
float theBCEnergy_
Definition: ConversionSeedFinder.h:107
runTauDisplay.gp
gp
Definition: runTauDisplay.py:431
edm::View
Definition: CaloClusterFwd.h:14
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
edm::ParameterSet
Definition: ParameterSet.h:36
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:102
OutInConversionSeedFinder::the2ndHitdzConst_
float the2ndHitdzConst_
Definition: OutInConversionSeedFinder.h:58
GeomDetEnumerators::endcap
Definition: GeomDetEnumerators.h:9
OutInConversionSeedFinder.h
ConversionSeedFinder::theSCPosition_
GlobalPoint theSCPosition_
Definition: ConversionSeedFinder.h:86
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:109
p1
double p1[4]
Definition: TauolaWrapper.h:89
ConversionSeedFinder::getMeasurementTracker
const MeasurementTracker * getMeasurementTracker() const
Definition: ConversionSeedFinder.h:64
MagneticField.h
ConversionSeedFinder::theSeeds_
TrajectorySeedCollection theSeeds_
Definition: ConversionSeedFinder.h:85
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:47
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:476
TrajectoryMeasurement::recHit
ConstRecHitPointer const & recHit() const
Definition: TrajectoryMeasurement.h:190
DDAxes::phi
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
ConversionSeedFinder::trackStateFromClusters
FreeTrajectoryState trackStateFromClusters(int aCharge, const GlobalPoint &gpOrigine, PropagationDirection dir, float scaleFactor) const
Definition: ConversionSeedFinder.cc:70
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
ConversionSeedFinder::theBCPosition_
GlobalPoint theBCPosition_
Definition: ConversionSeedFinder.h:106
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:50
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
GeometricSearchDet::surface
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
DetLayer.h
TrajectorySeed
Definition: TrajectorySeed.h:17
ConversionSeedFinder::theLayerList_
std::vector< const DetLayer * > theLayerList_
Definition: ConversionSeedFinder.h:104
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:303
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:39
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
deltaPhi.h
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