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  brem.calculatePositionREP();
549  pftrack.addBrem(brem);
550  iTrajPos++;
551  }
552  return true;
553 }
554 
555 bool PFTrackTransformer::addPointsAndBrems(reco::GsfPFRecTrack& pftrack,
556  const reco::GsfTrack& track,
557  const MultiTrajectoryStateTransform& mtjstate) const {
558  // float PT= track.pt();
559  unsigned int iTrajPos = 0;
560  unsigned int iid = 0; // not anymore saved
561 
562  // ***************************** INNER State *************************************
563  TrajectoryStateOnSurface inTSOS = mtjstate.innerStateOnSurface((track));
564  TrajectoryStateOnSurface outTSOS = mtjstate.outerStateOnSurface((track));
565 
566  if (!inTSOS.isValid() || !outTSOS.isValid()) {
567  if (!inTSOS.isValid())
568  LogWarning("PFTrackTransformer") << " INNER TSOS NOT VALID ";
569  if (!outTSOS.isValid())
570  LogWarning("PFTrackTransformer") << " OUTER TSOS NOT VALID ";
571  return false;
572  }
573 
574  GlobalVector InMom;
575  GlobalPoint InPos;
576  if (inTSOS.isValid()) {
577  multiTrajectoryStateMode::momentumFromModeCartesian(inTSOS, InMom);
578  multiTrajectoryStateMode::positionFromModeCartesian(inTSOS, InPos);
579  } else {
580  InMom = GlobalVector(track.pxMode(), track.pyMode(), track.pzMode());
581  InPos = GlobalPoint(0., 0., 0.);
582  }
583 
584  // float pfmass= (pftrack.algoType()==reco::PFRecTrack::KF_ELCAND) ? 0.0005 : 0.139;
585  float ptot = sqrt((InMom.x() * InMom.x()) + (InMom.y() * InMom.y()) + (InMom.z() * InMom.z()));
586  float pfenergy = ptot;
587 
588  math::XYZTLorentzVector momClosest = math::XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), ptot);
589  const math::XYZPoint& posClosest = track.vertex();
590  pftrack.addPoint(PFTrajectoryPoint(-1, PFTrajectoryPoint::ClosestApproach, posClosest, momClosest));
591 
592  BaseParticlePropagator theInnerParticle = BaseParticlePropagator(
593  RawParticle(XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), pfenergy),
594  XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.),
595  track.charge()), //DANIELE Same thing v.x(),v.y(),v.()?
596  0.,
597  0.,
598  B_.z());
599  //BEAMPIPE
600  theInnerParticle.setPropagationConditions(
601  pfGeometry_.outerRadius(PFGeometry::BeamPipe), pfGeometry_.outerZ(PFGeometry::BeamPipe), false);
602  theInnerParticle.propagate();
603  if (theInnerParticle.getSuccess() != 0)
604  pftrack.addPoint(PFTrajectoryPoint(-1,
605  PFTrajectoryPoint::BeamPipeOrEndVertex,
606  math::XYZPoint(theInnerParticle.particle().vertex()),
607  math::XYZTLorentzVector(theInnerParticle.particle().momentum())));
608  else {
609  PFTrajectoryPoint dummyBeam;
610  pftrack.addPoint(dummyBeam);
611  }
612 
613  // first tjpoint
614  pftrack.addPoint(PFTrajectoryPoint(iid,
615  -1,
616  math::XYZPoint(InPos.x(), InPos.y(), InPos.z()),
617  math::XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), InMom.mag())));
618 
619  //######### Photon at INNER State ##########
620 
621  unsigned int iTrajPoint = iTrajPos + 2;
622  double dp_tang = ptot;
623  double sdp_tang = track.ptModeError() * (track.pMode() / track.ptMode());
624  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
625  constexpr int gamma_charge = 0;
626  BaseParticlePropagator theBremParticle =
627  BaseParticlePropagator(RawParticle(XYZTLorentzVector(InMom.x(), InMom.y(), InMom.z(), dp_tang),
628  XYZTLorentzVector(InPos.x(), InPos.y(), InPos.z(), 0.),
629  gamma_charge),
630  0.,
631  0.,
632  B_.z());
633  // add TrajectoryPoint for Brem, PS, ECAL, ECALShowMax, HCAL
634  // Brem Entrance PS Layer1
635  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
636  brem.addPoint(dummyClosest);
637 
638  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
639  brem.addPoint(dummyBeamPipe);
640 
641  bool isBelowPS = false;
642  theBremParticle.propagateToPreshowerLayer1(false);
643  if (theBremParticle.getSuccess() != 0)
644  brem.addPoint(PFTrajectoryPoint(-1,
645  PFTrajectoryPoint::PS1,
646  math::XYZPoint(theBremParticle.particle().vertex()),
647  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
648  else {
649  PFTrajectoryPoint dummyPS1;
650  brem.addPoint(dummyPS1);
651  }
652 
653  // Brem Entrance PS Layer 2
654 
655  theBremParticle.propagateToPreshowerLayer2(false);
656  if (theBremParticle.getSuccess() != 0) {
657  brem.addPoint(PFTrajectoryPoint(-1,
658  PFTrajectoryPoint::PS2,
659  math::XYZPoint(theBremParticle.particle().vertex()),
660  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
661  isBelowPS = true;
662  } else {
663  PFTrajectoryPoint dummyPS2;
664  brem.addPoint(dummyPS2);
665  }
666 
667  //ECAL entrance
668  theBremParticle.propagateToEcalEntrance(false);
669  if (theBremParticle.getSuccess() != 0) {
670  brem.addPoint(PFTrajectoryPoint(-1,
671  PFTrajectoryPoint::ECALEntrance,
672  math::XYZPoint(theBremParticle.particle().vertex()),
673  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
674 
675  // for the first brem give a low default DP of 100 MeV.
676  double EDepthCorr = 0.01;
677  double ecalShowerDepth = PFCluster::getDepthCorrection(EDepthCorr, isBelowPS, false);
678 
679  math::XYZPoint meanShower =
680  math::XYZPoint(theBremParticle.particle().vertex()) +
681  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
682 
683  brem.addPoint(PFTrajectoryPoint(-1,
684  PFTrajectoryPoint::ECALShowerMax,
685  meanShower,
686  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
687  } else {
688  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
689  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
690  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
691  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
692  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
693  << " " << theBremParticle.particle().cos2ThetaV();
694  PFTrajectoryPoint dummyECAL;
695  brem.addPoint(dummyECAL);
696  PFTrajectoryPoint dummyMaxSh;
697  brem.addPoint(dummyMaxSh);
698  }
699 
700  //HCAL entrance
701  theBremParticle.propagateToHcalEntrance(false);
702  if (theBremParticle.getSuccess() != 0)
703  brem.addPoint(PFTrajectoryPoint(-1,
704  PFTrajectoryPoint::HCALEntrance,
705  math::XYZPoint(theBremParticle.particle().vertex()),
706  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
707  else {
708  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
709  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
710  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
711  PFTrajectoryPoint dummyHCALentrance;
712  brem.addPoint(dummyHCALentrance);
713  }
714 
715  //HCAL exit
716  theBremParticle.propagateToHcalExit(false);
717  if (theBremParticle.getSuccess() != 0)
718  brem.addPoint(PFTrajectoryPoint(-1,
719  PFTrajectoryPoint::HCALExit,
720  math::XYZPoint(theBremParticle.particle().vertex()),
721  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
722  else {
723  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
724  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
725  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
726  PFTrajectoryPoint dummyHCALexit;
727  brem.addPoint(dummyHCALexit);
728  }
729 
730  //HO Layer0
731  if (abs(theBremParticle.particle().vertex().z()) < 550) {
732  theBremParticle.setMagneticField(0);
733  theBremParticle.propagateToHOLayer(false);
734  if (theBremParticle.getSuccess() != 0)
735  brem.addPoint(PFTrajectoryPoint(-1,
736  PFTrajectoryPoint::HOLayer,
737  math::XYZPoint(theBremParticle.particle().vertex()),
738  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
739  else {
740  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
741  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
742  PFTrajectoryPoint dummyHOLayer;
743  brem.addPoint(dummyHOLayer);
744  }
745  }
746 
747  brem.calculatePositionREP();
748  pftrack.addBrem(brem);
749  iTrajPos++;
750 
751  // ***************************** INTERMIDIATE State *************************************
752  //From the new Wolfgang code
753 
754  // To think if the cout should be removed.
755  if (track.gsfExtra()->tangentsSize() == 0)
756  LogError("PFTrackTransformer")
757  << "BE CAREFUL: Gsf Tangents not stored in the event. You need to re-reco the particle-flow with "
758  "RecoToDisplay_cfg.py and not RecoToDisplay_NoTracking_cfg.py ";
759 
760  vector<GsfTangent> gsftang = track.gsfExtra()->tangents();
761  for (unsigned int iTang = 0; iTang < track.gsfExtra()->tangentsSize(); iTang++) {
762  dp_tang = gsftang[iTang].deltaP().value();
763  sdp_tang = gsftang[iTang].deltaP().error();
764 
765  //check that the vertex of the brem is in the tracker volume
766  if ((sqrt(gsftang[iTang].position().x() * gsftang[iTang].position().x() +
767  gsftang[iTang].position().y() * gsftang[iTang].position().y()) > 110) ||
768  (fabs(gsftang[iTang].position().z()) > 280))
769  continue;
770 
771  iTrajPoint = iTrajPos + 2;
772  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
773 
774  GlobalVector p_tang =
775  GlobalVector(gsftang[iTang].momentum().x(), gsftang[iTang].momentum().y(), gsftang[iTang].momentum().z());
776 
777  // ###### track tj points
778  pftrack.addPoint(PFTrajectoryPoint(
779  iid,
780  -1,
781  math::XYZPoint(gsftang[iTang].position().x(), gsftang[iTang].position().y(), gsftang[iTang].position().z()),
782  math::XYZTLorentzVector(p_tang.x(), p_tang.y(), p_tang.z(), p_tang.mag())));
783 
784  //rescale
785  GlobalVector p_gamma = p_tang * (fabs(dp_tang) / p_tang.mag());
786 
787  // GlobalVector
788 
789  double e_gamma = fabs(dp_tang); // DP = pout-pin so could be negative
790  theBremParticle = BaseParticlePropagator(
791  RawParticle(
792  XYZTLorentzVector(p_gamma.x(), p_gamma.y(), p_gamma.z(), e_gamma),
793  XYZTLorentzVector(
794  gsftang[iTang].position().x(), gsftang[iTang].position().y(), gsftang[iTang].position().z(), 0.),
795  gamma_charge),
796  0.,
797  0.,
798  B_.z());
799 
800  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
801  brem.addPoint(dummyClosest);
802 
803  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
804  brem.addPoint(dummyBeamPipe);
805 
806  isBelowPS = false;
807  theBremParticle.propagateToPreshowerLayer1(false);
808  if (theBremParticle.getSuccess() != 0)
809  brem.addPoint(PFTrajectoryPoint(-1,
810  PFTrajectoryPoint::PS1,
811  math::XYZPoint(theBremParticle.particle().vertex()),
812  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
813  else {
814  PFTrajectoryPoint dummyPS1;
815  brem.addPoint(dummyPS1);
816  }
817 
818  // Brem Entrance PS Layer 2
819 
820  theBremParticle.propagateToPreshowerLayer2(false);
821  if (theBremParticle.getSuccess() != 0) {
822  brem.addPoint(PFTrajectoryPoint(-1,
823  PFTrajectoryPoint::PS2,
824  math::XYZPoint(theBremParticle.particle().vertex()),
825  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
826  isBelowPS = true;
827  } else {
828  PFTrajectoryPoint dummyPS2;
829  brem.addPoint(dummyPS2);
830  }
831 
832  //ECAL entrance
833  theBremParticle.propagateToEcalEntrance(false);
834  if (theBremParticle.getSuccess() != 0) {
835  brem.addPoint(PFTrajectoryPoint(-1,
836  PFTrajectoryPoint::ECALEntrance,
837  math::XYZPoint(theBremParticle.particle().vertex()),
838  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
839 
840  double ecalShowerDepth =
841  PFCluster::getDepthCorrection(theBremParticle.particle().momentum().E(), isBelowPS, false);
842 
843  math::XYZPoint meanShower =
844  math::XYZPoint(theBremParticle.particle().vertex()) +
845  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
846 
847  brem.addPoint(PFTrajectoryPoint(-1,
848  PFTrajectoryPoint::ECALShowerMax,
849  meanShower,
850  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
851  } else {
852  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
853  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
854  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
855  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
856  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
857  << " " << theBremParticle.particle().cos2ThetaV();
858  PFTrajectoryPoint dummyECAL;
859  brem.addPoint(dummyECAL);
860  PFTrajectoryPoint dummyMaxSh;
861  brem.addPoint(dummyMaxSh);
862  }
863 
864  //HCAL entrance
865  theBremParticle.propagateToHcalEntrance(false);
866  if (theBremParticle.getSuccess() != 0)
867  brem.addPoint(PFTrajectoryPoint(-1,
868  PFTrajectoryPoint::HCALEntrance,
869  math::XYZPoint(theBremParticle.particle().vertex()),
870  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
871  else {
872  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
873  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
874  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
875  PFTrajectoryPoint dummyHCALentrance;
876  brem.addPoint(dummyHCALentrance);
877  }
878 
879  //HCAL exit
880  theBremParticle.propagateToHcalExit(false);
881  if (theBremParticle.getSuccess() != 0)
882  brem.addPoint(PFTrajectoryPoint(-1,
883  PFTrajectoryPoint::HCALExit,
884  math::XYZPoint(theBremParticle.particle().vertex()),
885  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
886  else {
887  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
888  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
889  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
890  PFTrajectoryPoint dummyHCALexit;
891  brem.addPoint(dummyHCALexit);
892  }
893 
894  //HO Layer0
895  if (abs(theBremParticle.particle().vertex().z()) < 550) {
896  theBremParticle.setMagneticField(0);
897  theBremParticle.propagateToHOLayer(false);
898  if (theBremParticle.getSuccess() != 0)
899  brem.addPoint(PFTrajectoryPoint(-1,
900  PFTrajectoryPoint::HOLayer,
901  math::XYZPoint(theBremParticle.particle().vertex()),
902  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
903  else {
904  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
905  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
906  PFTrajectoryPoint dummyHOLayer;
907  brem.addPoint(dummyHOLayer);
908  }
909  }
910 
911  brem.calculatePositionREP();
912  pftrack.addBrem(brem);
913  iTrajPos++;
914  }
915 
916  // ***************************** OUTER State *************************************
917 
918  if (outTSOS.isValid()) {
919  GlobalVector OutMom;
920  GlobalPoint OutPos;
921 
922  // DANIELE ????? if the out is not valid maybe take the last tangent?
923  // From Wolfgang. It should be always valid
924 
925  multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, OutMom);
926  multiTrajectoryStateMode::positionFromModeCartesian(outTSOS, OutPos);
927 
928  // last tjpoint
929  pftrack.addPoint(PFTrajectoryPoint(iid,
930  -1,
931  math::XYZPoint(OutPos.x(), OutPos.y(), OutPos.z()),
932  math::XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), OutMom.mag())));
933 
934  float ptot_out = sqrt((OutMom.x() * OutMom.x()) + (OutMom.y() * OutMom.y()) + (OutMom.z() * OutMom.z()));
935  float pTtot_out = sqrt((OutMom.x() * OutMom.x()) + (OutMom.y() * OutMom.y()));
936  float pfenergy_out = ptot_out;
937  BaseParticlePropagator theOutParticle =
938  BaseParticlePropagator(RawParticle(XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), pfenergy_out),
939  XYZTLorentzVector(OutPos.x(), OutPos.y(), OutPos.z(), 0.),
940  track.charge()),
941  0.,
942  0.,
943  B_.z());
944  isBelowPS = false;
945  theOutParticle.propagateToPreshowerLayer1(false);
946  if (theOutParticle.getSuccess() != 0)
947  pftrack.addPoint(PFTrajectoryPoint(-1,
948  PFTrajectoryPoint::PS1,
949  math::XYZPoint(theOutParticle.particle().vertex()),
950  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
951  else {
952  PFTrajectoryPoint dummyPS1;
953  pftrack.addPoint(dummyPS1);
954  }
955 
956  theOutParticle.propagateToPreshowerLayer2(false);
957  if (theOutParticle.getSuccess() != 0) {
958  pftrack.addPoint(PFTrajectoryPoint(-1,
959  PFTrajectoryPoint::PS2,
960  math::XYZPoint(theOutParticle.particle().vertex()),
961  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
962  isBelowPS = true;
963  } else {
964  PFTrajectoryPoint dummyPS2;
965  pftrack.addPoint(dummyPS2);
966  }
967 
968  //ECAL entrance
969  theOutParticle.propagateToEcalEntrance(false);
970  if (theOutParticle.getSuccess() != 0) {
971  pftrack.addPoint(PFTrajectoryPoint(-1,
972  PFTrajectoryPoint::ECALEntrance,
973  math::XYZPoint(theOutParticle.particle().vertex()),
974  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
975  double EDepthCorr = 0.01;
976  double ecalShowerDepth = PFCluster::getDepthCorrection(EDepthCorr, isBelowPS, false);
977 
978  math::XYZPoint meanShower =
979  math::XYZPoint(theOutParticle.particle().vertex()) +
980  math::XYZTLorentzVector(theOutParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
981 
982  pftrack.addPoint(PFTrajectoryPoint(-1,
983  PFTrajectoryPoint::ECALShowerMax,
984  meanShower,
985  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
986  } else {
987  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
988  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE ECAL HAS FAILED\n"
989  << "theOutParticle.particle() pt,eta: " << theOutParticle.particle().pt()
990  << " " << theOutParticle.particle().eta();
991 
992  PFTrajectoryPoint dummyECAL;
993  pftrack.addPoint(dummyECAL);
994  PFTrajectoryPoint dummyMaxSh;
995  pftrack.addPoint(dummyMaxSh);
996  }
997 
998  //HCAL entrance
999  theOutParticle.propagateToHcalEntrance(false);
1000  if (theOutParticle.getSuccess() != 0)
1001  pftrack.addPoint(PFTrajectoryPoint(-1,
1002  PFTrajectoryPoint::HCALEntrance,
1003  math::XYZPoint(theOutParticle.particle().vertex()),
1004  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1005  else {
1006  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1007  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
1008  PFTrajectoryPoint dummyHCALentrance;
1009  pftrack.addPoint(dummyHCALentrance);
1010  }
1011  //HCAL exit
1012  theOutParticle.propagateToHcalExit(false);
1013  if (theOutParticle.getSuccess() != 0)
1014  pftrack.addPoint(PFTrajectoryPoint(-1,
1015  PFTrajectoryPoint::HCALExit,
1016  math::XYZPoint(theOutParticle.particle().vertex()),
1017  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1018  else {
1019  if (pTtot_out > 5. && theOutParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1020  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
1021  PFTrajectoryPoint dummyHCALexit;
1022  pftrack.addPoint(dummyHCALexit);
1023  }
1024 
1025  //HO Layer0
1026  if (abs(theOutParticle.particle().vertex().z()) < 550) {
1027  theOutParticle.setMagneticField(0);
1028  theOutParticle.propagateToHOLayer(false);
1029  if (theOutParticle.getSuccess() != 0)
1030  pftrack.addPoint(PFTrajectoryPoint(-1,
1031  PFTrajectoryPoint::HOLayer,
1032  math::XYZPoint(theOutParticle.particle().vertex()),
1033  math::XYZTLorentzVector(theOutParticle.particle().momentum())));
1034  else {
1035  if (pTtot_out > 5. && abs(theOutParticle.particle().Z()) < 700.25)
1036  LogWarning("PFTrackTransformer") << "GSF TRACK " << pftrack << " PROPAGATION TO THE HO HAS FAILED";
1037  PFTrajectoryPoint dummyHOLayer;
1038  pftrack.addPoint(dummyHOLayer);
1039  }
1040  }
1041  //######## Photon at the OUTER State ##########
1042 
1043  dp_tang = OutMom.mag();
1044  // for the moment same inner error just for semplicity
1045  sdp_tang = track.ptModeError() * (track.pMode() / track.ptMode());
1046  iTrajPoint = iTrajPos + 2;
1047  PFBrem brem(dp_tang, sdp_tang, iTrajPoint);
1048 
1049  theBremParticle = BaseParticlePropagator(RawParticle(XYZTLorentzVector(OutMom.x(), OutMom.y(), OutMom.z(), dp_tang),
1050  XYZTLorentzVector(OutPos.x(), OutPos.y(), OutPos.z(), 0.),
1051  gamma_charge),
1052  0.,
1053  0.,
1054  B_.z());
1055 
1056  PFTrajectoryPoint dummyClosest; // Added just to have the right number order in PFTrack.cc
1057  brem.addPoint(dummyClosest);
1058 
1059  PFTrajectoryPoint dummyBeamPipe; // Added just to have the right number order in PFTrack.cc
1060  brem.addPoint(dummyBeamPipe);
1061 
1062  isBelowPS = false;
1063  theBremParticle.propagateToPreshowerLayer1(false);
1064  if (theBremParticle.getSuccess() != 0)
1065  brem.addPoint(PFTrajectoryPoint(-1,
1066  PFTrajectoryPoint::PS1,
1067  math::XYZPoint(theBremParticle.particle().vertex()),
1068  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1069  else {
1070  PFTrajectoryPoint dummyPS1;
1071  brem.addPoint(dummyPS1);
1072  }
1073 
1074  // Brem Entrance PS Layer 2
1075 
1076  theBremParticle.propagateToPreshowerLayer2(false);
1077  if (theBremParticle.getSuccess() != 0) {
1078  brem.addPoint(PFTrajectoryPoint(-1,
1079  PFTrajectoryPoint::PS2,
1080  math::XYZPoint(theBremParticle.particle().vertex()),
1081  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1082  isBelowPS = true;
1083  } else {
1084  PFTrajectoryPoint dummyPS2;
1085  brem.addPoint(dummyPS2);
1086  }
1087 
1088  //ECAL entrance
1089  theBremParticle.propagateToEcalEntrance(false);
1090  if (theBremParticle.getSuccess() != 0) {
1091  brem.addPoint(PFTrajectoryPoint(-1,
1092  PFTrajectoryPoint::ECALEntrance,
1093  math::XYZPoint(theBremParticle.particle().vertex()),
1094  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1095  double ecalShowerDepth =
1096  PFCluster::getDepthCorrection(theBremParticle.particle().momentum().E(), isBelowPS, false);
1097 
1098  math::XYZPoint meanShower =
1099  math::XYZPoint(theBremParticle.particle().vertex()) +
1100  math::XYZTLorentzVector(theBremParticle.particle().momentum()).Vect().Unit() * ecalShowerDepth;
1101 
1102  brem.addPoint(PFTrajectoryPoint(-1,
1103  PFTrajectoryPoint::ECALShowerMax,
1104  meanShower,
1105  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1106  } else {
1107  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1108  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1109  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE ECAL HAS FAILED\n"
1110  << "theBremParticle.particle() pt,eta,cos2ThetaV: "
1111  << theBremParticle.particle().pt() << " " << theBremParticle.particle().eta()
1112  << " " << theBremParticle.particle().cos2ThetaV();
1113  PFTrajectoryPoint dummyECAL;
1114  brem.addPoint(dummyECAL);
1115  PFTrajectoryPoint dummyMaxSh;
1116  brem.addPoint(dummyMaxSh);
1117  }
1118 
1119  //HCAL entrance
1120  theBremParticle.propagateToHcalEntrance(false);
1121  if (theBremParticle.getSuccess() != 0)
1122  brem.addPoint(PFTrajectoryPoint(-1,
1123  PFTrajectoryPoint::HCALEntrance,
1124  math::XYZPoint(theBremParticle.particle().vertex()),
1125  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1126  else {
1127  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1128  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1129  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL ENTRANCE HAS FAILED";
1130  PFTrajectoryPoint dummyHCALentrance;
1131  brem.addPoint(dummyHCALentrance);
1132  }
1133  //HCAL exit
1134  theBremParticle.propagateToHcalExit(false);
1135  if (theBremParticle.getSuccess() != 0)
1136  brem.addPoint(PFTrajectoryPoint(-1,
1137  PFTrajectoryPoint::HCALExit,
1138  math::XYZPoint(theBremParticle.particle().vertex()),
1139  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1140  else {
1141  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) &&
1142  theBremParticle.particle().cos2ThetaV() < cos2ThetaV_Endcap_HiEnd_)
1143  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE HCAL EXIT HAS FAILED";
1144  PFTrajectoryPoint dummyHCALexit;
1145  brem.addPoint(dummyHCALexit);
1146  }
1147 
1148  //HO Layer0
1149  if (abs(theBremParticle.particle().vertex().z()) < 550) {
1150  theBremParticle.setMagneticField(0);
1151  theBremParticle.propagateToHOLayer(false);
1152  if (theBremParticle.getSuccess() != 0)
1153  brem.addPoint(PFTrajectoryPoint(-1,
1154  PFTrajectoryPoint::HOLayer,
1155  math::XYZPoint(theBremParticle.particle().vertex()),
1156  math::XYZTLorentzVector(theBremParticle.particle().momentum())));
1157  else {
1158  if ((dp_tang > 5.) && ((dp_tang / sdp_tang) > 3) && abs(theBremParticle.particle().Z()) < 700.25)
1159  LogWarning("PFTrackTransformer") << "BREM " << brem << " PROPAGATION TO THE H0 HAS FAILED";
1160  PFTrajectoryPoint dummyHOLayer;
1161  brem.addPoint(dummyHOLayer);
1162  }
1163  }
1164  brem.calculatePositionREP();
1165  pftrack.addBrem(brem);
1166  iTrajPos++;
1167  }
1168 
1169  return true;
1170 }
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
BaseParticlePropagator::propagateToEcalEntrance
bool propagateToEcalEntrance(bool first=true)
Definition: BaseParticlePropagator.cc:450
reco::PFTrack::calculatePositionREP
void calculatePositionREP()
Definition: PFTrack.cc:52
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
edm::LogInfo
Definition: MessageLogger.h:254
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
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
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:29
reco::PFRecTrack::algoType
unsigned int algoType() const
Definition: PFRecTrack.h:41
edm::LogWarning
Definition: MessageLogger.h:141
PFGsfHelper::sigmafittedDP
double sigmafittedDP() const
Definition: PFGsfHelper.cc:131
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::LogError
Definition: MessageLogger.h:183
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
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
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:17
reco::PFRecTrack
reconstructed track used as an input to particle flow
Definition: PFRecTrack.h:22
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
HLT_2018_cff.track
track
Definition: HLT_2018_cff.py:10352
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