RootCore Packages: /cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.4.31/ElectronPhotonShowerShapeFudgeTool/ElectronPhotonShowerShapeFudgeTool/TPhotonMCShifterTool.h Source File

00001 #ifndef TPhotonMCShifterTool_h
00002 #define TPhotonMCShifterTool_h
00003 
00004 /* @author Giovanni Marchiori <giovanni.marchiori@cern.ch>
00005  *
00006  *  Reimplementation of FudgeMCTool with photon fudge factors 
00007  *  stored in ROOT files instead of a C++ file
00008  *
00009  */
00010 
00011 #include <iostream>
00012 #include <stdio.h>
00013 #include <cmath>
00014 #include "TH2D.h"
00015 #include "TGraphErrors.h"
00016 #include "TSystem.h"
00017 
00018 namespace IDVAR
00019 {
00020     const int RHAD1  = 0;
00021     const int RHAD   = 1;
00022     const int E277   = 2;
00023     const int RETA   = 3;
00024     const int RPHI   = 4;
00025     const int WETA2  = 5;
00026     const int F1     = 6;
00027     const int FSIDE  = 7;
00028     const int WTOT   = 8;
00029     const int W1     = 9;
00030     const int DE     = 10;
00031     const int ERATIO = 11;
00032 }
00033 
00034 
00035 class TPhotonMCShifterTool 
00036 {
00037 
00038  public:
00039   
00040   TPhotonMCShifterTool();
00041   ~TPhotonMCShifterTool();
00042 
00043   //fudge shower shapes for preselection chosen.
00044   void FudgeShowers( double pt     ,
00045              double eta2   ,
00046              double& rhad1  ,
00047              double& rhad   ,
00048              double& e277   ,
00049              double& reta   ,
00050              double& rphi   ,
00051              double& weta2  ,
00052              double& f1     ,
00053              double& fside  ,
00054              double& wtot   ,
00055              double& w1     ,
00056              double& deltae ,
00057              double& eratio ,
00058              int    isConv  ,
00059              int    preselection=-999);
00060 
00061   //fudge shower shapes for preselection chosen.
00062   void FudgeShowers( float pt     ,
00063              float eta2   ,
00064              float& rhad1  ,
00065              float& rhad   ,
00066              float& e277   ,
00067              float& reta   ,
00068              float& rphi   ,
00069              float& weta2  ,
00070              float& f1     ,
00071              float& fside  ,
00072              float& wtot   ,
00073              float& w1     ,
00074              float& deltae ,
00075              float& eratio ,
00076              int    isConv  ,
00077              int    preselection=-999);
00078     
00079   // fudge showers using D3PD vectors (except for eratio)
00080   void FudgeShowers( std::vector<float> clE,
00081              std::vector<float> eta2   ,
00082              std::vector<float>& rhad1  ,
00083              std::vector<float>& rhad   ,
00084              std::vector<float>& e277   ,
00085              std::vector<float>& reta   ,
00086              std::vector<float>& rphi   ,
00087              std::vector<float>& weta2  ,
00088              std::vector<float>& f1     ,
00089              std::vector<float>& fside  ,
00090              std::vector<float>& wtot   ,
00091              std::vector<float>& w1     ,
00092              std::vector<float>& deltae ,
00093              std::vector<float>& eratio ,
00094              std::vector<int> isConv  ,
00095              int    preselection=-999);
00096 
00100   static std::vector<float>* getEratio(std::vector<float> emaxs1, std::vector<float> Emax2)
00101   {
00102       std::vector<float> *eratio = new std::vector<float>();
00103       for (unsigned int i = 0; i < emaxs1.size(); ++i)
00104     eratio->push_back(emaxs1[i] + Emax2[i] == 0 ? 0 : 
00105               (emaxs1[i] - Emax2[i])/(emaxs1[i] + Emax2[i]));
00106       return eratio;
00107   }
00108   
00109 
00110   //---KM ADD for interpolation 
00111   //No variable width arrays in Athena
00112   double centers[100];
00113   int NBINS;
00114   double temp_FF; // final FF passed 
00115   double temp_FF_err; // final FF passed 
00116   double FF_bin;  // FF of the bin
00117   double FF_int_bin; // FF to interpolate to
00118   double pt_ratio;   // ratio of center of pt bins for interpolation
00119   void fillBins(); // Function to get bin edges and widths
00120 
00121 
00122   double interpolate(double pt, double eta2, TH2D* hist){
00123     FF_bin = (hist->GetBinContent(hist->FindBin(pt, fabs(eta2))));
00124 
00125     //Edge cases - No including high pt to infinite
00126     if(pt >= centers[NBINS-1])
00127       return (hist->GetBinContent(hist->FindBin((centers[NBINS-1] -1), fabs(eta2))));
00128 
00129     if(pt < centers[0])
00130       return FF_bin;
00131 
00132     temp_FF = FF_bin;
00133 
00134     int bin_loc = (hist->FindBin(pt, fabs(eta2)) % (NBINS+2)) - 1;  // NBins+2 = including under and overflow, -1 is getting the bin array without underflow (see fillBins())
00135     
00136     if( pt < centers[bin_loc]){
00137       FF_int_bin = hist->GetBinContent(hist->FindBin(pt, fabs(eta2))-1);
00138       pt_ratio = (pt - centers[bin_loc-1]) / (centers[bin_loc] - centers[bin_loc-1]);
00139       temp_FF = FF_int_bin + (FF_bin - FF_int_bin)*pt_ratio;
00140     }
00141     else if(pt > centers[bin_loc])
00142     {
00143       FF_int_bin = hist->GetBinContent(hist->FindBin(pt, fabs(eta2))+1);
00144       pt_ratio = (pt - centers[bin_loc]) / (centers[bin_loc+1] - centers[bin_loc]);
00145       temp_FF = FF_bin + (FF_int_bin - FF_bin)*pt_ratio;
00146     }
00147     //testing printouts
00148     //std::cout << "pt : " << pt << " center of bin "<< centers[bin_loc] << " FF_bin: " << FF_bin << std::endl; 
00149     //std::cout << "lower bin FF: " << hist->GetBinContent(hist->FindBin(pt, fabs(eta2))-1) << " Upper bin FF: " << hist->GetBinContent(hist->FindBin(pt, fabs(eta2))+1) <<std::endl;
00150     //if( pt < centers[bin_loc]) std::cout << "using lower interpolation " << std::endl;
00151     //else std::cout << "using upper interpolation " << std::endl;
00152     //std::cout << "interpolated FF: " << temp_FF << std::endl;
00153     return temp_FF;
00154   }
00155 
00156   double interpolateErrors(double pt, double eta2, TH2D* hist){
00157     FF_bin = (hist->GetBinError(hist->FindBin(pt, fabs(eta2))));
00158 
00159     //Edge cases - No including high pt to infinite
00160     if(pt >= centers[NBINS-1])
00161       return (hist->GetBinError(hist->FindBin((centers[NBINS-1] -1), fabs(eta2))));
00162 
00163     if(pt < centers[0])
00164       return FF_bin;
00165 
00166     temp_FF_err = FF_bin;
00167 
00168     int bin_loc = (hist->FindBin(pt, fabs(eta2)) % (NBINS+2)) - 1;  // NBins+2 = including under and overflow, -1 is getting the bin array without underflow (see fillBins())
00169     
00170     if( pt < centers[bin_loc]){
00171       FF_int_bin = hist->GetBinError(hist->FindBin(pt, fabs(eta2))-1);
00172       pt_ratio = (pt - centers[bin_loc-1]) / (centers[bin_loc] - centers[bin_loc-1]);
00173       temp_FF = (FF_int_bin*FF_int_bin*(1-pt_ratio)*(1-pt_ratio) + FF_bin*FF_bin)/2;
00174       temp_FF_err = sqrt(temp_FF);
00175     }
00176     else if(pt > centers[bin_loc])
00177     {
00178       FF_int_bin = hist->GetBinError(hist->FindBin(pt, fabs(eta2))+1);
00179       pt_ratio = (pt - centers[bin_loc]) / (centers[bin_loc+1] - centers[bin_loc]);
00180       temp_FF = (FF_bin*FF_bin*(1-pt_ratio)*(1-pt_ratio) + FF_int_bin*FF_int_bin)/2;
00181       temp_FF_err = sqrt(temp_FF);
00182     }
00183     //testing printouts
00184     //std::cout << "pt : " << pt << " center of bin "<< centers[bin_loc] << " FF_bin: " << FF_bin << std::endl; 
00185     //std::cout << "lower bin FF: " << hist->GetBinError(hist->FindBin(pt, fabs(eta2))-1) << " Upper bin FF: " << hist->GetBinError(hist->FindBin(pt, fabs(eta2))+1) <<std::endl;
00186     //if( pt < centers[bin_loc]) std::cout << "using lower interpolation " << std::endl;
00187     //else std::cout << "using upper interpolation " << std::endl;
00188     //std::cout << "interpolated FF: " << temp_FF_err << std::endl;
00189     return temp_FF_err;
00190   }
00191 
00192 
00194 
00195   // get fudge factor for predefined preselection
00196   double GetFF_Rhad1 (double pt, double eta2, int conv){
00197     if (conv) return (interpolate(pt, eta2, h_c_rhad1));
00198     else      return (interpolate(pt, eta2, h_u_rhad1));};
00199 
00200   double GetFF_Rhad  (double pt, double eta2, int conv){
00201     if (conv) return (interpolate(pt, eta2, h_c_rhad));
00202     else      return (interpolate(pt, eta2, h_u_rhad));};
00203 
00204   double GetFF_E277  (double pt, double eta2, int conv){
00205     if (conv) return (interpolate(pt, eta2, h_c_e277));
00206     else      return (interpolate(pt, eta2, h_u_e277));};
00207 
00208   double GetFF_Reta  (double pt, double eta2, int conv){
00209     if (conv) return (interpolate(pt, eta2, h_c_reta));
00210     else      return (interpolate(pt, eta2, h_u_reta));};
00211 
00212   double GetFF_Rphi  (double pt, double eta2, int conv){
00213     if (conv) return (interpolate(pt, eta2, h_c_rphi));
00214     else      return (interpolate(pt, eta2, h_u_rphi));};
00215 
00216   double GetFF_Weta2 (double pt, double eta2, int conv){
00217     if (conv) return (interpolate(pt, eta2, h_c_weta2));
00218     else      return (interpolate(pt, eta2, h_u_weta2));};
00219 
00220   double GetFF_F1    (double pt, double eta2, int conv){
00221     if (conv) return (interpolate(pt, eta2, h_c_f1));
00222     else      return (interpolate(pt, eta2, h_u_f1));};
00223   
00224   double GetFF_DE    (double pt, double eta2, int conv){
00225     if (conv) return (interpolate(pt, eta2, h_c_de));
00226     else      return (interpolate(pt, eta2, h_u_de));};
00227   
00228   double GetFF_Eratio(double pt, double eta2, int conv){
00229     if (conv) return (interpolate(pt, eta2, h_c_eratio));
00230     else      return (interpolate(pt, eta2, h_u_eratio));};
00231 
00232   double GetFF_Fside (double pt, double eta2, int conv){
00233     if (conv) return (interpolate(pt, eta2, h_c_fside));
00234     else      return (interpolate(pt, eta2, h_u_fside));};
00235 
00236   double GetFF_Wtot  (double pt, double eta2, int conv){
00237     if (conv) return (interpolate(pt, eta2, h_c_wtot));
00238     else      return (interpolate(pt, eta2, h_u_wtot));};
00239 
00240   double GetFF_W1    (double pt, double eta2, int conv){
00241     if (conv) return (interpolate(pt, eta2, h_c_w1));
00242     else      return (interpolate(pt, eta2, h_u_w1));};
00243 
00244   double GetFFerr_Rhad1 (double pt, double eta2, int conv){
00245     if (conv) return (interpolateErrors(pt, eta2, h_c_rhad1));
00246     else      return (interpolateErrors(pt, eta2, h_u_rhad1));};
00247 
00248   double GetFFerr_Rhad  (double pt, double eta2, int conv){
00249     if (conv) return (interpolateErrors(pt, eta2, h_c_rhad));
00250     else      return (interpolateErrors(pt, eta2, h_u_rhad));};
00251 
00252   double GetFFerr_E277  (double pt, double eta2, int conv){
00253     if (conv) return (interpolateErrors(pt, eta2, h_c_e277));
00254     else      return (interpolateErrors(pt, eta2, h_u_e277));};
00255 
00256   double GetFFerr_Reta  (double pt, double eta2, int conv){
00257     if (conv) return (interpolateErrors(pt, eta2, h_c_reta));
00258     else      return (interpolateErrors(pt, eta2, h_u_reta));};
00259 
00260   double GetFFerr_Rphi  (double pt, double eta2, int conv){
00261     if (conv) return (interpolateErrors(pt, eta2, h_c_rphi));
00262     else      return (interpolateErrors(pt, eta2, h_u_rphi));};
00263 
00264   double GetFFerr_Weta2 (double pt, double eta2, int conv){
00265     if (conv) return (interpolateErrors(pt, eta2, h_c_weta2));
00266     else      return (interpolateErrors(pt, eta2, h_u_weta2));};
00267 
00268   double GetFFerr_F1    (double pt, double eta2, int conv){
00269     if (conv) return (interpolateErrors(pt, eta2, h_c_f1));
00270     else      return (interpolateErrors(pt, eta2, h_u_f1));};
00271   
00272   double GetFFerr_DE    (double pt, double eta2, int conv){
00273     if (conv) return (interpolateErrors(pt, eta2, h_c_de));
00274     else      return (interpolateErrors(pt, eta2, h_u_de));};
00275   
00276   double GetFFerr_Eratio(double pt, double eta2, int conv){
00277     if (conv) return (interpolateErrors(pt, eta2, h_c_eratio));
00278     else      return (interpolateErrors(pt, eta2, h_u_eratio));};
00279 
00280   double GetFFerr_Fside (double pt, double eta2, int conv){
00281     if (conv) return (interpolateErrors(pt, eta2, h_c_fside));
00282     else      return (interpolateErrors(pt, eta2, h_u_fside));};
00283 
00284   double GetFFerr_Wtot  (double pt, double eta2, int conv){
00285     if (conv) return (interpolateErrors(pt, eta2, h_c_wtot));
00286     else      return (interpolateErrors(pt, eta2, h_u_wtot));};
00287 
00288   double GetFFerr_W1    (double pt, double eta2, int conv){
00289     if (conv) return (interpolateErrors(pt, eta2, h_c_w1));
00290     else      return (interpolateErrors(pt, eta2, h_u_w1));};
00291 
00292   double GetFF (int var, double pt, double eta2, int conv){
00293     switch (var) {
00294     case IDVAR::RHAD1: return GetFF_Rhad1( pt, eta2, conv );
00295     case IDVAR::RHAD: return GetFF_Rhad( pt, eta2, conv );
00296     case IDVAR::E277: return GetFF_E277( pt, eta2, conv );
00297     case IDVAR::RETA: return GetFF_Reta( pt, eta2, conv );
00298     case IDVAR::RPHI: return GetFF_Rphi( pt, eta2, conv );
00299     case IDVAR::WETA2: return GetFF_Weta2( pt, eta2, conv );
00300     case IDVAR::F1: return GetFF_F1( pt, eta2, conv );
00301     case IDVAR::FSIDE: return GetFF_Fside( pt, eta2, conv );
00302     case IDVAR::WTOT: return GetFF_Wtot( pt, eta2, conv );
00303     case IDVAR::W1: return GetFF_W1( pt, eta2, conv );
00304     case IDVAR::DE: return GetFF_DE( pt, eta2, conv );
00305     case IDVAR::ERATIO: return GetFF_Eratio( pt, eta2, conv );
00306     default: return 0.0;
00307     }
00308   }
00309         
00310   // fudge a specific variable                                      
00311   double Fudge_Rhad1 ( double rhad1,  double pt, double eta2, int conv){ return ( rhad1  + GetFF_Rhad1  ( pt, eta2, conv ) ); }
00312   double Fudge_Rhad  ( double rhad,   double pt, double eta2, int conv){ return ( rhad   + GetFF_Rhad   ( pt, eta2, conv ) ); }
00313   double Fudge_E277  ( double e277,   double pt, double eta2, int conv){ return ( e277   + GetFF_E277   ( pt, eta2, conv ) ); }
00314   double Fudge_Reta  ( double reta,   double pt, double eta2, int conv){ return ( reta   + GetFF_Reta   ( pt, eta2, conv ) ); }
00315   double Fudge_Rphi  ( double rphi,   double pt, double eta2, int conv){ return ( rphi   + GetFF_Rphi   ( pt, eta2, conv ) ); }
00316   double Fudge_Weta2 ( double weta2,  double pt, double eta2, int conv){ return ( weta2  + GetFF_Weta2  ( pt, eta2, conv ) ); }
00317   double Fudge_F1    ( double f1,     double pt, double eta2, int conv){ return ( f1     + GetFF_F1     ( pt, eta2, conv ) ); }
00318   double Fudge_DE    ( double deltae, double pt, double eta2, int conv){ return ( deltae + GetFF_DE     ( pt, eta2, conv ) ); }
00319   double Fudge_Eratio( double eratio, double pt, double eta2, int conv){ return ( eratio + GetFF_Eratio ( pt, eta2, conv ) ); }
00320   double Fudge_Fside ( double fside,  double pt, double eta2, int conv){ return ( fside  + GetFF_Fside  ( pt, eta2, conv ) ); }
00321   double Fudge_Wtot  ( double wtot,   double pt, double eta2, int conv){ return ( wtot   + GetFF_Wtot   ( pt, eta2, conv ) ); }
00322   double Fudge_W1    ( double w1,     double pt, double eta2, int conv){ return ( w1     + GetFF_W1     ( pt, eta2, conv ) ); }
00323 
00324   // set shower preselection cuts
00325   // *** 2010 
00326   // 0 = tight  isolated
00327   // 1 = loose  isolated
00328   // 2 = tightx isolated
00329   // 3 = tighty isolated
00330   // 4 = tight  isolated . Distorted material 
00331   // *** 2011 
00332   // 5 = tight  isolated
00333   // 6 = tight  non-isolated
00334   // ***  
00335   // 10= tight  isolated. Old menu rel15 (tune 3 in PhotonIDTool).
00336   // ..
00337   void SetVerbose( bool verbose=true ){ m_verbose=verbose; }
00338   int  GetPreselection(){ return m_preselection; };
00339 
00340   //drawing methods
00341   TGraphErrors* GetFFmap       (int var, double eta, int isConv, int preselection);
00342 
00343   TGraphErrors* GetFFmap_Rhad1 (double eta, int isConv, int preselection);
00344   TGraphErrors* GetFFmap_Rhad  (double eta, int isConv, int preselection);
00345   TGraphErrors* GetFFmap_E277  (double eta, int isConv, int preselection);
00346   TGraphErrors* GetFFmap_Reta  (double eta, int isConv, int preselection);
00347   TGraphErrors* GetFFmap_Rphi  (double eta, int isConv, int preselection);
00348   TGraphErrors* GetFFmap_Weta2 (double eta, int isConv, int preselection);
00349   TGraphErrors* GetFFmap_F1    (double eta, int isConv, int preselection);
00350   TGraphErrors* GetFFmap_Fside (double eta, int isConv, int preselection);
00351   TGraphErrors* GetFFmap_Wtot  (double eta, int isConv, int preselection);
00352   TGraphErrors* GetFFmap_W1    (double eta, int isConv, int preselection);
00353   TGraphErrors* GetFFmap_DE    (double eta, int isConv, int preselection);
00354   TGraphErrors* GetFFmap_Eratio(double eta, int isConv, int preselection);
00355   //  TH2D* GetFFTH2D(int var, int isConv, int preselection);
00356   void LoadFFs(int preselection);
00357 
00358  private:
00359   bool m_verbose;
00360 
00361   //shower preselection to extract FFs
00362   int m_preselection;
00363 
00364   // collections of fudge factors
00365   TH2D* h_u_rhad1;
00366   TH2D* h_u_rhad;
00367   TH2D* h_u_e277;
00368   TH2D* h_u_reta;
00369   TH2D* h_u_rphi;
00370   TH2D* h_u_weta2;
00371   TH2D* h_u_f1;
00372   TH2D* h_u_fside;
00373   TH2D* h_u_wtot;
00374   TH2D* h_u_w1;
00375   TH2D* h_u_de;
00376   TH2D* h_u_eratio;
00377   TH2D* h_c_rhad1;
00378   TH2D* h_c_rhad;
00379   TH2D* h_c_e277;
00380   TH2D* h_c_reta;
00381   TH2D* h_c_rphi;
00382   TH2D* h_c_weta2;
00383   TH2D* h_c_f1;
00384   TH2D* h_c_fside;
00385   TH2D* h_c_wtot;
00386   TH2D* h_c_w1;
00387   TH2D* h_c_de;
00388   TH2D* h_c_eratio;
00389 
00390   void SafeDeleteHistograms();
00391 };
00392 
00393 #endif // #ifdef TPhotonMCShifterTool_h
00394 
00395 #ifdef TPhotonMCShifterTool_cxx
00396 
00397 TPhotonMCShifterTool::TPhotonMCShifterTool()
00398 {
00399   m_preselection = -1;
00400   h_u_rhad1 = 0;
00401   h_u_rhad = 0;
00402   h_u_e277 = 0;
00403   h_u_reta = 0;
00404   h_u_rphi = 0;
00405   h_u_weta2 = 0;
00406   h_u_f1 = 0;
00407   h_u_fside = 0;
00408   h_u_wtot = 0;
00409   h_u_w1 = 0;
00410   h_u_de = 0;
00411   h_u_eratio = 0;
00412   h_c_rhad1 = 0;
00413   h_c_rhad = 0;
00414   h_c_e277 = 0;
00415   h_c_reta = 0;
00416   h_c_rphi = 0;
00417   h_c_weta2 = 0;
00418   h_c_f1 = 0;
00419   h_c_fside = 0;
00420   h_c_wtot = 0;
00421   h_c_w1 = 0;
00422   h_c_de = 0;
00423   h_c_eratio = 0;
00424 
00425   m_verbose = false;
00426   //  std::cout << "Initiliazing the tool ---->>> preselection " << preselection << std::endl; 
00427   
00428 }
00429 
00430 TPhotonMCShifterTool::~TPhotonMCShifterTool()
00431 {
00432   SafeDeleteHistograms();
00433 }
00434 
00435 void TPhotonMCShifterTool::SafeDeleteHistograms() {
00436   if (h_u_rhad1) {delete h_u_rhad1; h_u_rhad1 = 0;}
00437   if (h_u_rhad) {delete h_u_rhad; h_u_rhad = 0;}
00438   if (h_u_e277) {delete h_u_e277; h_u_e277 = 0;}
00439   if (h_u_reta) {delete h_u_reta; h_u_reta = 0;}
00440   if (h_u_rphi) {delete h_u_rphi; h_u_rphi = 0;}
00441   if (h_u_weta2) {delete h_u_weta2; h_u_weta2 = 0;}
00442   if (h_u_f1) {delete h_u_f1; h_u_f1 = 0;}
00443   if (h_u_fside) {delete h_u_fside; h_u_fside = 0;}
00444   if (h_u_wtot) {delete h_u_wtot; h_u_wtot = 0;}
00445   if (h_u_w1) {delete h_u_w1; h_u_w1 = 0;}
00446   if (h_u_de) {delete h_u_de; h_u_de = 0;}
00447   if (h_u_eratio) {delete h_u_eratio; h_u_eratio = 0;}
00448   if (h_c_rhad1) {delete h_c_rhad1; h_c_rhad1 = 0;}
00449   if (h_c_rhad) {delete h_c_rhad; h_c_rhad = 0;}
00450   if (h_c_e277) {delete h_c_e277; h_c_e277 = 0;}
00451   if (h_c_reta) {delete h_c_reta; h_c_reta = 0;}
00452   if (h_c_rphi) {delete h_c_rphi; h_c_rphi = 0;}
00453   if (h_c_weta2) {delete h_c_weta2; h_c_weta2 = 0;}
00454   if (h_c_f1) {delete h_c_f1; h_c_f1 = 0;}
00455   if (h_c_fside) {delete h_c_fside; h_c_fside = 0;}
00456   if (h_c_wtot) {delete h_c_wtot; h_c_wtot = 0;}
00457   if (h_c_w1) {delete h_c_w1; h_c_w1 = 0;}
00458   if (h_c_de) {delete h_c_de; h_c_de = 0;}
00459   if (h_c_eratio) {delete h_c_eratio; h_c_eratio = 0;}
00460 }
00461 
00462 #endif// #ifdef TPhotonMCShifterTool_cxx
00463