Util/src/ElectronEpCombination.cc

#include "ElectronEpCombination.h"
#include "ElectronMomentumRegression.h"
#include "TLorentzVector.h"

//these classifications are taken from:
//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DataFormats/EgammaCandidates/interface/GsfElectron.h?revision=1.59&view=markup
enum Classification { GSF_ELECTRON_UNKNOWN=-1, GSF_ELECTRON_GOLDEN=0, GSF_ELECTRON_BIGBREM=1, GSF_ELECTRON_BADTRACK=2, GSF_ELECTRON_SHOWERING=3, GSF_ELECTRON_GAP=4 } ;

//this is based on:
//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/EgammaAnalysis/ElectronTools/src/PatElectronEnergyCalibrator.cc?revision=1.2.2.3&view=markup

// void do_ep_combination(ControlFlags &ctrl, mithep::Electron * ele, const float regression_energy, const float regression_energy_error)
// {
//   float newEnergyError_ = regression_energy_error;
//   float regressionMomentum = regression_energy;

//   int elClass = ele->Classification();
//   float trackMomentum  = ele->PIn() ;
//   float errorTrackMomentum_ = 999. ;
  
//   // the electron's track momentum error was not available in bambu versions less than 029
//   if(ele->TrackMomentumError() > 0)
//     errorTrackMomentum_ = ele->TrackMomentumError();
//   else if ( ele->GsfTrk()->PtErr() > 0)
//     errorTrackMomentum_ = ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta());
//   else
//     assert(0);

//   float finalMomentum = ele->E(); // initial
//   float finalMomentumError = 999.;
  
//   // first check for large errors
//   if (errorTrackMomentum_/trackMomentum > 0.5 && newEnergyError_/regressionMomentum <= 0.5) {
//     finalMomentum = regressionMomentum;    finalMomentumError = newEnergyError_;
//   }
//   else if (errorTrackMomentum_/trackMomentum <= 0.5 && newEnergyError_/regressionMomentum > 0.5){  
//     finalMomentum = trackMomentum;  finalMomentumError = errorTrackMomentum_;
//   }
//   else if (errorTrackMomentum_/trackMomentum > 0.5 && newEnergyError_/regressionMomentum > 0.5){
//     if (errorTrackMomentum_/trackMomentum < newEnergyError_/regressionMomentum) {
//       finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum_;
//     }
//     else{
//       finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_;
//     }
//   }
//   // then apply the combination algorithm
//   else {
//      // calculate E/p and corresponding error
//     float eOverP = regressionMomentum / trackMomentum;
//     float errorEOverP = sqrt(
//                           (newEnergyError_/trackMomentum)*(newEnergyError_/trackMomentum) +
//                           (regressionMomentum*errorTrackMomentum_/trackMomentum/trackMomentum)*
//                           (regressionMomentum*errorTrackMomentum_/trackMomentum/trackMomentum));


//     bool eleIsNotInCombination = false ;
//     if ( (eOverP  > 1 + 2.5*errorEOverP) || (eOverP  < 1 - 2.5*errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3) )
//       { eleIsNotInCombination = true ; }
//      if (eleIsNotInCombination)
//        {
//       if (eOverP > 1)
//         { finalMomentum = regressionMomentum ; finalMomentumError = newEnergyError_ ; }
//       else
//         {
//           if (elClass == GSF_ELECTRON_GOLDEN)
//             { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_; }
//           if (elClass == GSF_ELECTRON_BIGBREM)
//             {
//               if (regressionMomentum<36)
//                 { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_ ; }
//               else
//                 { finalMomentum = regressionMomentum ; finalMomentumError = newEnergyError_ ; }
//             }
//           if (elClass == GSF_ELECTRON_BADTRACK)
//             { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_ ; }
//           if (elClass == GSF_ELECTRON_SHOWERING)
//             {
//               if (regressionMomentum<30)
//                 { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_; }
//               else
//                 { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_;}
//             }
//           if (elClass == GSF_ELECTRON_GAP)
//             {
//               if (regressionMomentum<60)
//                 { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_ ; }
//               else
//                 { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_ ; }
//             }
//         }
//        }     
//      else 
//        {
//       // combination
//       finalMomentum = (regressionMomentum/newEnergyError_/newEnergyError_ + trackMomentum/errorTrackMomentum_/errorTrackMomentum_) /
//         (1/newEnergyError_/newEnergyError_ + 1/errorTrackMomentum_/errorTrackMomentum_);
//       float finalMomentumVariance = 1 / (1/newEnergyError_/newEnergyError_ + 1/errorTrackMomentum_/errorTrackMomentum_);
//       finalMomentumError = sqrt(finalMomentumVariance);
//        }
//   }  

