CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions
ElectronMomentumCorrector Class Reference

#include <ElectronMomentumCorrector.h>

Public Member Functions

void correct (reco::GsfElectron &, TrajectoryStateOnSurface &)
 
 ElectronMomentumCorrector ()
 
 ~ElectronMomentumCorrector ()
 

Private Member Functions

float energyError (float E, float *par) const
 

Detailed Description

Definition at line 19 of file ElectronMomentumCorrector.h.

Constructor & Destructor Documentation

ElectronMomentumCorrector::ElectronMomentumCorrector ( )
inline

Definition at line 25 of file ElectronMomentumCorrector.h.

25 {}
ElectronMomentumCorrector::~ElectronMomentumCorrector ( )
inline

Definition at line 27 of file ElectronMomentumCorrector.h.

27 {}

Member Function Documentation

void ElectronMomentumCorrector::correct ( reco::GsfElectron &  electron,
TrajectoryStateOnSurface &  vtxTsos 
)

Definition at line 28 of file ElectronMomentumCorrector.cc.

References reco::GsfElectron::BADTRACK, reco::GsfElectron::BIGBREM, reco::GsfElectron::classification(), reco::GsfElectron::correctedEcalEnergy(), reco::GsfElectron::correctedEcalEnergyError(), reco::GsfElectron::correctMomentum(), edm::false, reco::GsfElectron::GAP, reco::GsfElectron::GOLDEN, reco::GsfElectron::isEB(), reco::GsfElectron::isEE(), GaussianSumUtilities1D::mode(), MultiGaussianStateTransform::multiState1D(), reco::GsfElectron::p4(), reco::GsfElectron::P4_COMBINATION, reco::GsfElectron::p4Error(), reco::return(), pileupReCalc_HLTpaths::scale, reco::GsfElectron::SHOWERING, mathSSE::sqrt(), reco::btau::trackMomentum, reco::GsfElectron::trackMomentumAtVtx(), funct::true, reco::GsfElectron::UNKNOWN, and SingleGaussianState1D::variance().

Referenced by GsfElectronAlgo::createElectron().

