CMS 3D CMS Logo
PFTrackTransformer.cc
Go to the documentation of this file.
1 //
2 // -*- C++ -*-
3 // Package: PFTracking
4 // Class: PFTrackTransformer
5 //
6 // Original Author: Michele Pioppi
7 // Other Author: Daniele Benedetti
8 
9 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
10 
11 #include "FWCore/MessageLogger/interface/MessageLogger.h"
12 #include "FWCore/Framework/interface/ESHandle.h"
13 
14 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
15 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
16 #include "DataFormats/TrackReco/interface/Track.h"
17 
18 #include "TrackingTools/PatternTools/interface/Trajectory.h"
19 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
20 
21 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
22 // Add by Daniele
23 #include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
24 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
25 #include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"
26 #include "RecoParticleFlow/PFTracking/interface/PFGsfHelper.h"
27 
28 using namespace std;
29 using namespace reco;
30 using namespace edm;
31 
32 PFTrackTransformer::PFTrackTransformer(const math::XYZVector& B) : B_(B) {
33  LogInfo("PFTrackTransformer") << "PFTrackTransformer built";
34 
35  onlyprop_ = false;
36 }
37 
38 PFTrackTransformer::~PFTrackTransformer() {}
39 
40 bool PFTrackTransformer::addPoints(reco::PFRecTrack& pftrack,
41  const reco::Track& track,
42  const Trajectory& traj,
43  bool msgwarning) const {
44  LogDebug("PFTrackTransformer") << "Trajectory propagation started";
45  using namespace reco;
46  using namespace std;
47 
48  float PT = track.pt();
49  float pfmass = (pftrack.algoType() == reco::PFRecTrack::KF_ELCAND) ? 0.0005 : 0.139;
50  float pfenergy = sqrt((pfmass * pfmass) + (track.p() * track.p()));
51  // closest approach
52  BaseParticlePropagator theParticle = BaseParticlePropagator(
53  RawParticle(XYZTLorentzVector(track.px(), track.py(), track.pz(), pfenergy),
54  XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.),
55  track.charge()),
56  0.,
57  0.,
58  B_.z());
59 
60  float pfoutenergy = sqrt((pfmass * pfmass) + track.outerMomentum().Mag2());
61  BaseParticlePropagator theOutParticle;
62  if (track.outerOk())
63  theOutParticle = BaseParticlePropagator(
64  RawParticle(
65  XYZTLorentzVector(
66  track.outerMomentum().x(), track.outerMomentum().y(), track.outerMomentum().z(), pfoutenergy),
67  XYZTLorentzVector(track.outerPosition().x(), track.outerPosition().y(), track.outerPosition().z(), 0.),
68  track.charge()),
69  0.,
70  0.,
71  B_.z());
72  else
73  theOutParticle = theParticle; // if outer state is not available, use the information at the closest approach
74 
75  math::XYZTLorentzVector momClosest = math::XYZTLorentzVector(track.px(), track.py(), track.pz(), track.p());
76  const math::XYZPoint& posClosest = track.vertex();
77 
78  pftrack.addPoint(PFTrajectoryPoint(-1, PFTrajectoryPoint::ClosestApproach, posClosest, momClosest));
79 
80  //BEAMPIPE
81  theParticle.setPropagationConditions(
82  pfGeometry_.outerRadius(PFGeometry::BeamPipe), pfGeometry_.outerZ(PFGeometry::BeamPipe), false);
83  theParticle.propagate();
84  if (theParticle.getSuccess() != 0)
85  pftrack.addPoint(PFTrajectoryPoint(-1,
86  PFTrajectoryPoint::BeamPipeOrEndVertex,
87  math::XYZPoint(theParticle.particle().vertex()),
88  math::XYZTLorentzVector(theParticle.particle().momentum())));
89  else {
90  PFTrajectoryPoint dummyMaxSh;
91  pftrack.addPoint(dummyMaxSh);
92  }
93 
94  //trajectory points
95 
96  if (!onlyprop_) {
97  bool direction = (traj.direction() == alongMomentum);
98  vector<TrajectoryMeasurement> measurements = traj.measurements();
99  int iTrajFirst = (direction) ? 0 : measurements.size() - 1;
100  int increment = (direction) ? +1 : -1;
101  int iTrajLast = (direction) ? int(measurements.size()) : -1;
102 
103  for (int iTraj = iTrajFirst; iTraj != iTrajLast; iTraj += increment) {
104  GlobalPoint v = measurements[iTraj].updatedState().globalPosition();
105  GlobalVector p = measurements[iTraj].updatedState().globalMomentum();
106  unsigned int iid = measurements[iTraj].recHit()->det()->geographicalId().rawId();
107  pftrack.addPoint(PFTrajectoryPoint(
108  iid, -1, math::XYZPoint(v.x(), v.y(), v.z()), math::XYZTLorentzVector(p.x(), p.y(), p.z(), p.mag())));
109  }
110  }
111 
112  // ES1
113  bool isBelowPS = false;
114  theOutParticle.propagateToPreshowerLayer1(false);
115  if (theOutParticle.getSuccess() != 0)
116  pftrack.addPoint(PFTrajectoryPoint(-1,
117  PFTrajectoryPoint::PS1,
118  math::XYZPoint(theOutParticle.particle().vertex()),
119  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
120  else {
121  PFTrajectoryPoint dummyPS1;
122  pftrack.addPoint(dummyPS1);
123  }
124 
125  // ES1
126  theOutParticle.propagateToPreshowerLayer2(false);
127  if (theOutParticle.getSuccess() != 0) {
128  pftrack.addPoint(PFTrajectoryPoint(-1,
129  PFTrajectoryPoint::PS2,
130  math::XYZPoint(theOutParticle.particle().vertex()),
131  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
132  isBelowPS = true;
133  } else {
134  PFTrajectoryPoint dummyPS2;
135  pftrack.addPoint(dummyPS2);
136  }
137 
138  //ECAL entrance
139  theOutParticle.propagateToEcalEntrance(false);
140  if (theOutParticle.getSuccess() != 0) {
141  pftrack.addPoint(PFTrajectoryPoint(-1,
142  PFTrajectoryPoint::ECALEntrance,
143  math::XYZPoint(theOutParticle.particle().vertex()),
144  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
145  double ecalShowerDepth = PFCluster::getDepthCorrection(theOutParticle.particle().momentum().E(), isBelowPS, false);
146 
147  math::XYZPoint meanShower =
148  math::XYZPoint(theOutParticle.particle().vertex()) +
149  math::XYZTLorentzVector(theOutParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
150 
151  pftrack.addPoint(PFTrajectoryPoint(-1,
152  PFTrajectoryPoint::ECALShowerMax,
153  meanShower,
154  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
155  } else {
156  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_ && msgwarning)
157  // Check consistent boundary as in CommonTools/BaseParticlePropagator/src/BaseParticlePropagator.cc
158  // out of the HB/HE acceptance
159  // eta = 3.0 -> cos^2(theta) = 0.99014 - cos**2(theta) is faster to determine than eta
160  LogWarning("PFTrackTransformer") << "KF TRACK " << pftrack << " PROPAGATION TO THE ECAL HAS FAILED\n"
161  << "theOutParticle.particle() pt,eta: " << theOutParticle.particle().pt() << " "
162  << theOutParticle.particle().eta();
163  PFTrajectoryPoint dummyECAL;
164  pftrack.addPoint(dummyECAL);
165  PFTrajectoryPoint dummyMaxSh;
166  pftrack.addPoint(dummyMaxSh);
167  }
168 
169  //HCAL entrance
170  theOutParticle.propagateToHcalEntrance(false);
171  if (theOutParticle.getSuccess() != 0)
172  pftrack.addPoint(PFTrajectoryPoint(-1,
173  PFTrajectoryPoint::HCALEntrance,
174  math::XYZPoint(theOutParticle.particle().vertex()),
175  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
176  else {
177  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_ && msgwarning)
178  LogWarning("PFTrackTransformer") << "KF TRACK " << pftrack << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
179  PFTrajectoryPoint dummyHCALentrance;
180  pftrack.addPoint(dummyHCALentrance);
181  }
182 
183  //HCAL exit
184  // theOutParticle.setMagneticField(0); //Show we propagate as straight line inside HCAL ?
185  theOutParticle.propagateToHcalExit(false);
186  if (theOutParticle.getSuccess() != 0)
187  pftrack.addPoint(PFTrajectoryPoint(-1,
188  PFTrajectoryPoint::HCALExit,
189  math::XYZPoint(theOutParticle.particle().vertex()),
190  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
191  else {
192  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_ && msgwarning)
193  LogWarning("PFTrackTransformer") << "KF TRACK " << pftrack << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
194  PFTrajectoryPoint dummyHCALexit;
195  pftrack.addPoint(dummyHCALexit);
196  }
197 
198  //HO layer0
199  // if (abs(theOutParticle.particle().vertex().z())<550) {
200  if (PT > 3.0) { //Same value is used in PFBlockAlgo::link( case PFBlockLink::TRACKandHO:
201  theOutParticle.setMagneticField(0);
202  theOutParticle.propagateToHOLayer(false);
203  if (theOutParticle.getSuccess() != 0) {
204  pftrack.addPoint(PFTrajectoryPoint(-1,
205  PFTrajectoryPoint::HOLayer,
206  math::XYZPoint(theOutParticle.particle().vertex()),
207  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
208  } else {
209  if (PT > 5. && abs(theOutParticle.particle().Z()) < 700.25 && msgwarning)
210  LogWarning("PFTrackTransformer") << "KF TRACK " << pftrack << " PROPAGATION TO THE HO HAS FAILED";
211  PFTrajectoryPoint dummyHOLayer;
212  pftrack.addPoint(dummyHOLayer);
213  }
214  } else {
215  PFTrajectoryPoint dummyHOLayer;
216  pftrack.addPoint(dummyHOLayer);
217  }
218 
219  //VFcal(HF) entrance
220  theOutParticle.setMagneticField(0); // assume B=0 works ok for extrapolation to HF
221  theOutParticle.propagateToVFcalEntrance(false);
222  if (theOutParticle.getSuccess() != 0) {
223  pftrack.addPoint(PFTrajectoryPoint(-1,
224  PFTrajectoryPoint::VFcalEntrance,
225  math::XYZPoint(theOutParticle.particle().vertex()),
226  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
227  } else {
228  PFTrajectoryPoint dummyVFcalentrance;
229  pftrack.addPoint(dummyVFcalentrance);
230  }
231 
232  return true;
233 }
234 bool PFTrackTransformer::addPointsAndBrems(reco::GsfPFRecTrack& pftrack,
235  const reco::Track& track,
236  const Trajectory& traj,
237  const bool& GetMode) const {
238  float PT = track.pt();
239  // Trajectory for each trajectory point
240 
241  bool direction = (traj.direction() == alongMomentum);
242  vector<TrajectoryMeasurement> measurements = traj.measurements();
243  int iTrajFirst = (direction) ? 0 : measurements.size() - 1;
244  int increment = (direction) ? +1 : -1;
245  int iTrajLast = (direction) ? int(measurements.size()) : -1;
246 
247  unsigned int iTrajPos = 0;
248  for (int iTraj = iTrajFirst; iTraj != iTrajLast; iTraj += increment) {
249  GlobalPoint v = measurements[iTraj].updatedState().globalPosition();
250  PFGsfHelper* PFGsf = new PFGsfHelper(measurements[iTraj]);
251  //if (PFGsf->isValid()){
252  bool ComputeMODE = GetMode;
253  GlobalVector p = PFGsf->computeP(ComputeMODE);
254  double DP = PFGsf->fittedDP();
255  double SigmaDP = PFGsf->sigmafittedDP();
256  unsigned int iid = measurements[iTraj].recHit()->det()->geographicalId().rawId();
257  delete PFGsf;
258 
259  // -------------------------- Fill GSF Track -------------------------------------
260 
261  // float pfmass= (pftrack.algoType()==reco::PFRecTrack::KF_ELCAND) ? 0.0005 : 0.139;
262  float ptot = sqrt((p.x() * p.x()) + (p.y() * p.y()) + (p.z() * p.z()));
263  float pfenergy = ptot;
264 
265  if (iTraj == iTrajFirst) {
266  math::XYZTLorentzVector momClosest = math::XYZTLorentzVector(p.x(), p.y(), p.z(), ptot);
267  const math::XYZPoint& posClosest = track.vertex();
268  pftrack.addPoint(PFTrajectoryPoint(-1, PFTrajectoryPoint::ClosestApproach, posClosest, momClosest));
269 
270  BaseParticlePropagator theInnerParticle = BaseParticlePropagator(
271  RawParticle(XYZTLorentzVector(p.x(), p.y(), p.z(), pfenergy),
272  XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.),
273  track.charge()), //DANIELE Same thing v.x(),v.y(),v.()?
274  0.,
275  0.,
276  B_.z());
277 
278  //BEAMPIPE
279  theInnerParticle.setPropagationConditions(
280  pfGeometry_.outerRadius(PFGeometry::BeamPipe), pfGeometry_.outerZ(PFGeometry::BeamPipe), false);
281  theInnerParticle.propagate();
282  if (theInnerParticle.getSuccess() != 0)
283  pftrack.addPoint(PFTrajectoryPoint(-1,
284  PFTrajectoryPoint::BeamPipeOrEndVertex,
285  math::XYZPoint(theInnerParticle.particle().vertex()),
286  math::XYZTLorentzVector(theInnerParticle.particle().momentum())));
287  else {
288  PFTrajectoryPoint dummyMaxSh;
289  pftrack.addPoint(dummyMaxSh);
290  }
291 
292  // First Point for the trajectory == Vertex ??
293  pftrack.addPoint(PFTrajectoryPoint(
294  iid, -1, math::XYZPoint(v.x(), v.y(), v.z()), math::XYZTLorentzVector(p.x(), p.y(), p.z(), p.mag())));
295  }
296  if (iTraj != iTrajFirst && iTraj != (abs(iTrajLast) - 1)) {
297  pftrack.addPoint(PFTrajectoryPoint(
298  iid, -1, math::XYZPoint(v.x(), v.y(), v.z()), math::XYZTLorentzVector(p.x(), p.y(), p.z(), p.mag())));
299  }
300  if (iTraj == (abs(iTrajLast) - 1)) {
301  // Last Trajectory Meas
302  pftrack.addPoint(PFTrajectoryPoint(
303  iid, -1, math::XYZPoint(v.x(), v.y(), v.z()), math::XYZTLorentzVector(p.x(), p.y(), p.z(), p.mag())));
304 
305  BaseParticlePropagator theOutParticle =
306  BaseParticlePropagator(RawParticle(XYZTLorentzVector(p.x(), p.y(), p.z(), pfenergy),
307  XYZTLorentzVector(v.x(), v.y(), v.z(), 0.),
308  track.charge()),
309  0.,
310  0.,
311  B_.z());
312 
313  // ES1
314  bool isBelowPS = false;
315  theOutParticle.propagateToPreshowerLayer1(false);
316  if (theOutParticle.getSuccess() != 0)
317  pftrack.addPoint(PFTrajectoryPoint(-1,
318  PFTrajectoryPoint::PS1,
319  math::XYZPoint(theOutParticle.particle().vertex()),
320  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
321  else {
322  PFTrajectoryPoint dummyPS1;
323  pftrack.addPoint(dummyPS1);
324  }
325 
326  // ES2
327  theOutParticle.propagateToPreshowerLayer2(false);
328  if (theOutParticle.getSuccess() != 0) {
329  pftrack.addPoint(PFTrajectoryPoint(-1,
330  PFTrajectoryPoint::PS2,
331  math::XYZPoint(theOutParticle.particle().vertex()),
332  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
333  isBelowPS = true;
334  } else {
335  PFTrajectoryPoint dummyPS2;
336  pftrack.addPoint(dummyPS2);
337  }
338 
339  theOutParticle.propagateToEcalEntrance(false);
340 
341  if (theOutParticle.getSuccess() != 0) {
342  pftrack.addPoint(PFTrajectoryPoint(-1,
343  PFTrajectoryPoint::ECALEntrance,
344  math::XYZPoint(theOutParticle.particle().vertex()),
345  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
346  double ecalShowerDepth =
347  PFCluster::getDepthCorrection(theOutParticle.particle().momentum().E(), isBelowPS, false);
348 
349  math::XYZPoint meanShower =
350  math::XYZPoint(theOutParticle.particle().vertex()) +
351  math::XYZTLorentzVector(theOutParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
352 
353  pftrack.addPoint(PFTrajectoryPoint(-1,
354  PFTrajectoryPoint::ECALShowerMax,
355  meanShower,
356  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
357  } else {
358  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
359  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE ECAL HAS FAILED\n"
360  << "theOutParticle.particle() pt,eta: " << theOutParticle.particle().pt()
361  << " " << theOutParticle.particle().eta();
362  PFTrajectoryPoint dummyECAL;
363  pftrack.addPoint(dummyECAL);
364  PFTrajectoryPoint dummyMaxSh;
365  pftrack.addPoint(dummyMaxSh);
366  }
367 
368  //HCAL entrance
369  theOutParticle.propagateToHcalEntrance(false);
370  if (theOutParticle.getSuccess() != 0)
371  pftrack.addPoint(PFTrajectoryPoint(-1,
372  PFTrajectoryPoint::HCALEntrance,
373  math::XYZPoint(theOutParticle.particle().vertex()),
374  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
375  else {
376  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
377  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
378  PFTrajectoryPoint dummyHCALentrance;
379  pftrack.addPoint(dummyHCALentrance);
380  }
381  //HCAL exit
382  theOutParticle.propagateToHcalExit(false);
383  if (theOutParticle.getSuccess() != 0)
384  pftrack.addPoint(PFTrajectoryPoint(-1,
385  PFTrajectoryPoint::HCALExit,
386  math::XYZPoint(theOutParticle.particle().vertex()),
387  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
388  else {
389  if (PT > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
390  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
391  PFTrajectoryPoint dummyHCALexit;
392  pftrack.addPoint(dummyHCALexit);
393  }
394 
395  //HO Layer0
396  if (abs(theOutParticle.particle().vertex().z()) < 550) {
397  theOutParticle.setMagneticField(0);
398  theOutParticle.propagateToHOLayer(false);
399  if (theOutParticle.getSuccess() != 0)
400  pftrack.addPoint(PFTrajectoryPoint(-1,
401  PFTrajectoryPoint::HOLayer,
402  math::XYZPoint(theOutParticle.particle().vertex()),
403  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
404  else {
405  if (PT > 5. && abs(theOutParticle.particle().Z()) < 700.25)
406  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HO HAS FAILED";
407  PFTrajectoryPoint dummyHOLayer;
408  pftrack.addPoint(dummyHOLayer);
409  }
410  }
411  }
412 
413  // -------------------------- END GSF Track -------------------------------------
414 
415  // -------------------------- Fill Brem "Track" ---------------------------------
416  // Fill the brem for each traj point
417 
418  //check that the vertex of the brem is in the tracker volume
419  if ((v.perp() > 110) || (fabs(v.z()) > 280))
420  continue;
421  unsigned int iTrajPoint = iTrajPos + 2;
422  if (iid % 2 == 1)
423  iTrajPoint = 99;
424 
425  PFBrem brem(DP, SigmaDP, iTrajPoint);
426 
427  GlobalVector p_gamma = p * (fabs(DP) / p.mag()); // Direction from the electron (tangent), DP without any sign!;
428  float e_gamma = fabs(DP); // DP = pout-pin so could be negative
429  constexpr int gamma_charge = 0;
430  BaseParticlePropagator theBremParticle =
431  BaseParticlePropagator(RawParticle(XYZTLorentzVector(p_gamma.x(), p_gamma.y(), p_gamma.z(), e_gamma),
432  XYZTLorentzVector(v.x(), v.y(), v.z(), 0.),
433  gamma_charge),
434  0.,
435  0.,
436  B_.z());
437 
438  // add TrajectoryPoint for Brem, PS, ECAL, ECALShowMax, HCAL
439  // Brem Entrance PS Layer1
440 
441  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
442  brem.addPoint(dummyClosest);
443 
444  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
445  brem.addPoint(dummyBeamPipe);
446 
447  // ES1
448  bool isBelowPS = false;
449  theBremParticle.propagateToPreshowerLayer1(false);
450  if (theBremParticle.getSuccess() != 0)
451  brem.addPoint(PFTrajectoryPoint(-1,
452  PFTrajectoryPoint::PS1,
453  math::XYZPoint(theBremParticle.particle().vertex()),
454  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
455  else {
456  PFTrajectoryPoint dummyPS1;
457  brem.addPoint(dummyPS1);
458  }
459 
460  // ES2
461  // Brem Entrance PS Layer 2
462  theBremParticle.propagateToPreshowerLayer2(false);
463  if (theBremParticle.getSuccess() != 0) {
464  brem.addPoint(PFTrajectoryPoint(-1,
465  PFTrajectoryPoint::PS2,
466  math::XYZPoint(theBremParticle.particle().vertex()),
467  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
468  isBelowPS = true;
469  } else {
470  PFTrajectoryPoint dummyPS2;
471  brem.addPoint(dummyPS2);
472  }
473 
474  //ECAL entrance
475  theBremParticle.propagateToEcalEntrance(false);
476  if (theBremParticle.getSuccess() != 0) {
477  brem.addPoint(PFTrajectoryPoint(-1,
478  PFTrajectoryPoint::ECALEntrance,
479  math::XYZPoint(theBremParticle.particle().vertex()),
480  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
481  double ecalShowerDepth =
482  PFCluster::getDepthCorrection(theBremParticle.particle().momentum().E(), isBelowPS, false);
483 
484  math::XYZPoint meanShower =
485  math::XYZPoint(theBremParticle.particle().vertex()) +
486  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
487 
488  brem.addPoint(PFTrajectoryPoint(-1,
489  PFTrajectoryPoint::ECALShowerMax,
490  meanShower,
491  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
492  } else {
493  if ((DP > 5.) && ((DP / SigmaDP) > 3) && theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
494  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
495  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
496  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
497  << " " << theBremParticle.particle().cos2ThetaV();
498  PFTrajectoryPoint dummyECAL;
499  brem.addPoint(dummyECAL);
500  PFTrajectoryPoint dummyMaxSh;
501  brem.addPoint(dummyMaxSh);
502  }
503 
504  //HCAL entrance
505  theBremParticle.propagateToHcalEntrance(false);
506  if (theBremParticle.getSuccess() != 0)
507  brem.addPoint(PFTrajectoryPoint(-1,
508  PFTrajectoryPoint::HCALEntrance,
509  math::XYZPoint(theBremParticle.particle().vertex()),
510  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
511  else {
512  if ((DP > 5.) && ((DP / SigmaDP) > 3) && theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
513  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
514  PFTrajectoryPoint dummyHCALentrance;
515  brem.addPoint(dummyHCALentrance);
516  }
517 
518  //HCAL exit
519  theBremParticle.propagateToHcalExit(false);
520  if (theBremParticle.getSuccess() != 0)
521  brem.addPoint(PFTrajectoryPoint(-1,
522  PFTrajectoryPoint::HCALExit,
523  math::XYZPoint(theBremParticle.particle().vertex()),
524  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
525  else {
526  if ((DP > 5.) && ((DP / SigmaDP) > 3) && theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
527  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
528  PFTrajectoryPoint dummyHCALexit;
529  brem.addPoint(dummyHCALexit);
530  }
531 
532  //HO Layer0
533  if (abs(theBremParticle.particle().vertex().z()) < 550.0) {
534  theBremParticle.setMagneticField(0);
535  theBremParticle.propagateToHOLayer(false);
536  if (theBremParticle.getSuccess() != 0)
537  brem.addPoint(PFTrajectoryPoint(-1,
538  PFTrajectoryPoint::HOLayer,
539  math::XYZPoint(theBremParticle.particle().vertex()),
540  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
541  else {
542  if ((DP > 5.) && ((DP / SigmaDP) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
543  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
544  PFTrajectoryPoint dummyHOLayer;
545  brem.addPoint(dummyHOLayer);
546  }
547  }
548  pftrack.addBrem(brem);
549  iTrajPos++;
550  }
551  return true;
552 }
553 
554 bool PFTrackTransformer::addPointsAndBrems(reco::GsfPFRecTrack& pftrack,
555  const reco::GsfTrack& track,
556  const MultiTrajectoryStateTransform& mtjstate) const {
557  // float PT= track.pt();
558  unsigned int iTrajPos = 0;
559  unsigned int iid = 0; // not anymore saved
560 
561  // ***************************** INNER State *************************************
562  TrajectoryStateOnSurface inTSOS = mtjstate.innerStateOnSurface((track));
563  TrajectoryStateOnSurface outTSOS = mtjstate.outerStateOnSurface((track));
564 
565  if (!inTSOS.isValid() || !outTSOS.isValid()) {
566  if (!inTSOS.isValid())
567  LogWarning("PFTrackTransformer") << " INNER TSOS NOT VALID ";
568  if (!outTSOS.isValid())
569  LogWarning("PFTrackTransformer") << " OUTER TSOS NOT VALID ";
570  return false;
571  }
572 
573  GlobalVector InMom;
574  GlobalPoint InPos;
575  if (inTSOS.isValid()) {
576  multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS, InMom);
577  multiTrajectoryStateMode::positionFromModeCartesian(inTSOS, InPos);
578  } else {
579  InMom = GlobalVector(track.pxMode(), track.pyMode(), track.pzMode());
580  InPos = GlobalPoint(0., 0., 0.);
581  }
582 
583  // float pfmass= (pftrack.algoType()==reco::PFRecTrack::KF_ELCAND) ? 0.0005 : 0.139;
584  float ptot = sqrt((InMom.x() * InMom.x()) + (InMom.y() * InMom.y()) + (InMom.z() * InMom.z()));
585  float pfenergy = ptot;
586 
587  math::XYZTLorentzVector momClosest = math::XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), ptot);
588  const math::XYZPoint& posClosest = track.vertex();
589  pftrack.addPoint(PFTrajectoryPoint(-1, PFTrajectoryPoint::ClosestApproach, posClosest, momClosest));
590 
591  BaseParticlePropagator theInnerParticle = BaseParticlePropagator(
592  RawParticle(XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), pfenergy),
593  XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.),
594  track.charge()), //DANIELE Same thing v.x(),v.y(),v.()?
595  0.,
596  0.,
597  B_.z());
598  //BEAMPIPE
599  theInnerParticle.setPropagationConditions(
600  pfGeometry_.outerRadius(PFGeometry::BeamPipe), pfGeometry_.outerZ(PFGeometry::BeamPipe), false);
601  theInnerParticle.propagate();
602  if (theInnerParticle.getSuccess() != 0)
603  pftrack.addPoint(PFTrajectoryPoint(-1,
604  PFTrajectoryPoint::BeamPipeOrEndVertex,
605  math::XYZPoint(theInnerParticle.particle().vertex()),
606  math::XYZTLorentzVector(theInnerParticle.particle().momentum())));
607  else {
608  PFTrajectoryPoint dummyBeam;
609  pftrack.addPoint(dummyBeam);
610  }
611 
612  // first tjpoint
613  pftrack.addPoint(PFTrajectoryPoint(iid,
614  -1,
615  math::XYZPoint(InPos.x(), InPos.y(), InPos.z()),
616  math::XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), InMom.mag())));
617 
618  //######### Photon at INNER State ##########
619 
620  unsigned int iTrajPoint = iTrajPos + 2;
621  double dp_tang = ptot;
622  double sdp_tang = track.ptModeError() * (track.pMode() / track.ptMode());
623  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
624  constexpr int gamma_charge = 0;
625  BaseParticlePropagator theBremParticle =
626  BaseParticlePropagator(RawParticle(XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), dp_tang),
627  XYZTLorentzVector(InPos.x(), InPos.y(), InPos.z(), 0.),
628  gamma_charge),
629  0.,
630  0.,
631  B_.z());
632  // add TrajectoryPoint for Brem, PS, ECAL, ECALShowMax, HCAL
633  // Brem Entrance PS Layer1
634  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
635  brem.addPoint(dummyClosest);
636 
637  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
638  brem.addPoint(dummyBeamPipe);
639 
640  bool isBelowPS = false;
641  theBremParticle.propagateToPreshowerLayer1(false);
642  if (theBremParticle.getSuccess() != 0)
643  brem.addPoint(PFTrajectoryPoint(-1,
644  PFTrajectoryPoint::PS1,
645  math::XYZPoint(theBremParticle.particle().vertex()),
646  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
647  else {
648  PFTrajectoryPoint dummyPS1;
649  brem.addPoint(dummyPS1);
650  }
651 
652  // Brem Entrance PS Layer 2
653 
654  theBremParticle.propagateToPreshowerLayer2(false);
655  if (theBremParticle.getSuccess() != 0) {
656  brem.addPoint(PFTrajectoryPoint(-1,
657  PFTrajectoryPoint::PS2,
658  math::XYZPoint(theBremParticle.particle().vertex()),
659  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
660  isBelowPS = true;
661  } else {
662  PFTrajectoryPoint dummyPS2;
663  brem.addPoint(dummyPS2);
664  }
665 
666  //ECAL entrance
667  theBremParticle.propagateToEcalEntrance(false);
668  if (theBremParticle.getSuccess() != 0) {
669  brem.addPoint(PFTrajectoryPoint(-1,
670  PFTrajectoryPoint::ECALEntrance,
671  math::XYZPoint(theBremParticle.particle().vertex()),
672  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
673 
674  // for the first brem give a low default DP of 100 MeV.
675  double EDepthCorr = 0.01;
676  double ecalShowerDepth = PFCluster::getDepthCorrection(EDepthCorr, isBelowPS, false);
677 
678  math::XYZPoint meanShower =
679  math::XYZPoint(theBremParticle.particle().vertex()) +
680  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
681 
682  brem.addPoint(PFTrajectoryPoint(-1,
683  PFTrajectoryPoint::ECALShowerMax,
684  meanShower,
685  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
686  } else {
687  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
688  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
689  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
690  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
691  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
692  << " " << theBremParticle.particle().cos2ThetaV();
693  PFTrajectoryPoint dummyECAL;
694  brem.addPoint(dummyECAL);
695  PFTrajectoryPoint dummyMaxSh;
696  brem.addPoint(dummyMaxSh);
697  }
698 
699  //HCAL entrance
700  theBremParticle.propagateToHcalEntrance(false);
701  if (theBremParticle.getSuccess() != 0)
702  brem.addPoint(PFTrajectoryPoint(-1,
703  PFTrajectoryPoint::HCALEntrance,
704  math::XYZPoint(theBremParticle.particle().vertex()),
705  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
706  else {
707  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
708  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
709  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
710  PFTrajectoryPoint dummyHCALentrance;
711  brem.addPoint(dummyHCALentrance);
712  }
713 
714  //HCAL exit
715  theBremParticle.propagateToHcalExit(false);
716  if (theBremParticle.getSuccess() != 0)
717  brem.addPoint(PFTrajectoryPoint(-1,
718  PFTrajectoryPoint::HCALExit,
719  math::XYZPoint(theBremParticle.particle().vertex()),
720  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
721  else {
722  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
723  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
724  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
725  PFTrajectoryPoint dummyHCALexit;
726  brem.addPoint(dummyHCALexit);
727  }
728 
729  //HO Layer0
730  if (abs(theBremParticle.particle().vertex().z()) < 550) {
731  theBremParticle.setMagneticField(0);
732  theBremParticle.propagateToHOLayer(false);
733  if (theBremParticle.getSuccess() != 0)
734  brem.addPoint(PFTrajectoryPoint(-1,
735  PFTrajectoryPoint::HOLayer,
736  math::XYZPoint(theBremParticle.particle().vertex()),
737  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
738  else {
739  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
740  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
741  PFTrajectoryPoint dummyHOLayer;
742  brem.addPoint(dummyHOLayer);
743  }
744  }
745 
746  pftrack.addBrem(brem);
747  iTrajPos++;
748 
749  // ***************************** INTERMIDIATE State *************************************
750  //From the new Wolfgang code
751 
752  // To think if the cout should be removed.
753  if (track.gsfExtra()->tangentsSize() == 0)
754  LogError("PFTrackTransformer")
755  << "BE CAREFUL: Gsf Tangents not stored in the event. You need to re-reco the particle-flow with "
756  "RecoToDisplay_cfg.py and not RecoToDisplay_NoTracking_cfg.py ";
757 
758  vector<GsfTangent> gsftang = track.gsfExtra()->tangents();
759  for (unsigned int iTang = 0; iTang < track.gsfExtra()->tangentsSize(); iTang++) {
760  dp_tang = gsftang[iTang].deltaP().value();
761  sdp_tang = gsftang[iTang].deltaP().error();
762 
763  //check that the vertex of the brem is in the tracker volume
764  if ((sqrt(gsftang[iTang].position().x() * gsftang[iTang].position().x() +
765  gsftang[iTang].position().y() * gsftang[iTang].position().y()) > 110) ||
766  (fabs(gsftang[iTang].position().z()) > 280))
767  continue;
768 
769  iTrajPoint = iTrajPos + 2;
770  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
771 
772  GlobalVector p_tang =
773  GlobalVector(gsftang[iTang].momentum().x(), gsftang[iTang].momentum().y(), gsftang[iTang].momentum().z());
774 
775  // ###### track tj points
776  pftrack.addPoint(PFTrajectoryPoint(
777  iid,
778  -1,
779  math::XYZPoint(gsftang[iTang].position().x(), gsftang[iTang].position().y(), gsftang[iTang].position().z()),
780  math::XYZTLorentzVector(p_tang.x(), p_tang.y(), p_tang.z(), p_tang.mag())));
781 
782  //rescale
783  GlobalVector p_gamma = p_tang * (fabs(dp_tang) / p_tang.mag());
784 
785  // GlobalVector
786 
787  double e_gamma = fabs(dp_tang); // DP = pout-pin so could be negative
788  theBremParticle = BaseParticlePropagator(
789  RawParticle(
790  XYZTLorentzVector(p_gamma.x(), p_gamma.y(), p_gamma.z(), e_gamma),
791  XYZTLorentzVector(
792  gsftang[iTang].position().x(), gsftang[iTang].position().y(), gsftang[iTang].position().z(), 0.),
793  gamma_charge),
794  0.,
795  0.,
796  B_.z());
797 
798  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
799  brem.addPoint(dummyClosest);
800 
801  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
802  brem.addPoint(dummyBeamPipe);
803 
804  isBelowPS = false;
805  theBremParticle.propagateToPreshowerLayer1(false);
806  if (theBremParticle.getSuccess() != 0)
807  brem.addPoint(PFTrajectoryPoint(-1,
808  PFTrajectoryPoint::PS1,
809  math::XYZPoint(theBremParticle.particle().vertex()),
810  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
811  else {
812  PFTrajectoryPoint dummyPS1;
813  brem.addPoint(dummyPS1);
814  }
815 
816  // Brem Entrance PS Layer 2
817 
818  theBremParticle.propagateToPreshowerLayer2(false);
819  if (theBremParticle.getSuccess() != 0) {
820  brem.addPoint(PFTrajectoryPoint(-1,
821  PFTrajectoryPoint::PS2,
822  math::XYZPoint(theBremParticle.particle().vertex()),
823  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
824  isBelowPS = true;
825  } else {
826  PFTrajectoryPoint dummyPS2;
827  brem.addPoint(dummyPS2);
828  }
829 
830  //ECAL entrance
831  theBremParticle.propagateToEcalEntrance(false);
832  if (theBremParticle.getSuccess() != 0) {
833  brem.addPoint(PFTrajectoryPoint(-1,
834  PFTrajectoryPoint::ECALEntrance,
835  math::XYZPoint(theBremParticle.particle().vertex()),
836  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
837 
838  double ecalShowerDepth =
839  PFCluster::getDepthCorrection(theBremParticle.particle().momentum().E(), isBelowPS, false);
840 
841  math::XYZPoint meanShower =
842  math::XYZPoint(theBremParticle.particle().vertex()) +
843  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
844 
845  brem.addPoint(PFTrajectoryPoint(-1,
846  PFTrajectoryPoint::ECALShowerMax,
847  meanShower,
848  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
849  } else {
850  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
851  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
852  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
853  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
854  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
855  << " " << theBremParticle.particle().cos2ThetaV();
856  PFTrajectoryPoint dummyECAL;
857  brem.addPoint(dummyECAL);
858  PFTrajectoryPoint dummyMaxSh;
859  brem.addPoint(dummyMaxSh);
860  }
861 
862  //HCAL entrance
863  theBremParticle.propagateToHcalEntrance(false);
864  if (theBremParticle.getSuccess() != 0)
865  brem.addPoint(PFTrajectoryPoint(-1,
866  PFTrajectoryPoint::HCALEntrance,
867  math::XYZPoint(theBremParticle.particle().vertex()),
868  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
869  else {
870  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
871  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
872  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
873  PFTrajectoryPoint dummyHCALentrance;
874  brem.addPoint(dummyHCALentrance);
875  }
876 
877  //HCAL exit
878  theBremParticle.propagateToHcalExit(false);
879  if (theBremParticle.getSuccess() != 0)
880  brem.addPoint(PFTrajectoryPoint(-1,
881  PFTrajectoryPoint::HCALExit,
882  math::XYZPoint(theBremParticle.particle().vertex()),
883  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
884  else {
885  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
886  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
887  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
888  PFTrajectoryPoint dummyHCALexit;
889  brem.addPoint(dummyHCALexit);
890  }
891 
892  //HO Layer0
893  if (abs(theBremParticle.particle().vertex().z()) < 550) {
894  theBremParticle.setMagneticField(0);
895  theBremParticle.propagateToHOLayer(false);
896  if (theBremParticle.getSuccess() != 0)
897  brem.addPoint(PFTrajectoryPoint(-1,
898  PFTrajectoryPoint::HOLayer,
899  math::XYZPoint(theBremParticle.particle().vertex()),
900  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
901  else {
902  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
903  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
904  PFTrajectoryPoint dummyHOLayer;
905  brem.addPoint(dummyHOLayer);
906  }
907  }
908 
909  pftrack.addBrem(brem);
910  iTrajPos++;
911  }
912 
913  // ***************************** OUTER State *************************************
914 
915  if (outTSOS.isValid()) {
916  GlobalVector OutMom;
917  GlobalPoint OutPos;
918 
919  // DANIELE ????? if the out is not valid maybe take the last tangent?
920  // From Wolfgang. It should be always valid
921 
922  multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, OutMom);
923  multiTrajectoryStateMode::positionFromModeCartesian(outTSOS, OutPos);
924 
925  // last tjpoint
926  pftrack.addPoint(PFTrajectoryPoint(iid,
927  -1,
928  math::XYZPoint(OutPos.x(), OutPos.y(), OutPos.z()),
929  math::XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), OutMom.mag())));
930 
931  float ptot_out = sqrt((OutMom.x() * OutMom.x()) + (OutMom.y() * OutMom.y()) + (OutMom.z() * OutMom.z()));
932  float pTtot_out = sqrt((OutMom.x() * OutMom.x()) + (OutMom.y() * OutMom.y()));
933  float pfenergy_out = ptot_out;
934  BaseParticlePropagator theOutParticle =
935  BaseParticlePropagator(RawParticle(XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), pfenergy_out),
936  XYZTLorentzVector(OutPos.x(), OutPos.y(), OutPos.z(), 0.),
937  track.charge()),
938  0.,
939  0.,
940  B_.z());
941  isBelowPS = false;
942  theOutParticle.propagateToPreshowerLayer1(false);
943  if (theOutParticle.getSuccess() != 0)
944  pftrack.addPoint(PFTrajectoryPoint(-1,
945  PFTrajectoryPoint::PS1,
946  math::XYZPoint(theOutParticle.particle().vertex()),
947  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
948  else {
949  PFTrajectoryPoint dummyPS1;
950  pftrack.addPoint(dummyPS1);
951  }
952 
953  theOutParticle.propagateToPreshowerLayer2(false);
954  if (theOutParticle.getSuccess() != 0) {
955  pftrack.addPoint(PFTrajectoryPoint(-1,
956  PFTrajectoryPoint::PS2,
957  math::XYZPoint(theOutParticle.particle().vertex()),
958  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
959  isBelowPS = true;
960  } else {
961  PFTrajectoryPoint dummyPS2;
962  pftrack.addPoint(dummyPS2);
963  }
964 
965  //ECAL entrance
966  theOutParticle.propagateToEcalEntrance(false);
967  if (theOutParticle.getSuccess() != 0) {
968  pftrack.addPoint(PFTrajectoryPoint(-1,
969  PFTrajectoryPoint::ECALEntrance,
970  math::XYZPoint(theOutParticle.particle().vertex()),
971  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
972  double EDepthCorr = 0.01;
973  double ecalShowerDepth = PFCluster::getDepthCorrection(EDepthCorr, isBelowPS, false);
974 
975  math::XYZPoint meanShower =
976  math::XYZPoint(theOutParticle.particle().vertex()) +
977  math::XYZTLorentzVector(theOutParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
978 
979  pftrack.addPoint(PFTrajectoryPoint(-1,
980  PFTrajectoryPoint::ECALShowerMax,
981  meanShower,
982  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
983  } else {
984  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
985  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE ECAL HAS FAILED\n"
986  << "theOutParticle.particle() pt,eta: " << theOutParticle.particle().pt()
987  << " " << theOutParticle.particle().eta();
988 
989  PFTrajectoryPoint dummyECAL;
990  pftrack.addPoint(dummyECAL);
991  PFTrajectoryPoint dummyMaxSh;
992  pftrack.addPoint(dummyMaxSh);
993  }
994 
995  //HCAL entrance
996  theOutParticle.propagateToHcalEntrance(false);
997  if (theOutParticle.getSuccess() != 0)
998  pftrack.addPoint(PFTrajectoryPoint(-1,
999  PFTrajectoryPoint::HCALEntrance,
1000  math::XYZPoint(theOutParticle.particle().vertex()),
1001  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1002  else {
1003  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1004  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
1005  PFTrajectoryPoint dummyHCALentrance;
1006  pftrack.addPoint(dummyHCALentrance);
1007  }
1008  //HCAL exit
1009  theOutParticle.propagateToHcalExit(false);
1010  if (theOutParticle.getSuccess() != 0)
1011  pftrack.addPoint(PFTrajectoryPoint(-1,
1012  PFTrajectoryPoint::HCALExit,
1013  math::XYZPoint(theOutParticle.particle().vertex()),
1014  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1015  else {
1016  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1017  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
1018  PFTrajectoryPoint dummyHCALexit;
1019  pftrack.addPoint(dummyHCALexit);
1020  }
1021 
1022  //HO Layer0
1023  if (abs(theOutParticle.particle().vertex().z()) < 550) {
1024  theOutParticle.setMagneticField(0);
1025  theOutParticle.propagateToHOLayer(false);
1026  if (theOutParticle.getSuccess() != 0)
1027  pftrack.addPoint(PFTrajectoryPoint(-1,
1028  PFTrajectoryPoint::HOLayer,
1029  math::XYZPoint(theOutParticle.particle().vertex()),
1030  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1031  else {
1032  if (pTtot_out > 5. && abs(theOutParticle.particle().Z()) < 700.25)
1033  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HO HAS FAILED";
1034  PFTrajectoryPoint dummyHOLayer;
1035  pftrack.addPoint(dummyHOLayer);
1036  }
1037  }
1038  //######## Photon at the OUTER State ##########
1039 
1040  dp_tang = OutMom.mag();
1041  // for the moment same inner error just for semplicity
1042  sdp_tang = track.ptModeError() * (track.pMode() / track.ptMode());
1043  iTrajPoint = iTrajPos + 2;
1044  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
1045 
1046  theBremParticle = BaseParticlePropagator(RawParticle(XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), dp_tang),
1047  XYZTLorentzVector(OutPos.x(), OutPos.y(), OutPos.z(), 0.),
1048  gamma_charge),
1049  0.,
1050  0.,
1051  B_.z());
1052 
1053  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
1054  brem.addPoint(dummyClosest);
1055 
1056  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
1057  brem.addPoint(dummyBeamPipe);
1058 
1059  isBelowPS = false;
1060  theBremParticle.propagateToPreshowerLayer1(false);
1061  if (theBremParticle.getSuccess() != 0)
1062  brem.addPoint(PFTrajectoryPoint(-1,
1063  PFTrajectoryPoint::PS1,
1064  math::XYZPoint(theBremParticle.particle().vertex()),
1065  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1066  else {
1067  PFTrajectoryPoint dummyPS1;
1068  brem.addPoint(dummyPS1);
1069  }
1070 
1071  // Brem Entrance PS Layer 2
1072 
1073  theBremParticle.propagateToPreshowerLayer2(false);
1074  if (theBremParticle.getSuccess() != 0) {
1075  brem.addPoint(PFTrajectoryPoint(-1,
1076  PFTrajectoryPoint::PS2,
1077  math::XYZPoint(theBremParticle.particle().vertex()),
1078  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1079  isBelowPS = true;
1080  } else {
1081  PFTrajectoryPoint dummyPS2;
1082  brem.addPoint(dummyPS2);
1083  }
1084 
1085  //ECAL entrance
1086  theBremParticle.propagateToEcalEntrance(false);
1087  if (theBremParticle.getSuccess() != 0) {
1088  brem.addPoint(PFTrajectoryPoint(-1,
1089  PFTrajectoryPoint::ECALEntrance,
1090  math::XYZPoint(theBremParticle.particle().vertex()),
1091  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1092  double ecalShowerDepth =
1093  PFCluster::getDepthCorrection(theBremParticle.particle().momentum().E(), isBelowPS, false);
1094 
1095  math::XYZPoint meanShower =
1096  math::XYZPoint(theBremParticle.particle().vertex()) +
1097  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
1098 
1099  brem.addPoint(PFTrajectoryPoint(-1,
1100  PFTrajectoryPoint::ECALShowerMax,
1101  meanShower,
1102  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1103  } else {
1104  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1105  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1106  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
1107  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
1108  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
1109  << " " << theBremParticle.particle().cos2ThetaV();
1110  PFTrajectoryPoint dummyECAL;
1111  brem.addPoint(dummyECAL);
1112  PFTrajectoryPoint dummyMaxSh;
1113  brem.addPoint(dummyMaxSh);
1114  }
1115 
1116  //HCAL entrance
1117  theBremParticle.propagateToHcalEntrance(false);
1118  if (theBremParticle.getSuccess() != 0)
1119  brem.addPoint(PFTrajectoryPoint(-1,
1120  PFTrajectoryPoint::HCALEntrance,
1121  math::XYZPoint(theBremParticle.particle().vertex()),
1122  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1123  else {
1124  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1125  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1126  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
1127  PFTrajectoryPoint dummyHCALentrance;
1128  brem.addPoint(dummyHCALentrance);
1129  }
1130  //HCAL exit
1131  theBremParticle.propagateToHcalExit(false);
1132  if (theBremParticle.getSuccess() != 0)
1133  brem.addPoint(PFTrajectoryPoint(-1,
1134  PFTrajectoryPoint::HCALExit,
1135  math::XYZPoint(theBremParticle.particle().vertex()),
1136  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1137  else {
1138  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1139  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1140  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
1141  PFTrajectoryPoint dummyHCALexit;
1142  brem.addPoint(dummyHCALexit);
1143  }
1144 
1145  //HO Layer0
1146  if (abs(theBremParticle.particle().vertex().z()) < 550) {
1147  theBremParticle.setMagneticField(0);
1148  theBremParticle.propagateToHOLayer(false);
1149  if (theBremParticle.getSuccess() != 0)
1150  brem.addPoint(PFTrajectoryPoint(-1,
1151  PFTrajectoryPoint::HOLayer,
1152  math::XYZPoint(theBremParticle.particle().vertex()),
1153  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1154  else {
1155  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
1156  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
1157  PFTrajectoryPoint dummyHOLayer;
1158  brem.addPoint(dummyHOLayer);
1159  }
1160  }
1161  pftrack.addBrem(brem);
1162  iTrajPos++;
1163  }
1164 
1165  return true;
1166 }
Vector3DBase
Definition: Vector3DBase.h:8
PFTrackTransformer::onlyprop_
bool onlyprop_
Definition: PFTrackTransformer.h:55
PFGsfHelper.h
BaseParticlePropagator::propagateToHcalExit
bool propagateToHcalExit(bool first=true)
Definition: BaseParticlePropagator.cc:528
DDAxes::y
TrajectoryStateOnSurface.h
RawParticle
Definition: RawParticle.h:37
HLT_FULL_cff.track
track
Definition: HLT_FULL_cff.py:11713
BaseParticlePropagator::propagateToEcalEntrance
bool propagateToEcalEntrance(bool first=true)
Definition: BaseParticlePropagator.cc:450
MessageLogger.h
BaseParticlePropagator::particle
RawParticle const & particle() const
The particle being propagated.
Definition: BaseParticlePropagator.h:164
MultiGaussianStateTransform.h
reco::PFTrack::addPoint
void addPoint(const reco::PFTrajectoryPoint &trajPt)
Definition: PFTrack.cc:28
ESHandle.h
RawParticle::momentum
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:133
RawParticle::vertex
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
edm
HLT enums.
Definition: AlignableModifier.h:19
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
PT
#define PT
Definition: PFJetBenchmark.cc:41
BaseParticlePropagator.h
PFGeometry::outerRadius
float outerRadius(PFGeometry::Layers_t layer) const
return outer radius of a given layer
Definition: PFGeometry.h:59
BaseParticlePropagator::getSuccess
int getSuccess() const
Has propagation been performed and was barrel or endcap reached ?
Definition: BaseParticlePropagator.h:296
multiTrajectoryStateMode::momentumFromModeCartesian
bool momentumFromModeCartesian(TrajectoryStateOnSurface const &tsos, GlobalVector &momentum)
Definition: MultiTrajectoryStateMode.cc:16
DDAxes::x
XYZTLorentzVector
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
GaussianSumUtilities1D.h
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
PFTrackTransformer::cos2ThetaV_Endcap_HiEnd_
static constexpr double cos2ThetaV_Endcap_HiEnd_
Definition: PFTrackTransformer.h:61
findQualityFiles.v
v
Definition: findQualityFiles.py:179
reco::PFBrem
Definition: PFBrem.h:10
PFTrackTransformer.h
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
BaseParticlePropagator
Definition: BaseParticlePropagator.h:82
reco::GsfPFRecTrack
Definition: GsfPFRecTrack.h:24
MultiTrajectoryStateTransform::outerStateOnSurface
TrajectoryStateOnSurface outerStateOnSurface(const reco::GsfTrack &tk) const
Definition: MultiTrajectoryStateTransform.cc:19
MultiGaussianState1D.h
PFTrackTransformer::PFTrackTransformer
PFTrackTransformer(const math::XYZVector &)
Definition: PFTrackTransformer.cc:32
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
RawParticle::Z
double Z() const
z of vertex
Definition: RawParticle.h:288
BaseParticlePropagator::propagateToHcalEntrance
bool propagateToHcalEntrance(bool first=true)
Definition: BaseParticlePropagator.cc:476
reco::GsfTrack
Definition: GsfTrack.h:12
TrajectoryStateOnSurface
Definition: TrajectoryStateOnSurface.h:16
PFTrackTransformer::addPoints
bool addPoints(reco::PFRecTrack &pftrack, const reco::Track &track, const Trajectory &traj, bool msgwarning=true) const
Add points to a PFTrack. return false if a TSOS is invalid.
Definition: PFTrackTransformer.cc:40
Track.h
BaseParticlePropagator::propagateToPreshowerLayer2
bool propagateToPreshowerLayer2(bool first=true)
Definition: BaseParticlePropagator.cc:431
RawParticle::eta
double eta() const
Definition: RawParticle.h:279
RawParticle::pt
double pt() const
transverse momentum
Definition: RawParticle.h:309
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
DDAxes::z
reco::Track
Definition: Track.h:27
PFCluster.h
PFGeometry::outerZ
float outerZ(PFGeometry::Layers_t layer) const
return outer position along z axis of a given layer
Definition: PFGeometry.h:65
BaseParticlePropagator::propagateToHOLayer
bool propagateToHOLayer(bool first=true)
Definition: BaseParticlePropagator.cc:545
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
Point3DBase< float, GlobalTag >
reco::PFRecTrack::KF_ELCAND
Definition: PFRecTrack.h:27
reco::PFRecTrack::algoType
unsigned int algoType() const
Definition: PFRecTrack.h:39
PFGsfHelper::sigmafittedDP
double sigmafittedDP() const
Definition: PFGsfHelper.cc:131
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:223
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
PFTrackTransformer::~PFTrackTransformer
~PFTrackTransformer()
Definition: PFTrackTransformer.cc:38
BaseParticlePropagator::propagate
bool propagate()
Definition: BaseParticlePropagator.cc:34
RawParticle::cos2ThetaV
double cos2ThetaV() const
Definition: RawParticle.h:281
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
PFGsfHelper::fittedDP
double fittedDP() const
Definition: PFGsfHelper.cc:130
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
multiTrajectoryStateMode::positionFromModeCartesian
bool positionFromModeCartesian(TrajectoryStateOnSurface const &tsos, GlobalPoint &position)
Definition: MultiTrajectoryStateMode.cc:63
BaseParticlePropagator::propagateToPreshowerLayer1
bool propagateToPreshowerLayer1(bool first=true)
Definition: BaseParticlePropagator.cc:412
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
BaseParticlePropagator::setMagneticField
void setMagneticField(double b)
Set the magnetic field.
Definition: BaseParticlePropagator.h:299
Trajectory::measurements
DataContainer const & measurements() const
Definition: Trajectory.h:178
BaseParticlePropagator::setPropagationConditions
void setPropagationConditions(double r, double z, bool firstLoop=true)
Set the propagation characteristics (rCyl, zCyl and first loop only)
Definition: BaseParticlePropagator.cc:318
MultiTrajectoryStateTransform
Definition: MultiTrajectoryStateTransform.h:18
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
reco::PFRecTrack::eta
auto eta() const
Definition: PFRecTrack.h:54
Trajectory.h
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
std
Definition: JetResolutionObject.h:76
PFTrackTransformer::addPointsAndBrems
bool addPointsAndBrems(reco::GsfPFRecTrack &pftrack, const reco::Track &track, const Trajectory &traj, const bool &GetMode) const
Definition: PFTrackTransformer.cc:234
reco::GsfPFRecTrack::addBrem
void addBrem(const reco::PFBrem &brem)
add a Bremsstrahlung photon
Definition: GsfPFRecTrack.cc:16
reco::PFRecTrack
reconstructed track used as an input to particle flow
Definition: PFRecTrack.h:20
Trajectory
Definition: Trajectory.h:38
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
PFGeometry::BeamPipe
Definition: PFGeometry.h:26
reco::PFTrajectoryPoint
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
Definition: PFTrajectoryPoint.h:26
PFGsfHelper::computeP
GlobalVector computeP(bool ComputeMode) const
Definition: PFGsfHelper.cc:122
MultiTrajectoryStateTransform::innerStateOnSurface
TrajectoryStateOnSurface innerStateOnSurface(const reco::GsfTrack &tk) const
Definition: MultiTrajectoryStateTransform.cc:23
PFRecTrackFwd.h
PFGsfHelper
Definition: PFGsfHelper.h:29
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
alongMomentum
Definition: PropagationDirection.h:4
TrajectoryStateOnSurface::isValid
bool isValid() const
Definition: TrajectoryStateOnSurface.h:54
PFTrackTransformer::pfGeometry_
PFGeometry pfGeometry_
Definition: PFTrackTransformer.h:60
PFTrackTransformer::B_
math::XYZVector B_
B field.
Definition: PFTrackTransformer.h:59
Generated for CMSSW Reference Manual by doxygen 1.8.16