//   TLorentzVector oldMomentum = TLorentzVector(ele->Px(),ele->Py(),ele->Pz(),ele->E()) ;

//   if(ctrl.debug) cout << "    e/p comb " << setw(12) << setprecision(7) << regression_energy << " -> " << setw(12) << finalMomentum << setprecision(5) << endl;
  
//   TLorentzVector newMomentum_ 
//     = TLorentzVector
//    ( oldMomentum.X()*finalMomentum/oldMomentum.T(),
//      oldMomentum.Y()*finalMomentum/oldMomentum.T(),
//      oldMomentum.Z()*finalMomentum/oldMomentum.T(),
//      finalMomentum ) ;
  
//   ele->SetPtEtaPhi(newMomentum_.Pt(),newMomentum_.Eta(),newMomentum_.Phi() );
  
// }
void do_ep_combination(ControlFlags &ctrl, mithep::Electron * ele, const float regression_energy, const float regression_energy_error)
//
// Basically a cut and paste job from UserCode/Mangano/WWAnalysis/AnalysisStep/src/PatElectronEnergyCalibrator
//
{
  float regressionMomentumError = regression_energy_error; // change names for consistency with Boris's code
  float regressionMomentum = regression_energy;

  int elClass = ele->Classification();
  float trackMomentum  = ele->PIn() ;
  // the electron's track momentum error was not available in bambu versions less than 029.
  // So, use an approximation
  float trackMomentumError(0);
  if(ele->TrackMomentumError() > 0)
    trackMomentumError = ele->TrackMomentumError();
  else if ( ele->GsfTrk()->PtErr() > 0) {
    trackMomentumError = ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta());
  } else
    assert(0);

  // NOTE: not friggin close at all
  // if(ctrl.debug) cout << "    tme: " << ele->TrackMomentumError() << " hack: " << ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta()) << endl;

  if(ctrl.debug) cout << setprecision(7) << "e/p comb:: regr E: " << setw(12) << regressionMomentum << " (" << setw(12) << regressionMomentumError << ") track: " << setw(12) << trackMomentum << setw(12) << " (" << trackMomentumError << ")" << setprecision(5) << endl;

  float finalMomentum = regressionMomentum;
  float finalMomentumError = regression_energy_error;
  
  //------- HERE BEGINS THE COMMON E-P COMBINATION PART  -------------------
  // first check for large errors
  if (trackMomentumError/trackMomentum > 0.5 && regressionMomentumError/regressionMomentum <= 0.5) {
    finalMomentum = regressionMomentum;    finalMomentumError = regressionMomentumError;
  }
  else if (trackMomentumError/trackMomentum <= 0.5 && regressionMomentumError/regressionMomentum > 0.5){  
    finalMomentum = trackMomentum;  finalMomentumError = trackMomentumError;
  }
  else if (trackMomentumError/trackMomentum > 0.5 && regressionMomentumError/regressionMomentum > 0.5){
    if (trackMomentumError/trackMomentum < regressionMomentumError/regressionMomentum) {
      finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;
    }
    else{
      finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
    }
  }
  // then apply the combination algorithm
  else {
    // calculate E/p and corresponding error
    float eOverP = regressionMomentum / trackMomentum;
    float errorEOverP = sqrt(
                             (regressionMomentumError/trackMomentum)*(regressionMomentumError/trackMomentum) +
                             (regressionMomentum*trackMomentumError/trackMomentum/trackMomentum)*
                             (regressionMomentum*trackMomentumError/trackMomentum/trackMomentum));


#if ROOT_VERSION_CODE >=  ROOT_VERSION(5,30,00)
    // ======== HERE IS THE 5XY (NEW) E-P COMBINATION
    bool eleIsNotInCombination = false ;
    if ( (eOverP  > 1 + 2.5*errorEOverP) || (eOverP  < 1 - 2.5*errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3) )
      { eleIsNotInCombination = true ; }
    if (eleIsNotInCombination)
      {
        if (eOverP > 1)
          { finalMomentum = regressionMomentum ; finalMomentumError = regressionMomentumError ; }
        else
          {
            if (elClass == GSF_ELECTRON_GOLDEN)
              { finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError; }
            if (elClass == GSF_ELECTRON_BIGBREM)
              {
                if (regressionMomentum<36)
                  { finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError ; }
                else
                  { finalMomentum = regressionMomentum ; finalMomentumError = regressionMomentumError ; }
              }
            if (elClass == GSF_ELECTRON_BADTRACK)  //for 53X
              { finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError ; }
            if (elClass == GSF_ELECTRON_SHOWERING)
              {
                if (regressionMomentum<30)
                  { finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError; }
                else
                  { finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;}
              }
            if (elClass == GSF_ELECTRON_GAP)
              {
                if (regressionMomentum<60)
                  { finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError ; }
                else
                  { finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError ; }
              }
          }
      }
 
    else
      {
        // combination
        finalMomentum = (regressionMomentum/regressionMomentumError/regressionMomentumError + trackMomentum/trackMomentumError/trackMomentumError) /
          (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
        float finalMomentumVariance = 1 / (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
        finalMomentumError = sqrt(finalMomentumVariance);
      }
    // ======== FINISH 5XY E-P COMBINATION
#else
    // ======== HERE IS THE 4XY (NEW) E-P COMBINATION
    //
    // NOTE: this will never get compiled the way we're running it now... but I'm not sure if we're supposed to be using it for 2011, so here it stays
    //
    if ( eOverP  > 1 + 2.5*errorEOverP )
      {
        finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
        if ((elClass==GSF_ELECTRON_GOLDEN) && ele->IsEB() && (eOverP<1.15))
          {
            if (regressionMomentum<15) {finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError;}
          }
      }
    else if ( eOverP < 1 - 2.5*errorEOverP )
      {
        finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
        if (elClass==GSF_ELECTRON_SHOWERING)
          {
            if (ele->IsEB())
              {
                if(regressionMomentum<18) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
              }
            else if (ele->IsEE())
              {
                if(regressionMomentum<13) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
              }
            else
              assert(0);
          }
        else if (ele->IsEBEEGap())
          {
            if(regressionMomentum<60) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
          }
      }
    else 
      {
        // combination
        finalMomentum = (regressionMomentum/regressionMomentumError/regressionMomentumError + trackMomentum/trackMomentumError/trackMomentumError) /
          (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
        float finalMomentumVariance = 1 / (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
        finalMomentumError = sqrt(finalMomentumVariance);
      }        
    // ======== FINISH 4XY E-P COMBINATION
#endif


  } // end block where neither error is super huge

  TLorentzVector oldMomentum = TLorentzVector(ele->Px(),ele->Py(),ele->Pz(),ele->E()) ;

  if(ctrl.debug) cout << "    e/p comb " << setw(12) << setprecision(7) << regression_energy << " -> " << setw(12) << finalMomentum << setprecision(5) << endl;
  
  TLorentzVector newMomentum_ 
    = TLorentzVector
    ( oldMomentum.X()*finalMomentum/oldMomentum.T(),
      oldMomentum.Y()*finalMomentum/oldMomentum.T(),
      oldMomentum.Z()*finalMomentum/oldMomentum.T(),
      finalMomentum ) ;
  
  ele->SetPtEtaPhi(newMomentum_.Pt(),newMomentum_.Eta(),newMomentum_.Phi() );
  
}
Revision:	1.5
Committed:	Mon Dec 17 17:23:06 2012 UTC (12 years, 5 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
CVS Tags:	compiled
Changes since 1.4:	+234 -70 lines
Log Message:	* empty log message *
#	Content
1	#include "ElectronEpCombination.h"
2	#include "ElectronMomentumRegression.h"
3	#include "TLorentzVector.h"
4
5	//these classifications are taken from:
6	//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DataFormats/EgammaCandidates/interface/GsfElectron.h?revision=1.59&view=markup
7	enum Classification { GSF_ELECTRON_UNKNOWN=-1, GSF_ELECTRON_GOLDEN=0, GSF_ELECTRON_BIGBREM=1, GSF_ELECTRON_BADTRACK=2, GSF_ELECTRON_SHOWERING=3, GSF_ELECTRON_GAP=4 } ;
8
9	//this is based on:
10	//http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/EgammaAnalysis/ElectronTools/src/PatElectronEnergyCalibrator.cc?revision=1.2.2.3&view=markup
11
12	// void do_ep_combination(ControlFlags &ctrl, mithep::Electron * ele, const float regression_energy, const float regression_energy_error)
13	// {
14	// float newEnergyError_ = regression_energy_error;
15	// float regressionMomentum = regression_energy;
16
17	// int elClass = ele->Classification();
18	// float trackMomentum = ele->PIn() ;
19	// float errorTrackMomentum_ = 999. ;
20
21	// // the electron's track momentum error was not available in bambu versions less than 029
22	// if(ele->TrackMomentumError() > 0)
23	// errorTrackMomentum_ = ele->TrackMomentumError();
24	// else if ( ele->GsfTrk()->PtErr() > 0)
25	// errorTrackMomentum_ = ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta());
26	// else
27	// assert(0);
28
29	// float finalMomentum = ele->E(); // initial
30	// float finalMomentumError = 999.;
31
32	// // first check for large errors
33	// if (errorTrackMomentum_/trackMomentum > 0.5 && newEnergyError_/regressionMomentum <= 0.5) {
34	// finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_;
35	// }
36	// else if (errorTrackMomentum_/trackMomentum <= 0.5 && newEnergyError_/regressionMomentum > 0.5){
37	// finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum_;
38	// }
39	// else if (errorTrackMomentum_/trackMomentum > 0.5 && newEnergyError_/regressionMomentum > 0.5){
40	// if (errorTrackMomentum_/trackMomentum < newEnergyError_/regressionMomentum) {
41	// finalMomentum = trackMomentum; finalMomentumError = errorTrackMomentum_;
42	// }
43	// else{
44	// finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_;
45	// }
46	// }
47	// // then apply the combination algorithm
48	// else {
49	// // calculate E/p and corresponding error
50	// float eOverP = regressionMomentum / trackMomentum;
51	// float errorEOverP = sqrt(
52	// (newEnergyError_/trackMomentum)*(newEnergyError_/trackMomentum) +
53	// (regressionMomentumerrorTrackMomentum_/trackMomentum/trackMomentum)
54	// (regressionMomentum*errorTrackMomentum_/trackMomentum/trackMomentum));
55
56
57	// bool eleIsNotInCombination = false ;
58	// if ( (eOverP > 1 + 2.5errorEOverP) \|\| (eOverP < 1 - 2.5errorEOverP) \|\| (eOverP < 0.8) \|\| (eOverP > 1.3) )
59	// { eleIsNotInCombination = true ; }
60	// if (eleIsNotInCombination)
61	// {
62	// if (eOverP > 1)
63	// { finalMomentum = regressionMomentum ; finalMomentumError = newEnergyError_ ; }
64	// else
65	// {
66	// if (elClass == GSF_ELECTRON_GOLDEN)
67	// { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_; }
68	// if (elClass == GSF_ELECTRON_BIGBREM)
69	// {
70	// if (regressionMomentum<36)
71	// { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_ ; }
72	// else
73	// { finalMomentum = regressionMomentum ; finalMomentumError = newEnergyError_ ; }
74	// }
75	// if (elClass == GSF_ELECTRON_BADTRACK)
76	// { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_ ; }
77	// if (elClass == GSF_ELECTRON_SHOWERING)
78	// {
79	// if (regressionMomentum<30)
80	// { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_; }
81	// else
82	// { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_;}
83	// }
84	// if (elClass == GSF_ELECTRON_GAP)
85	// {
86	// if (regressionMomentum<60)
87	// { finalMomentum = trackMomentum ; finalMomentumError = errorTrackMomentum_ ; }
88	// else
89	// { finalMomentum = regressionMomentum; finalMomentumError = newEnergyError_ ; }
90	// }
91	// }
92	// }
93	// else
94	// {
95	// // combination
96	// finalMomentum = (regressionMomentum/newEnergyError_/newEnergyError_ + trackMomentum/errorTrackMomentum_/errorTrackMomentum_) /
97	// (1/newEnergyError_/newEnergyError_ + 1/errorTrackMomentum_/errorTrackMomentum_);
98	// float finalMomentumVariance = 1 / (1/newEnergyError_/newEnergyError_ + 1/errorTrackMomentum_/errorTrackMomentum_);
99	// finalMomentumError = sqrt(finalMomentumVariance);
100	// }
101	// }
102
103	// TLorentzVector oldMomentum = TLorentzVector(ele->Px(),ele->Py(),ele->Pz(),ele->E()) ;
104
105	// if(ctrl.debug) cout << " e/p comb " << setw(12) << setprecision(7) << regression_energy << " -> " << setw(12) << finalMomentum << setprecision(5) << endl;
106
107	// TLorentzVector newMomentum_
108	// = TLorentzVector
109	// ( oldMomentum.X()*finalMomentum/oldMomentum.T(),
110	// oldMomentum.Y()*finalMomentum/oldMomentum.T(),
111	// oldMomentum.Z()*finalMomentum/oldMomentum.T(),
112	// finalMomentum ) ;
113
114	// ele->SetPtEtaPhi(newMomentum_.Pt(),newMomentum_.Eta(),newMomentum_.Phi() );
115
116	// }
117	void do_ep_combination(ControlFlags &ctrl, mithep::Electron * ele, const float regression_energy, const float regression_energy_error)
118	//
119	// Basically a cut and paste job from UserCode/Mangano/WWAnalysis/AnalysisStep/src/PatElectronEnergyCalibrator
120	//
121	{
122	float regressionMomentumError = regression_energy_error; // change names for consistency with Boris's code
123	float regressionMomentum = regression_energy;
124
125	int elClass = ele->Classification();
126	float trackMomentum = ele->PIn() ;
127	// the electron's track momentum error was not available in bambu versions less than 029.
128	// So, use an approximation
129	float trackMomentumError(0);
130	if(ele->TrackMomentumError() > 0)
131	trackMomentumError = ele->TrackMomentumError();
132	else if ( ele->GsfTrk()->PtErr() > 0) {
133	trackMomentumError = ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta());
134	} else
135	assert(0);
136
137	// NOTE: not friggin close at all
138	// if(ctrl.debug) cout << " tme: " << ele->TrackMomentumError() << " hack: " << ele->GsfTrk()->PtErr()*cosh(ele->GsfTrk()->Eta()) << endl;
139
140	if(ctrl.debug) cout << setprecision(7) << "e/p comb:: regr E: " << setw(12) << regressionMomentum << " (" << setw(12) << regressionMomentumError << ") track: " << setw(12) << trackMomentum << setw(12) << " (" << trackMomentumError << ")" << setprecision(5) << endl;
141
142	float finalMomentum = regressionMomentum;
143	float finalMomentumError = regression_energy_error;
144
145	//------- HERE BEGINS THE COMMON E-P COMBINATION PART -------------------
146	// first check for large errors
147	if (trackMomentumError/trackMomentum > 0.5 && regressionMomentumError/regressionMomentum <= 0.5) {
148	finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
149	}
150	else if (trackMomentumError/trackMomentum <= 0.5 && regressionMomentumError/regressionMomentum > 0.5){
151	finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;
152	}
153	else if (trackMomentumError/trackMomentum > 0.5 && regressionMomentumError/regressionMomentum > 0.5){
154	if (trackMomentumError/trackMomentum < regressionMomentumError/regressionMomentum) {
155	finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;
156	}
157	else{
158	finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
159	}
160	}
161	// then apply the combination algorithm
162	else {
163	// calculate E/p and corresponding error
164	float eOverP = regressionMomentum / trackMomentum;
165	float errorEOverP = sqrt(
166	(regressionMomentumError/trackMomentum)*(regressionMomentumError/trackMomentum) +
167	(regressionMomentumtrackMomentumError/trackMomentum/trackMomentum)
168	(regressionMomentum*trackMomentumError/trackMomentum/trackMomentum));
169
170
171	#if ROOT_VERSION_CODE >= ROOT_VERSION(5,30,00)
172	// ======== HERE IS THE 5XY (NEW) E-P COMBINATION
173	bool eleIsNotInCombination = false ;
174	if ( (eOverP > 1 + 2.5errorEOverP) \|\| (eOverP < 1 - 2.5errorEOverP) \|\| (eOverP < 0.8) \|\| (eOverP > 1.3) )
175	{ eleIsNotInCombination = true ; }
176	if (eleIsNotInCombination)
177	{
178	if (eOverP > 1)
179	{ finalMomentum = regressionMomentum ; finalMomentumError = regressionMomentumError ; }
180	else
181	{
182	if (elClass == GSF_ELECTRON_GOLDEN)
183	{ finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError; }
184	if (elClass == GSF_ELECTRON_BIGBREM)
185	{
186	if (regressionMomentum<36)
187	{ finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError ; }
188	else
189	{ finalMomentum = regressionMomentum ; finalMomentumError = regressionMomentumError ; }
190	}
191	if (elClass == GSF_ELECTRON_BADTRACK) //for 53X
192	{ finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError ; }
193	if (elClass == GSF_ELECTRON_SHOWERING)
194	{
195	if (regressionMomentum<30)
196	{ finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError; }
197	else
198	{ finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;}
199	}
200	if (elClass == GSF_ELECTRON_GAP)
201	{
202	if (regressionMomentum<60)
203	{ finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError ; }
204	else
205	{ finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError ; }
206	}
207	}
208	}
209
210	else
211	{
212	// combination
213	finalMomentum = (regressionMomentum/regressionMomentumError/regressionMomentumError + trackMomentum/trackMomentumError/trackMomentumError) /
214	(1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
215	float finalMomentumVariance = 1 / (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
216	finalMomentumError = sqrt(finalMomentumVariance);
217	}
218	// ======== FINISH 5XY E-P COMBINATION
219	#else
220	// ======== HERE IS THE 4XY (NEW) E-P COMBINATION
221	//
222	// NOTE: this will never get compiled the way we're running it now... but I'm not sure if we're supposed to be using it for 2011, so here it stays
223	//
224	if ( eOverP > 1 + 2.5*errorEOverP )
225	{
226	finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
227	if ((elClass==GSF_ELECTRON_GOLDEN) && ele->IsEB() && (eOverP<1.15))
228	{
229	if (regressionMomentum<15) {finalMomentum = trackMomentum ; finalMomentumError = trackMomentumError;}
230	}
231	}
232	else if ( eOverP < 1 - 2.5*errorEOverP )
233	{
234	finalMomentum = regressionMomentum; finalMomentumError = regressionMomentumError;
235	if (elClass==GSF_ELECTRON_SHOWERING)
236	{
237	if (ele->IsEB())
238	{
239	if(regressionMomentum<18) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
240	}
241	else if (ele->IsEE())
242	{
243	if(regressionMomentum<13) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
244	}
245	else
246	assert(0);
247	}
248	else if (ele->IsEBEEGap())
249	{
250	if(regressionMomentum<60) {finalMomentum = trackMomentum; finalMomentumError = trackMomentumError;}
251	}
252	}
253	else
254	{
255	// combination
256	finalMomentum = (regressionMomentum/regressionMomentumError/regressionMomentumError + trackMomentum/trackMomentumError/trackMomentumError) /
257	(1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
258	float finalMomentumVariance = 1 / (1/regressionMomentumError/regressionMomentumError + 1/trackMomentumError/trackMomentumError);
259	finalMomentumError = sqrt(finalMomentumVariance);
260	}
261	// ======== FINISH 4XY E-P COMBINATION
262	#endif
263
264
265	} // end block where neither error is super huge
266
267	TLorentzVector oldMomentum = TLorentzVector(ele->Px(),ele->Py(),ele->Pz(),ele->E()) ;
268
269	if(ctrl.debug) cout << " e/p comb " << setw(12) << setprecision(7) << regression_energy << " -> " << setw(12) << finalMomentum << setprecision(5) << endl;
270
271	TLorentzVector newMomentum_
272	= TLorentzVector
273	( oldMomentum.X()*finalMomentum/oldMomentum.T(),
274	oldMomentum.Y()*finalMomentum/oldMomentum.T(),
275	oldMomentum.Z()*finalMomentum/oldMomentum.T(),
276	finalMomentum ) ;
277
278	ele->SetPtEtaPhi(newMomentum_.Pt(),newMomentum_.Eta(),newMomentum_.Phi() );
279
280	}