29  {
30  if (electron.p4Error(reco::GsfElectron::P4_COMBINATION)!= 999.)
31  {
32  edm::LogWarning("ElectronMomentumCorrector::correct")<<"already done" ;
33  return ;
34  }
35 
36  //math::XYZTLorentzVector newMomentum = electron.p4() ; // default
37  int elClass = electron.classification() ;
38 
39  // irrelevant classification
40  if ( (elClass <= reco::GsfElectron::UNKNOWN) ||
41  (elClass>reco::GsfElectron::GAP) )
42  {
43  edm::LogWarning("ElectronMomentumCorrector::correct")<<"unexpected classification" ;
44  return ;
45  }
46 
47 
48  //=======================================================================================
49  // cluster energy
50  //=======================================================================================
51 
52  float scEnergy = electron.correctedEcalEnergy() ;
53  float errorEnergy = electron.correctedEcalEnergyError() ;
54 
55 
56  //=======================================================================================
57  // track momentum
58  //=======================================================================================
59 
60  // basic values
61  float trackMomentum = electron.trackMomentumAtVtx().R() ;
62  //float errorTrackMomentum = 999. ;
63 
64  // tracker momentum scale corrections (Mykhailo Dalchenko)
65  double scale=1.;
66  if (electron.isEB()){
67  if (elClass==0) {scale = 1./(0.00104*sqrt(trackMomentum)+1);}
68  if (elClass==1) {scale = 1./(0.0017*sqrt(trackMomentum)+0.9986);}
69  if (elClass==3) {scale = 1./(1.004 - 0.00021*trackMomentum);}
70  if (elClass==4) {scale = 0.995;}
71  } else if (electron.isEE()){
72  if (elClass==3){scale = 1./(1.01432-0.00201872*trackMomentum+0.0000142621*trackMomentum*trackMomentum);}
73  if (elClass==4){scale = 1./(0.996859-0.000345347*trackMomentum);}
74  }
75  if (scale<0.) scale = 1.; // CC added protection
76  trackMomentum = trackMomentum*scale ;
77 
78  // error (must be done after trackMomentum rescaling)
79  MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(vtxTsos,0)) ;
80  GaussianSumUtilities1D qpUtils(qpState) ;
81  float errorTrackMomentum = trackMomentum*trackMomentum*sqrt(qpUtils.mode().variance()) ;
82 
83 
84  //=======================================================================================
85  // combination
86  //=======================================================================================
87 
88  float finalMomentum = electron.p4().t(); // initial
89  float finalMomentumError = 999.;
90 
91  // first check for large errors
92 
93  if (errorTrackMomentum/trackMomentum > 0.5 && errorEnergy/scEnergy <= 0.5) {
94  finalMomentum = scEnergy; finalMomentumError = errorEnergy;
95  }
96  else if (errorTrackMomentum/trackMomentum <= 0.5 && errorEnergy/scEnergy > 0.5){
97  finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum;
98  }
99  else if (errorTrackMomentum/trackMomentum > 0.5 && errorEnergy/scEnergy > 0.5){
100  if (errorTrackMomentum/trackMomentum < errorEnergy/scEnergy) {
101  finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum;
102  }
103  else{
104  finalMomentum = scEnergy; finalMomentumError = errorEnergy;
105  }
106  }
107 
108  // then apply the combination algorithm
109  else {
110 
111  // calculate E/p and corresponding error
112  float eOverP = scEnergy / trackMomentum;
113  float errorEOverP = sqrt(
114  (errorEnergy/trackMomentum)*(errorEnergy/trackMomentum) +
115  (scEnergy*errorTrackMomentum/trackMomentum/trackMomentum)*
116  (scEnergy*errorTrackMomentum/trackMomentum/trackMomentum));
117 
118 // if ( eOverP > 1 + 2.5*errorEOverP )
119 // {
120 // finalMomentum = scEnergy; finalMomentumError = errorEnergy;
121 // if ((elClass==reco::GsfElectron::GOLDEN) && electron.isEB() && (eOverP<1.15))
122 // {
123 // if (scEnergy<15) {finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum;}
124 // }
125 // }
126 // else if ( eOverP < 1 - 2.5*errorEOverP )
127 // {
128 // finalMomentum = scEnergy; finalMomentumError = errorEnergy;
129 // if (elClass==reco::GsfElectron::SHOWERING)
130 // {
131 // if (electron.isEB())
132 // {
133 // if (scEnergy<18)
134 // { finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum; }
135 // }
136 // else if (electron.isEE())
137 // {
138 // if (scEnergy<13)
139 // {finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum;}
140 // }
141 // else
142 // { edm::LogWarning("ElectronMomentumCorrector::correct")<<"nor barrel neither endcap electron ?!" ; }
143 // }
144 // else if (electron.isGap())
145 // {
146 // if (scEnergy<60)
147 // { finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum ; }
148 // }
149 // }
150 
151  bool eleIsNotInCombination = false ;
152  if ( (eOverP > 1 + 2.5*errorEOverP) || (eOverP < 1 - 2.5*errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3) )
153  { eleIsNotInCombination = true ; }
154  if (eleIsNotInCombination)
155  {
156  if (eOverP > 1)
157  { finalMomentum = scEnergy ; finalMomentumError = errorEnergy ; }
158  else
159  {
160  if (elClass == reco::GsfElectron::GOLDEN)
161  { finalMomentum = scEnergy; finalMomentumError = errorEnergy; }
162  if (elClass == reco::GsfElectron::BIGBREM)
163  {
164  if (scEnergy<36)
165  { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum ; }
166  else
167  { finalMomentum = scEnergy ; finalMomentumError = errorEnergy ; }
168  }
169  if (elClass == reco::GsfElectron::BADTRACK)
170  { finalMomentum = scEnergy; finalMomentumError = errorEnergy ; }
171  if (elClass == reco::GsfElectron::SHOWERING)
172  {
173  if (scEnergy<30)
174  { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum; }
175  else
176  { finalMomentum = scEnergy; finalMomentumError = errorEnergy;}
177  }
178  if (elClass == reco::GsfElectron::GAP)
179  {
180  if (scEnergy<60)
181  { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum ; }
182  else
183  { finalMomentum = scEnergy; finalMomentumError = errorEnergy ; }
184  }
185  }
186  }
187 
188  else
189  {
190  // combination
191  finalMomentum = (scEnergy/errorEnergy/errorEnergy + trackMomentum/errorTrackMomentum/errorTrackMomentum) /
192  (1/errorEnergy/errorEnergy + 1/errorTrackMomentum/errorTrackMomentum);
193  float finalMomentumVariance = 1 / (1/errorEnergy/errorEnergy + 1/errorTrackMomentum/errorTrackMomentum);
194  finalMomentumError = sqrt(finalMomentumVariance);
195  }
196 
197  }
198 
199 
200  //=======================================================================================
201  // final set
202  //=======================================================================================
203 
204  math::XYZTLorentzVector oldMomentum = electron.p4() ;
205  math::XYZTLorentzVector newMomentum = math::XYZTLorentzVector
206  ( oldMomentum.x()*finalMomentum/oldMomentum.t(),
207  oldMomentum.y()*finalMomentum/oldMomentum.t(),
208  oldMomentum.z()*finalMomentum/oldMomentum.t(),
209  finalMomentum ) ;
210 
211  electron.correctMomentum(newMomentum,errorTrackMomentum,finalMomentumError);
212 
213  }
reco::GsfElectron::p4
const LorentzVector & p4(P4Kind kind) const
Definition: GsfElectron.cc:204
reco::GsfElectron::correctMomentum
void correctMomentum(const LorentzVector &p4, float trackMomentumError, float p4Error)
Definition: GsfElectron.h:741
reco::GsfElectron::trackMomentumAtVtx
math::XYZVectorF trackMomentumAtVtx() const
Definition: GsfElectron.h:272
reco::GsfElectron::p4Error
float p4Error(P4Kind kind) const
Definition: GsfElectron.cc:216
reco::GsfElectron::isEE
bool isEE() const
Definition: GsfElectron.h:331
reco::GsfElectron::isEB
bool isEB() const
Definition: GsfElectron.h:330
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
MultiGaussianStateTransform::multiState1D
MultiGaussianState1D multiState1D(const std::vector< MultiGaussianState< N >::Vector > &, const std::vector< MultiGaussianState< N >::Matrix > &, const std::vector< double > &, unsigned int)
Definition: MultiGaussianStateTransform.cc:71
pileupReCalc_HLTpaths.scale
int scale
Definition: pileupReCalc_HLTpaths.py:111
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
reco::GsfElectron::correctedEcalEnergy
float correctedEcalEnergy() const
Definition: GsfElectron.h:720
reco::GsfElectron::classification
Classification classification() const
Definition: GsfElectron.h:642
reco::return
return(e1-e2)*(e1-e2)+dp *dp
reco::GsfElectron::correctedEcalEnergyError
float correctedEcalEnergyError() const
Definition: GsfElectron.h:721
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
float ElectronMomentumCorrector::energyError ( float  E,
float *  par 
) const
private
Generated for CMSSW Reference Manual by doxygen 1.8.5