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