PWGLF/totEt/macros/hadEt/GetCorrections.C

   1 //Create by Christine Nattrass, Rebecca Scott, Irakli Martashvili
   2 //University of Tennessee at Knoxville
   3 //This macro takes an input file created by AliAnalysisHadEtMonteCarlo and creates an AliAnalysisHadEtCorrections for the determination of the corrected Et
   4 // #include "TFile.h"
   5 // #include <iostream>
   6 // #include "AliAnalysisHadEtCorrections.h"
   7
   8 // #include <iostream>
   9 // #include <TROOT.h>
  10 // #include <TSystem.h>
  11 // #include "TStopwatch.h"
  12
  13 Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool ispp = true, bool forSim = true, bool TPC, bool hadronic = false, float etacut = 0.7);
  14 TH1D *GetHistoCorrNeutral(float cut, char *name, bool ispp, bool forSim, int mycase, bool eta, int color, int marker, bool hadronic = false);
  15
  16 Float_t CorrPtCut(float ptcut, char *prodname = "Enter Production Name", char *shortprodname = "EnterProductionName", bool ispp = true, bool forSim = true, int mycase = 0);
  17 TH1D *GetHistoCorrPtCut(float ptcut = 0.15, char *name, bool ispp = true, bool forSim = true, int mycase);
  18
  19 TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, bool eta, bool ispp = true, bool forSim = true);
  20 TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, bool ispp = true, bool forSim = true);
  21 Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, bool ispp, bool forSim);
  22
  23 TH1D *GetHistoNoID(float etacut, char *name, bool eta, bool TPC, bool ispp, bool forSim);
  24 TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, bool ispp, bool forSim);
  25 Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, bool ispp, bool forSim);
  26
  27 TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name);
  28 TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC, bool ITS, int cb = -1, int cblast = -1);
  29 void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname);
  30
  31 TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC,bool ispp,bool forSim);
  32 void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC,bool ispp,bool forSim);
  33
  34 //Some variables that we'll use multiple times.  We'll declare them here since they don't seem to delete right in the functions
  35 char prefix[100];
  36 char histoname[100];
  37 char epsname[100];
  38 char pngname[100];
  39 TFile *file = NULL;//initiated in main function
  40 const char *mynameTPC = "TPC";
  41 const char *mynameITS = "ITS";
  42 const char *mynameTPCITS = "TPCITS";
  43 const char *detectorEMCAL = "EMCAL";
  44 const char *detectorPHOS = "PHOS";
  45 const char *reweightedNo = "";
  46 const char *reweightedYes = "Reweighted";
  47
  48 //===========================================================================================
  49
  50 void GetCorrections(char *prodname = "Enter Production Name", char *shortprodname = "EnterProductionName", bool ispp = true, bool forSim = true, bool TPC = true, char *infilename="Et.ESD.new.sim.merged.root", int dataset = 2009){
  51     TStopwatch timer;
  52     timer.Start();
  53     gSystem->Load("libTree.so");
  54     gSystem->Load("libGeom.so");
  55     gSystem->Load("libVMC.so");
  56     gSystem->Load("libXMLIO.so");
  57
  58     gSystem->Load("libSTEERBase.so");
  59     gSystem->Load("libESD.so");
  60     gSystem->Load("libAOD.so");
  61
  62     gSystem->Load("libANALYSIS");
  63     gSystem->Load("libANALYSISalice");
  64
  65     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  66    gROOT->ProcessLine(".L AliAnalysisEtCuts.cxx+g");
  67    gROOT->ProcessLine(".L AliAnalysisHadEtCorrections.cxx+g");
  68    file = new TFile(infilename);
  69
  70    char outfilename[200];
  71    char *sim = "ForData";
  72    if(forSim) sim = "ForSimulations";
  73    char *system = "PbPb";
  74    if(ispp) system = "pp";
  75    sprintf(outfilename,"rootFiles/corrections/corrections.%s.%s.%s.root",shortprodname,system,sim);
  76    TFile *outfile = new TFile(outfilename,"RECREATE");
  77    AliAnalysisHadEtCorrections *hadCorrectionEMCAL = new AliAnalysisHadEtCorrections();
  78    hadCorrectionEMCAL->SetName("hadCorrectionEMCAL");
  79    float etacut = 0.7;
  80    hadCorrectionEMCAL->SetEtaCut(etacut);
  81    hadCorrectionEMCAL->IsData(!forSim);
  82    hadCorrectionEMCAL->IsEMCal(kTRUE);
  83    hadCorrectionEMCAL->SetProduction(shortprodname);
  84    hadCorrectionEMCAL->SetProductionDescription(prodname);
  85    hadCorrectionEMCAL->SetDataSet(dataset);
  86    //float etacut = hadCorrectionEMCAL->GetEtaCut();
  87    //cout<<"eta cut is "<<etacut<<endl;
  88    hadCorrectionEMCAL->SetAcceptanceCorrectionFull(1.0);
  89    cout<<"Warning:  Acceptance corrections will have to be updated to include real acceptance maps of the EMCAL"<<endl;
  90    hadCorrectionEMCAL->SetAcceptanceCorrectionPHOS(360.0/60.0);
  91    hadCorrectionEMCAL->SetAcceptanceCorrectionEMCAL(360.0/60.0);
  92
  93    float ptcut = 0.1;
  94    float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,false,etacut);
  95    hadCorrectionEMCAL->SetNeutralCorrection(neutralCorr);
  96    //Using error from data, see analysis note for details
  97    if(ispp){
  98      hadCorrectionEMCAL->SetNeutralCorrectionLowBound(neutralCorr*(1.0-0.013));
  99      hadCorrectionEMCAL->SetNeutralCorrectionHighBound(neutralCorr*(1.0+0.013));
 100    }
 101    else{
 102      hadCorrectionEMCAL->SetNeutralCorrectionLowBound(neutralCorr*(1.0-0.049));
 103      hadCorrectionEMCAL->SetNeutralCorrectionHighBound(neutralCorr*(1.0+0.049));
 104    }
 105
 106    float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,true,etacut);
 107    hadCorrectionEMCAL->SetNotHadronicCorrection(hadronicCorr);
 108    if(ispp){
 109      hadCorrectionEMCAL->SetNotHadronicCorrectionLowBound(hadronicCorr*(1.0-0.008));
 110      hadCorrectionEMCAL->SetNotHadronicCorrectionHighBound(hadronicCorr*(1.0+0.008));
 111    }
 112    else{
 113      hadCorrectionEMCAL->SetNotHadronicCorrectionLowBound(hadronicCorr*(1.0-0.023));
 114      hadCorrectionEMCAL->SetNotHadronicCorrectionHighBound(hadronicCorr*(1.0+0.023));
 115    }
 116
 117    float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim);
 118    hadCorrectionEMCAL->SetpTCutCorrectionITS(ptcutITS);
 119    float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim);
 120    hadCorrectionEMCAL->SetpTCutCorrectionTPC(ptcutTPC);
 121    float ptcutITSLow = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim,-1);
 122    float ptcutTPCLow = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim,-1);
 123    hadCorrectionEMCAL->SetpTCutCorrectionITSLowBound(ptcutITSLow);
 124    hadCorrectionEMCAL->SetpTCutCorrectionTPCLowBound(ptcutTPCLow);
 125    float ptcutITSHigh = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim,1);
 126    float ptcutTPCHigh = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim,1);
 127    hadCorrectionEMCAL->SetpTCutCorrectionITSHighBound(ptcutITSHigh);
 128    hadCorrectionEMCAL->SetpTCutCorrectionTPCHighBound(ptcutTPCHigh);
 129
 130    TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDEMCALTPC",prodname,shortprodname,true,ispp,forSim);
 131    TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDEMCALITS",prodname,shortprodname,false,ispp,forSim);
 132    hadCorrectionEMCAL->SetNotIDCorrectionTPC(NotIDTPC);
 133    hadCorrectionEMCAL->SetNotIDCorrectionITS(NotIDITS);
 134
 135    Float_t NotIDConstTPC = CorrNotIDConst(0.15,etacut,"CorrNotIDEMCALTPC2",prodname,shortprodname,true,ispp,forSim);
 136    Float_t NotIDConstITS = CorrNotIDConst(0.10,etacut,"CorrNotIDEMCALTPC2",prodname,shortprodname,true,ispp,forSim);
 137    hadCorrectionEMCAL->SetNotIDConstCorrectionTPC(1.0/NotIDConstTPC);
 138    hadCorrectionEMCAL->SetNotIDConstCorrectionITS(1.0/NotIDConstITS);
 139    if(ispp){
 140      hadCorrectionEMCAL->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*(1.0-0.010));
 141      hadCorrectionEMCAL->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*(1.0-0.010));
 142      hadCorrectionEMCAL->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*(1.0+0.010));
 143      hadCorrectionEMCAL->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*(1.0+0.010));
 144    }
 145    else{
 146      hadCorrectionEMCAL->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*(1.0-0.022));
 147      hadCorrectionEMCAL->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*(1.0-0.022));
 148      hadCorrectionEMCAL->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*(1.0+0.022));
 149      hadCorrectionEMCAL->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*(1.0+0.022));
 150    }
 151
 152    TH1D *NoID = CorrNoID(etacut,"CorrNoIDEMCAL",prodname,shortprodname,ispp,forSim);
 153    hadCorrectionEMCAL->SetNotIDCorrectionNoPID(NoID);
 154
 155    Float_t NoIDTPC = CorrNoIDConst(etacut,0.15,"CorrNoIDEMCAL2",prodname,shortprodname,ispp,forSim);
 156    Float_t NoIDITS = CorrNoIDConst(etacut,0.1,"CorrNoIDEMCAL2",prodname,shortprodname,ispp,forSim);
 157    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoID(1./NoIDTPC);
 158    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoID(1./NoIDITS);
 159    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoIDLowBound(1./NoIDTPC*.98);
 160    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoIDLowBound(1./NoIDITS*.98);
 161    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoIDHighBound(1./NoIDTPC*1.02);
 162    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoIDHighBound(1./NoIDITS*1.02);
 163
 164    TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,true);
 165    hadCorrectionEMCAL->SetEfficiencyPionTPC(efficiencyPionTPC);
 166    TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,true);
 167    hadCorrectionEMCAL->SetEfficiencyKaonTPC(efficiencyKaonTPC);
 168    TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,true);
 169    hadCorrectionEMCAL->SetEfficiencyProtonTPC(efficiencyProtonTPC);
 170    TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,true);
 171    hadCorrectionEMCAL->SetEfficiencyHadronTPC(efficiencyHadronTPC);
 172    TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,true);
 173    hadCorrectionEMCAL->SetEfficiencyPionITS(efficiencyPionITS);
 174    TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,true);
 175    hadCorrectionEMCAL->SetEfficiencyKaonITS(efficiencyKaonITS);
 176    TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,true);
 177    hadCorrectionEMCAL->SetEfficiencyProtonITS(efficiencyProtonITS);
 178    TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,true);
 179    hadCorrectionEMCAL->SetEfficiencyHadronITS(efficiencyHadronITS);
 180
 181    if(!ispp){
 182      TH1D *efficiencyPionTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB0",1,1,20,true,true,0,4);
 183      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB0->Clone(Form("Test%i",i)),i);
 184      TH1D *efficiencyPionTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB5",1,1,20,true,true,5,9);
 185      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB5->Clone(Form("Test%i",i)),i);
 186      TH1D *efficiencyPionTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB10",1,1,20,true,true,10,15);
 187      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB10->Clone(Form("Test%i",i)),i);
 188
 189      TH1D *efficiencyKaonTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB0",2,1,20,true,true,0,4);
 190      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB0->Clone(Form("Test%i",i)),i);
 191      TH1D *efficiencyKaonTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB5",2,1,20,true,true,5,9);
 192      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB5->Clone(Form("Test%i",i)),i);
 193      TH1D *efficiencyKaonTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB10",2,1,20,true,true,10,15);
 194      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB10->Clone(Form("Test%i",i)),i);//Kaon
 195
 196      TH1D *efficiencyProtonTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB0",3,1,20,true,true,0,4);
 197      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB0->Clone(Form("Test%i",i)),i);
 198      TH1D *efficiencyProtonTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB5",3,1,20,true,true,5,9);
 199      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB5->Clone(Form("Test%i",i)),i);
 200      TH1D *efficiencyProtonTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB10",3,1,20,true,true,10,15);
 201      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB10->Clone(Form("Test%i",i)),i);//Proton
 202
 203      TH1D *efficiencyHadronTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB0",0,1,20,true,true,0,4);
 204      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB0->Clone(Form("Test%i",i)),i);
 205      TH1D *efficiencyHadronTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB5",0,1,20,true,true,5,9);
 206      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB5->Clone(Form("Test%i",i)),i);
 207      TH1D *efficiencyHadronTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB10",0,1,20,true,true,10,15);
 208      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB10->Clone(Form("Test%i",i)),i);//Hadron
 209
 210
 211      TH1D *efficiencyPionITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB0",1,1,20,false,true,0,4);
 212      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB0->Clone(Form("Test%i",i)),i);
 213      TH1D *efficiencyPionITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB5",1,1,20,false,true,5,9);
 214      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB5->Clone(Form("Test%i",i)),i);
 215      TH1D *efficiencyPionITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB10",1,1,20,false,true,10,15);
 216      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB10->Clone(Form("Test%i",i)),i);//Pion
 217
 218      TH1D *efficiencyKaonITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB0",2,1,20,false,true,0,4);
 219      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB0->Clone(Form("Test%i",i)),i);
 220      TH1D *efficiencyKaonITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB5",2,1,20,false,true,5,9);
 221      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB5->Clone(Form("Test%i",i)),i);
 222      TH1D *efficiencyKaonITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB10",2,1,20,false,true,10,15);
 223      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB10->Clone(Form("Test%i",i)),i);//Kaon
 224
 225      TH1D *efficiencyProtonITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB0",3,1,20,false,true,0,4);
 226      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB0->Clone(Form("Test%i",i)),i);
 227      TH1D *efficiencyProtonITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB5",3,1,20,false,true,5,9);
 228      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB5->Clone(Form("Test%i",i)),i);
 229      TH1D *efficiencyProtonITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB10",3,1,20,false,true,10,15);
 230      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB10->Clone(Form("Test%i",i)),i);//Proton
 231
 232      TH1D *efficiencyHadronITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB0",0,1,20,false,true,0,4);
 233      for(int i=0;i<=4;i++) hadCorrectionEMCAL->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB0->Clone(Form("Test%i",i)),i);
 234      TH1D *efficiencyHadronITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB5",0,1,20,false,true,5,9);
 235      for(int i=5;i<=9;i++) hadCorrectionEMCAL->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB5->Clone(Form("Test%i",i)),i);
 236      TH1D *efficiencyHadronITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB10",0,1,20,false,true,10,15);
 237      for(int i=10;i<=19;i++) hadCorrectionEMCAL->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB10->Clone(Form("Test%i",i)),i);//Hadron
 238    }//EMCAL
 239
 240    //CorrEfficiencyPlots(true,prodname,shortprodname);
 241    //CorrEfficiencyPlots(false,prodname,shortprodname,infilename);
 242
 243
 244    //hadCorrectionEMCAL->GetEfficiencyHadronTPC()->Draw();
 245    TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,ispp,forSim);
 246    TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,ispp,forSim);
 247    hadCorrectionEMCAL->SetBackgroundCorrectionTPC(backgroundTPC);
 248    hadCorrectionEMCAL->SetBackgroundCorrectionITS(backgroundITS);
 249    //CorrBkgdPlots(prodname,shortprodname,true,ispp,forSim);
 250    //CorrBkgdPlots(prodname,shortprodname,false,ispp,forSim);
 251
 252    hadCorrectionEMCAL->Report();
 253
 254    outfile->cd();
 255    hadCorrectionEMCAL->Write();
 256    outfile->Write();
 257    delete hadCorrectionEMCAL;
 258
 259    AliAnalysisHadEtCorrections *hadCorrectionPHOS = new AliAnalysisHadEtCorrections();
 260    hadCorrectionPHOS->SetName("hadCorrectionPHOS");
 261    float etacut = 0.12;
 262    hadCorrectionPHOS->SetEtaCut(etacut);
 263    hadCorrectionPHOS->IsData(!forSim);
 264    hadCorrectionPHOS->IsEMCal(kTRUE);
 265    hadCorrectionPHOS->SetProduction(shortprodname);
 266    hadCorrectionPHOS->SetProductionDescription(prodname);
 267    hadCorrectionPHOS->SetDataSet(dataset);
 268    //float etacut = hadCorrectionPHOS->GetEtaCut();
 269    //cout<<"eta cut is "<<etacut<<endl;
 270    hadCorrectionPHOS->SetAcceptanceCorrectionFull(1.0);
 271    cout<<"Warning:  Acceptance corrections will have to be updated to include real acceptance maps of the PHOS"<<endl;
 272    hadCorrectionPHOS->SetAcceptanceCorrectionPHOS(360.0/60.0);
 273    hadCorrectionPHOS->SetAcceptanceCorrectionEMCAL(360.0/60.0);
 274
 275    float ptcut = 0.1;
 276    float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,false,etacut);
 277    hadCorrectionPHOS->SetNeutralCorrection(neutralCorr);
 278    //Using error from data, see analysis note for details
 279    if(ispp){
 280      hadCorrectionPHOS->SetNeutralCorrectionLowBound(neutralCorr*(1.0-.013));
 281      hadCorrectionPHOS->SetNeutralCorrectionHighBound(neutralCorr*(1.0+.013));
 282    }
 283    else{
 284      hadCorrectionPHOS->SetNeutralCorrectionLowBound(neutralCorr*(1.0-0.049));
 285      hadCorrectionPHOS->SetNeutralCorrectionHighBound(neutralCorr*(1.0+0.049));
 286    }
 287
 288    float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,true,etacut);
 289    hadCorrectionPHOS->SetNotHadronicCorrection(hadronicCorr);
 290    if(ispp){
 291      hadCorrectionPHOS->SetNotHadronicCorrectionLowBound(neutralCorr*(1.0-0.008));
 292      hadCorrectionPHOS->SetNotHadronicCorrectionHighBound(neutralCorr*(1.0+0.008));
 293    }
 294    else{
 295      hadCorrectionPHOS->SetNotHadronicCorrectionLowBound(neutralCorr*(1.0-0.023));
 296      hadCorrectionPHOS->SetNotHadronicCorrectionHighBound(neutralCorr*(1.0+0.023));
 297    }
 298
 299    float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim);
 300    hadCorrectionPHOS->SetpTCutCorrectionITS(ptcutITS);
 301    float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim);
 302    hadCorrectionPHOS->SetpTCutCorrectionTPC(ptcutTPC);
 303
 304    float ptcutITSLow = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim,-1);
 305    float ptcutTPCLow = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim,-1);
 306    hadCorrectionPHOS->SetpTCutCorrectionITSLowBound(ptcutITSLow);
 307    hadCorrectionPHOS->SetpTCutCorrectionTPCLowBound(ptcutTPCLow);
 308    float ptcutITSHigh = CorrPtCut(0.1,prodname,shortprodname,ispp,forSim,1);
 309    float ptcutTPCHigh = CorrPtCut(0.15,prodname,shortprodname,ispp,forSim,1);
 310    hadCorrectionPHOS->SetpTCutCorrectionITSHighBound(ptcutITSHigh);
 311    hadCorrectionPHOS->SetpTCutCorrectionTPCHighBound(ptcutTPCHigh);
 312
 313    TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDPHOSTPC",prodname,shortprodname,true,ispp,forSim);
 314    TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDPHOSITS",prodname,shortprodname,false,ispp,forSim);
 315    hadCorrectionPHOS->SetNotIDCorrectionTPC(NotIDTPC);
 316    hadCorrectionPHOS->SetNotIDCorrectionITS(NotIDITS);
 317
 318    Float_t NotIDConstTPC = CorrNotIDConst(0.15,etacut,"CorrNotIDPHOSTPC2",prodname,shortprodname,true,ispp,forSim);
 319    Float_t NotIDConstITS = CorrNotIDConst(0.10,etacut,"CorrNotIDPHOSTPC2",prodname,shortprodname,true,ispp,forSim);
 320    hadCorrectionPHOS->SetNotIDConstCorrectionTPC(1./NotIDConstTPC);
 321    hadCorrectionPHOS->SetNotIDConstCorrectionITS(1./NotIDConstITS);
 322    if(ispp){
 323      hadCorrectionPHOS->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*(1.0-0.010));
 324      hadCorrectionPHOS->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*(1.0-0.010));
 325      hadCorrectionPHOS->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*(1.0+0.010));
 326      hadCorrectionPHOS->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*(1.0+0.010));
 327    }
 328    else{
 329      hadCorrectionPHOS->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*(1.0-0.022));
 330      hadCorrectionPHOS->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*(1.0-0.022));
 331      hadCorrectionPHOS->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*(1.0+0.022));
 332      hadCorrectionPHOS->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*(1.0+0.022));
 333    }
 334
 335
 336    TH1D *NoID = CorrNoID(etacut,"CorrNoIDPHOS",prodname,shortprodname,ispp,forSim);
 337    hadCorrectionPHOS->SetNotIDCorrectionNoPID(NoID);
 338
 339
 340    Float_t NoIDTPC = CorrNoIDConst(etacut,0.15,"CorrNoIDPHOS2",prodname,shortprodname,ispp,forSim);
 341    Float_t NoIDITS = CorrNoIDConst(etacut,0.1,"CorrNoIDPHOS2",prodname,shortprodname,ispp,forSim);
 342    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoID(1./NoIDTPC);
 343    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoID(1./NoIDITS);
 344    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoIDLowBound(1./NoIDTPC*.98);
 345    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoIDLowBound(1./NoIDITS*.98);
 346    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoIDHighBound(1./NoIDTPC*1.02);
 347    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoIDHighBound(1./NoIDITS*1.02);
 348
 349    TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,true);
 350    TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,true);
 351    TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,true);
 352    TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,true);
 353    TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,true);
 354    TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,true);
 355    TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,true);
 356    TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,true);
 357    //CorrEfficiencyPlots(true,prodname,shortprodname);
 358    //CorrEfficiencyPlots(false,prodname,shortprodname);
 359    hadCorrectionPHOS->SetEfficiencyPionTPC(efficiencyPionTPC);
 360    hadCorrectionPHOS->SetEfficiencyKaonTPC(efficiencyKaonTPC);
 361    hadCorrectionPHOS->SetEfficiencyProtonTPC(efficiencyProtonTPC);
 362    hadCorrectionPHOS->SetEfficiencyHadronTPC(efficiencyHadronTPC);
 363    hadCorrectionPHOS->SetEfficiencyPionITS(efficiencyPionITS);
 364    hadCorrectionPHOS->SetEfficiencyKaonITS(efficiencyKaonITS);
 365    hadCorrectionPHOS->SetEfficiencyProtonITS(efficiencyProtonITS);
 366    hadCorrectionPHOS->SetEfficiencyHadronITS(efficiencyHadronITS);
 367
 368
 369    if(!ispp){
 370      TH1D *efficiencyPionTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB0",1,1,20,true,true,0,4);
 371      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB0->Clone(Form("Test%i",i)),i);
 372      TH1D *efficiencyPionTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB5",1,1,20,true,true,5,9);
 373      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB5->Clone(Form("Test%i",i)),i);
 374      TH1D *efficiencyPionTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyPionTPCCB10",1,1,20,true,true,10,15);
 375      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyPionTPC((TH1D*)efficiencyPionTPCCB10->Clone(Form("Test%i",i)),i);
 376
 377      TH1D *efficiencyKaonTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB0",2,1,20,true,true,0,4);
 378      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB0->Clone(Form("Test%i",i)),i);
 379      TH1D *efficiencyKaonTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB5",2,1,20,true,true,5,9);
 380      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB5->Clone(Form("Test%i",i)),i);
 381      TH1D *efficiencyKaonTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyKaonTPCCB10",2,1,20,true,true,10,15);
 382      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyKaonTPC((TH1D*)efficiencyKaonTPCCB10->Clone(Form("Test%i",i)),i);//Kaon
 383
 384      TH1D *efficiencyProtonTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB0",3,1,20,true,true,0,4);
 385      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB0->Clone(Form("Test%i",i)),i);
 386      TH1D *efficiencyProtonTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB5",3,1,20,true,true,5,9);
 387      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB5->Clone(Form("Test%i",i)),i);
 388      TH1D *efficiencyProtonTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyProtonTPCCB10",3,1,20,true,true,10,15);
 389      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyProtonTPC((TH1D*)efficiencyProtonTPCCB10->Clone(Form("Test%i",i)),i);//Proton
 390
 391      TH1D *efficiencyHadronTPCCB0 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB0",0,1,20,true,true,0,4);
 392      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB0->Clone(Form("Test%i",i)),i);
 393      TH1D *efficiencyHadronTPCCB5 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB5",0,1,20,true,true,5,9);
 394      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB5->Clone(Form("Test%i",i)),i);
 395      TH1D *efficiencyHadronTPCCB10 = GetHistoEfficiency(etacut,"hEfficiencyHadronTPCCB10",0,1,20,true,true,10,15);
 396      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyHadronTPC((TH1D*)efficiencyHadronTPCCB10->Clone(Form("Test%i",i)),i);//Hadron
 397
 398
 399      TH1D *efficiencyPionITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB0",1,1,20,false,true,0,4);
 400      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB0->Clone(Form("Test%i",i)),i);
 401      TH1D *efficiencyPionITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB5",1,1,20,false,true,5,9);
 402      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB5->Clone(Form("Test%i",i)),i);
 403      TH1D *efficiencyPionITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyPionITSCB10",1,1,20,false,true,10,15);
 404      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyPionITS((TH1D*)efficiencyPionITSCB10->Clone(Form("Test%i",i)),i);//Pion
 405
 406      TH1D *efficiencyKaonITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB0",2,1,20,false,true,0,4);
 407      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB0->Clone(Form("Test%i",i)),i);
 408      TH1D *efficiencyKaonITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB5",2,1,20,false,true,5,9);
 409      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB5->Clone(Form("Test%i",i)),i);
 410      TH1D *efficiencyKaonITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyKaonITSCB10",2,1,20,false,true,10,15);
 411      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyKaonITS((TH1D*)efficiencyKaonITSCB10->Clone(Form("Test%i",i)),i);//Kaon
 412
 413      TH1D *efficiencyProtonITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB0",3,1,20,false,true,0,4);
 414      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB0->Clone(Form("Test%i",i)),i);
 415      TH1D *efficiencyProtonITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB5",3,1,20,false,true,5,9);
 416      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB5->Clone(Form("Test%i",i)),i);
 417      TH1D *efficiencyProtonITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyProtonITSCB10",3,1,20,false,true,10,15);
 418      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyProtonITS((TH1D*)efficiencyProtonITSCB10->Clone(Form("Test%i",i)),i);//Proton
 419
 420      TH1D *efficiencyHadronITSCB0 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB0",0,1,20,false,true,0,4);
 421      for(int i=0;i<=4;i++) hadCorrectionPHOS->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB0->Clone(Form("Test%i",i)),i);
 422      TH1D *efficiencyHadronITSCB5 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB5",0,1,20,false,true,5,9);
 423      for(int i=5;i<=9;i++) hadCorrectionPHOS->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB5->Clone(Form("Test%i",i)),i);
 424      TH1D *efficiencyHadronITSCB10 = GetHistoEfficiency(etacut,"hEfficiencyHadronITSCB10",0,1,20,false,true,10,15);
 425      for(int i=10;i<=19;i++) hadCorrectionPHOS->SetEfficiencyHadronITS((TH1D*)efficiencyHadronITSCB10->Clone(Form("Test%i",i)),i);//Hadron
 426    }//EMCAL
 427
 428    TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,ispp,forSim);
 429    TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,ispp,forSim);
 430    hadCorrectionPHOS->SetBackgroundCorrectionTPC(backgroundTPC);
 431    hadCorrectionPHOS->SetBackgroundCorrectionITS(backgroundITS);
 432    CorrBkgdPlots(prodname,shortprodname,true,ispp,forSim);
 433    CorrBkgdPlots(prodname,shortprodname,false,ispp,forSim);
 434
 435    hadCorrectionPHOS->Report();
 436    //Write the output
 437    outfile->cd();
 438    hadCorrectionPHOS->Write();
 439    outfile->Write();
 440    outfile->Close();
 441
 442
 443   timer.Stop();
 444   timer.Print();
 445 }
 446
 447 //==================================CorrNeutral==============================================
 448 Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool ispp, bool forSim, bool TPC, bool hadronic, float etacut){
 449   if(!forSim){//for data we have evaluated the neutral correction from ALICE data
 450     if(hadronic){//for tot et from had et
 451       if(ispp){
 452         return 1.0/0.571;
 453       }
 454       else{
 455         return 1.0/0.549;
 456       }
 457     }
 458     else{//for had et only
 459       if(ispp){
 460         return 1.0/0.736;
 461       }
 462       else{
 463         return 1.0/0.689;
 464       }
 465     }
 466   }
 467   gStyle->SetOptTitle(0);
 468   gStyle->SetOptStat(0);
 469   gStyle->SetOptFit(0);
 470   TCanvas *c = new TCanvas("c","c",400,400);
 471   c->SetTopMargin(0.0);
 472   c->SetRightMargin(0.0);
 473   c->SetBorderSize(0);
 474   c->SetFillColor(0);
 475   c->SetFillColor(0);
 476   c->SetBorderMode(0);
 477   c->SetFrameFillColor(0);
 478   c->SetFrameBorderMode(0);
 479
 480   TPad *ptpad = c->cd(1);
 481   ptpad->SetTopMargin(0.04);
 482   ptpad->SetRightMargin(0.04);
 483   ptpad->SetLeftMargin(0.149288);
 484   ptpad->SetBorderSize(0);
 485   ptpad->SetFillColor(0);
 486   ptpad->SetFillColor(0);
 487   ptpad->SetBorderMode(0);
 488   ptpad->SetFrameFillColor(0);
 489   ptpad->SetFrameBorderMode(0);
 490
 491   int phosmarker = 20;
 492
 493   sprintf(prefix,"%s%2.1f",shortprodname,ptcut);
 494
 495   sprintf(histoname,"%stotal",histoname);
 496   int colortotal = 1;
 497   int casetotal = 4;
 498   if(hadronic) casetotal = 8;
 499   TH1D *total = GetHistoCorrNeutral(ptcut,histoname,ispp,forSim,casetotal,false,colortotal,phosmarker,hadronic);
 500
 501   int colorallneutral = 2;
 502   TH1D *allneutral = GetHistoCorrNeutral(ptcut,"allneutral",ispp,forSim,3,false,colorallneutral,phosmarker,hadronic);
 503
 504   int colorchargedsecondary = TColor::kViolet-3;
 505   TH1D *chargedsecondary = GetHistoCorrNeutral(ptcut,"chargedsecondary",ispp,forSim,2,false,colorchargedsecondary,phosmarker,hadronic);
 506
 507   int colorneutralUndet = 4;
 508   TH1D *neutralUndet = GetHistoCorrNeutral(ptcut,"neutralUndet",ispp,forSim,1,false,colorneutralUndet,phosmarker,hadronic);
 509
 510   int colorv0 = TColor::kGreen+2;
 511   TH1D *v0 = GetHistoCorrNeutral(ptcut,"v0",ispp,forSim,0,false,colorv0,phosmarker,hadronic);
 512
 513   int colorem = TColor::kCyan;
 514   TH1D *em = GetHistoCorrNeutral(ptcut,"em",ispp,forSim,9,false,colorem,phosmarker,hadronic);
 515
 516   TF1 *func = new TF1("func","[0]",-.7,.7);
 517   func->SetParameter(0,0.2);
 518   total->Fit(func,"","",-etacut,etacut);
 519
 520   //total->SetAxisRange(0.0,4);
 521   total->GetXaxis()->SetLabelSize(0.05);
 522   total->GetYaxis()->SetLabelSize(0.045);
 523   total->GetXaxis()->SetTitleSize(0.05);
 524   total->GetYaxis()->SetTitleSize(0.06);
 525   if(hadronic){
 526     total->SetMaximum(0.6);
 527     total->SetMinimum(0.0);
 528   }
 529   else{
 530     total->SetMaximum(0.3);
 531     total->SetMinimum(0.0);
 532   }
 533   total->Draw();
 534   allneutral->Draw("same");
 535   chargedsecondary->Draw("same");
 536   neutralUndet->Draw("same");
 537   v0->Draw("same");
 538   if(hadronic) em->Draw("same");
 539
 540   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 541   tex->SetTextSize(0.0537634);
 542   tex->Draw();
 543   TLegend *leg2 = new TLegend(0.518321,0.746873,0.774812,0.955343);
 544   leg2->AddEntry(total,"Total");
 545   leg2->AddEntry(allneutral,"#Lambda,#bar{#Lambda},K^{0}_{S},K^{0}_{L},n,#bar{n}");
 546   leg2->AddEntry(neutralUndet,"K^{0}_{L},n,#bar{n}");
 547   if(hadronic) leg2->AddEntry(em,"e^{#pm},#gamma,#eta,#pi^{0},#omega");
 548   leg2->SetFillStyle(0);
 549   leg2->SetFillColor(0);
 550   leg2->SetBorderSize(0);
 551   leg2->SetTextSize(0.0548607);
 552   leg2->Draw();
 553   if(hadronic){
 554     sprintf(epsname,"pics/%s/fhadronic.eps",shortprodname);
 555     sprintf(pngname,"pics/%s/fhadronic.png",shortprodname);
 556   }
 557   else{
 558     sprintf(epsname,"pics/%s/fneutral.eps",shortprodname);
 559     sprintf(pngname,"pics/%s/fneutral.png",shortprodname);
 560   }
 561   c->SaveAs(epsname);
 562   c->SaveAs(pngname);
 563
 564   delete total;
 565   delete allneutral;
 566   delete chargedsecondary;
 567   delete neutralUndet;
 568   delete v0;
 569   delete em;
 570   delete c;
 571   float corr = func->GetParameter(0);
 572   delete func;
 573   delete tex;
 574   delete leg2;
 575   return 1.0/(1.0-corr);
 576
 577 }
 578 TH1D *GetHistoCorrNeutral(float cut, char *name, bool ispp, bool forSim, int mycase, bool eta, int color, int marker, bool hadronic){
 579   file->cd();
 580   char *reweightname = "";
 581   if(!forSim) reweightname = "Reweighted";
 582   TH2F *numeratorParent;
 583   switch(mycase){
 584   case 0:
 585     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("v0");
 586     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 587     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 588     break;
 589   case 1:
 590     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)))->Clone("Knnbar");
 591     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 592     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 593     break;
 594   case 2:
 595     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedOmega"))->Clone("ch2ndary");
 596     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 597     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 598     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 599     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 600     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 601     break;
 602   case 3:
 603     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 604     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 605     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 606     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 607     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 608     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 609     break;
 610   case 4:
 611     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 612     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 613     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 614     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 615     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 616     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 617     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 618     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 619     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 620     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 621     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 622     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 623     break;
 624   case 5:
 625     numeratorParent= (TH2F*)((TH2F*) out2->FindObject("EtSimulatedXi"))->Clone("allxi");
 626     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 627     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 628     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 629     break;
 630   case 6:
 631     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedOmega"))->Clone("allomega");
 632     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 633     break;
 634   case 7:
 635     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedSigma"))->Clone("allsigma");
 636     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiSigma"));
 637     break;
 638   case 8:
 639     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 640     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 641     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 642     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 643     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 644     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 645     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 646     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 647     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 648     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 649     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 650     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 651     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedGamma"));
 652     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 653     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 654     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 655     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 656     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 657     break;
 658   case 9:
 659     numeratorParent= (TH2F*)((TH2F*) out2->FindObject("EtSimulatedGamma"))->Clone("allem");
 660     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 661     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 662     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 663     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 664     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 665     break;
 666   }
 667
 668   TH2F *allhad;
 669   //allhad=(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedAllHadron"))->Clone("id");
 670   allhad=(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedPiPlus"))->Clone("id");
 671   allhad->Add((TH2F*) out2->FindObject("EtSimulatedPiMinus"));
 672   allhad->Add((TH2F*) out2->FindObject("EtSimulatedKMinus"));
 673   allhad->Add((TH2F*) out2->FindObject("EtSimulatedKPlus"));
 674   allhad->Add((TH2F*) out2->FindObject("EtSimulatedProton"));
 675   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiProton"));
 676   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)));
 677   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 678   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 679   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 680   allhad->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 681   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 682   allhad->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 683   allhad->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 684   allhad->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 685   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 686   allhad->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 687   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 688   allhad->Add((TH2F*) out2->FindObject("EtSimulatedSigma"));
 689   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiSigma"));
 690   allhad->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 691   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 692
 693   if(hadronic){//if we are getting the correction for the hadronic only case...
 694     allhad->Add((TH2F*) out2->FindObject("EtSimulatedGamma"));
 695     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 696     allhad->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 697     allhad->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 698     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 699     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 700   }
 701
 702   TH1D *denominator;
 703   TH1D *numerator;
 704   if(eta){
 705     int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
 706     int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
 707     //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 708     denominator = allhad->ProjectionX("name",lowbin,highbin);
 709     numerator = numeratorParent->ProjectionX("numerator",lowbin,highbin);
 710   }
 711   else{
 712     int lowbin = allhad->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
 713     int highbin = allhad->GetXaxis()->GetNbins();
 714     //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 715     numerator = numeratorParent->ProjectionY("name",lowbin,highbin);
 716     denominator = allhad->ProjectionY("denominator",lowbin,highbin);
 717   }
 718   numerator->Divide(denominator);
 719   //numerator->Rebin(2);
 720   //numerator->Scale(0.5);
 721   numerator->SetYTitle("E_{T}^{had,sample}/E_{T}^{had,total}");
 722   numerator->GetYaxis()->SetTitleOffset(1.2);
 723   numerator->SetMarkerColor(color);
 724   numerator->SetLineColor(color);
 725   numerator->SetMarkerStyle(marker);
 726   delete denominator;
 727   delete numeratorParent;
 728   delete allhad;
 729   //file->Close();
 730   numerator->SetName(name);
 731   return numerator;
 732
 733 }
 734
 735 //===============================CorrPtCut=========================================
 736 TH1D *GetHistoCorrPtCut(float ptcut, char *name, bool ispp, bool forSim, int mycase){
 737   file->cd();
 738   TH2F *allhad = ((TH2F*) out2->FindObject("EtSimulatedAllHadron"))->Clone("allhad");
 739   TH2F *ptlow = ((TH2F*) out2->FindObject("EtSimulatedChargedHadronAssumingNoPt"))->Clone("ptlow");
 740   TH2F *pthigh;
 741   if(ptcut>0.14){//TPC cut off
 742     (TH2F*)pthigh =(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedChargedHadronAssumingPtTPCCut"))->Clone("pthigh");
 743   }
 744   else{
 745     (TH2F*)pthigh =(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedChargedHadronAssumingPtITSCut"))->Clone("pthigh");
 746   }
 747
 748   int lowbin = allhad->GetXaxis()->FindBin(0.0);//make sure we don't accidentally get the wrong bin
 749   int highbin = allhad->GetXaxis()->FindBin(ptcut);
 750   int nbins = allhad->GetXaxis()->GetNbins();
 751   //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 752   //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(nbins)<<endl;
 753
 754   //allhad->Sumw2();
 755   TH1D *numerator;
 756   TH1D *denominator;
 757   switch(mycase){
 758   case -1:
 759     numerator = ptlow->ProjectionY("nameLow",lowbin,highbin);
 760     denominator = allhad->ProjectionY("denominatorLow",highbin,nbins);
 761     denominator->Add(ptlow);
 762     break;
 763   case 1:
 764     numerator = pthigh->ProjectionY("nameHigh",lowbin,highbin);
 765     denominator = allhad->ProjectionY("denominatorHigh",highbin,nbins);
 766     denominator->Add(pthigh);
 767     break;
 768   default:
 769     numerator = allhad->ProjectionY("name",lowbin,highbin);
 770     denominator = allhad->ProjectionY("denominator",lowbin,nbins);
 771   }
 772   numerator->Divide(denominator);
 773   numerator->SetYTitle("E_{T}^{had, p_{T}<cut-off}/E_{T}^{had, all p_{T}}");
 774   numerator->GetYaxis()->SetTitleOffset(1.);
 775   numerator->GetYaxis()->SetTitleSize(0.08);
 776   numerator->GetYaxis()->SetLabelSize(0.05);
 777   numerator->GetXaxis()->SetTitleSize(0.08);
 778   numerator->GetXaxis()->SetLabelSize(0.05);
 779   numerator->GetXaxis()->SetTitleOffset(.6);
 780   //numerator->Rebin(2);
 781   //numerator->Scale(0.5);
 782   //numerator->Draw("e");
 783   delete allhad;
 784   delete denominator;
 785   delete ptlow;
 786   delete pthigh;
 787   numerator->SetName(name);
 788   return numerator;
 789
 790 }
 791
 792 Float_t CorrPtCut(float ptcut, char *prodname, char *shortprodname, bool ispp, bool forSim, int mycase){
 793
 794   gStyle->SetOptTitle(0);
 795   gStyle->SetOptStat(0);
 796   gStyle->SetOptFit(0);
 797   TCanvas *c = new TCanvas("c","c",500,400);
 798   c->SetTopMargin(0.04);
 799   c->SetRightMargin(0.04);
 800   c->SetLeftMargin(0.181452);
 801   c->SetBottomMargin(0.134409);
 802   c->SetBorderSize(0);
 803   c->SetFillColor(0);
 804   c->SetFillColor(0);
 805   c->SetBorderMode(0);
 806   c->SetFrameFillColor(0);
 807   c->SetFrameBorderMode(0);
 808
 809
 810
 811   TH1D *High = GetHistoCorrPtCut(0.15-.001,"High",ispp,forSim,mycase);
 812   TH1D *Low = GetHistoCorrPtCut(0.1-.001,"Low",ispp,forSim,mycase);
 813   TH1D *Lowest = GetHistoCorrPtCut(0.05-.001,"Lowest",ispp,forSim,mycase);
 814
 815   TF1 *func = new TF1("func","[0]",-.7,.7);
 816   func->SetParameter(0,0.2);
 817   if(ptcut<.125){//its cuts
 818     Low->Fit(func);
 819   }
 820   else{//tpc cuts
 821     High->Fit(func);
 822   }
 823
 824   High->SetMaximum(0.04);
 825   High->SetMinimum(0.0);
 826   High->SetMarkerColor(2);
 827   Low->SetMarkerColor(4);
 828   High->SetLineColor(2);
 829   Low->SetLineColor(4);
 830   High->SetMinimum(0.0);
 831   High->SetMarkerStyle(20);
 832   Low->SetMarkerStyle(21);
 833   Lowest->SetMarkerStyle(22);
 834   High->Draw();
 835   Low->Draw("same");
 836   Lowest->Draw("same");
 837
 838
 839
 840   TLatex *tex = new TLatex(-0.723444,0.0373593,prodname);
 841   tex->SetTextSize(0.0537634);
 842   tex->Draw();
 843   TLegend *leg = new TLegend(0.217742,0.66129,0.477823,0.873656);
 844   leg->AddEntry(High,"p_{T} cut-off = 0.15 GeV/c");
 845   leg->AddEntry(Low,"p_{T} cut-off = 0.1 GeV/c");
 846   leg->AddEntry(Lowest,"p_{T} cut-off = 0.05 GeV/c");
 847   leg->SetFillStyle(0);
 848   leg->SetFillColor(0);
 849   leg->SetBorderSize(0);
 850   leg->SetTextSize(0.0537634);
 851   leg->Draw();
 852
 853   if(ptcut<.125){//its cuts
 854     c->SaveAs(Form("pics/%s/fptcutITS.eps",shortprodname));
 855     c->SaveAs(Form("pics/%s/fptcutITS.png",shortprodname));
 856   }
 857   else{//tpc cuts
 858     c->SaveAs(Form("pics/%s/fptcutTPC.eps",shortprodname));
 859     c->SaveAs(Form("pics/%s/fptcutTPC.png",shortprodname));
 860   }
 861
 862   float corr = func->GetParameter(0);
 863   //cout<<"Pt cut correction: "<<1.0/(1.0-corr)<<endl;
 864   delete High;
 865   delete Low;
 866   delete Lowest;
 867   delete func;
 868   delete c;
 869   delete tex;
 870   delete leg;
 871   return 1.0/(1.0-corr);
 872 }
 873
 874
 875
 876 //==================================CorrNotID=================================================
 877 TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, bool eta, bool ispp, bool forSim){
 878   file->cd();
 879   char *myname = mynameITS;
 880   if(TPC) myname = mynameTPC;
 881   TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentifiedAssumingPion",myname)))->Clone("notid");
 882   TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sUnidentified",myname)))->Clone("nNotid");
 883   if(!eta){
 884     //cout<<"Correction determined for all charged hadrons"<<endl;
 885     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)));
 886     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
 887     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
 888     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
 889     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
 890     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
 891     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sPiPlus",myname)));
 892     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sPiMinus",myname)));
 893     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sKPlus",myname)));
 894     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sKMinus",myname)));
 895     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sProton",myname)));
 896     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sAntiProton",myname)));
 897   }
 898
 899   TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentified",myname)))->Clone("id");
 900   if(!eta){
 901     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)));
 902     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
 903     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
 904     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
 905     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
 906     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
 907   }
 908
 909   TH1D *nNotidProj;
 910   TH1D *denominator;
 911   TH1D *numerator;
 912   if(eta){
 913     int lowbin = notid->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accv0entally get the wrong bin
 914     int highbin = notid->GetYaxis()->FindBin(etacut-.001);
 915     //cout<<"Projecting from "<<notid->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<notid->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 916     denominator = id->ProjectionX("name",lowbin,highbin);
 917     numerator = notid->ProjectionX("numerator",lowbin,highbin);
 918     nNotidProj = nNotid->ProjectionX("nNotidProj",lowbin,highbin);
 919   }
 920   else{
 921     //cout<<"Getting eta dependence"<<endl;
 922     int lowbin = id->GetXaxis()->FindBin(etacut);//make sure we don't accidentally get the wrong bin
 923     int highbin;
 924     if(etacut<0.15){//then we actually have ITS standalone tracks and we only want this to run from 0.1 to 0.15 because this will be used only for ITS standalone tracks
 925       highbin = id->GetXaxis()->FindBin(0.15);
 926     }
 927     else{
 928       highbin = id->GetXaxis()->GetNbins();
 929     }
 930     //cout<<"Projecting from "<<id->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<id->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 931     numerator = notid->ProjectionY("name",lowbin,highbin);
 932     denominator = id->ProjectionY("denominator",lowbin,highbin);
 933     nNotidProj = nNotid->ProjectionY("nNotidProj",lowbin,highbin);
 934   }
 935   TH1D *result = numerator;
 936   if(!denominator){
 937     cerr<<"Uh-oh!  Can't find denominator!!";
 938     return numerator;
 939   }
 940   else{result->Divide(denominator);}
 941   if(result->GetNbinsX() != nNotidProj->GetNbinsX()){
 942     cerr<<"Uh-oh!  Can't rescale errors! "<<result->GetNbinsX()<<"!="<<nNotidProj->GetNbinsX()<<endl;
 943     return result;
 944   }
 945   if(!nNotidProj){
 946     cerr<<"Uh-oh!  Can't find histogram!!";
 947     return numerator;
 948   }
 949   //fixing the errors
 950   for(int i=1;i<=result->GetNbinsX();i++){
 951     Float_t value = result->GetBinContent(i);
 952     Float_t valueerr = result->GetBinError(i);
 953     Float_t n = nNotidProj->GetBinContent(i);
 954     Float_t err;
 955     if(n<=0) err = 0.0;
 956     else{err= value/TMath::Power(n,0.5);}
 957     result->SetBinError(i,err);
 958     //cout<<"Was "<<valueerr<<", setting to "<<err<<endl;
 959   }
 960   result->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
 961   result->GetYaxis()->SetTitleOffset(1.2);
 962   delete denominator;
 963   delete nNotidProj;
 964   delete notid;
 965   delete nNotid;
 966   delete id;
 967   result->SetName(name);
 968   return result;
 969
 970 }
 971
 972 TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, bool ispp, bool forSim){
 973   gStyle->SetOptTitle(0);
 974   gStyle->SetOptStat(0);
 975   gStyle->SetOptFit(0);
 976   TCanvas *c = new TCanvas("c","c",500,400);
 977   c->SetTopMargin(0.04);
 978   c->SetRightMargin(0.04);
 979   c->SetBorderSize(0);
 980   c->SetFillColor(0);
 981   c->SetFillColor(0);
 982   c->SetBorderMode(0);
 983   c->SetFrameFillColor(0);
 984   c->SetFrameBorderMode(0);
 985
 986   TH1D *PHOS = GetHistoCorrNotID(etacut,name,TPC,true,ispp,forSim);
 987   PHOS->SetMarkerColor(2);
 988   PHOS->SetLineColor(2);
 989   PHOS->SetAxisRange(0.0,4);
 990   if(TPC){
 991     PHOS->SetMaximum(1.1);
 992     PHOS->SetMinimum(0.85);
 993   }
 994   else{
 995     PHOS->SetAxisRange(0.0,0.5);
 996   }
 997   PHOS->SetMarkerStyle(20);
 998   PHOS->Draw();
 999   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1000   tex->SetTextSize(0.0537634);
1001   tex->Draw();
1002   char *detector = detectorEMCAL;
1003   if(etacut<0.2) detector = detectorPHOS;
1004   if(TPC){
1005     sprintf(epsname,"pics/%s/fnotidTPC%s.eps",shortprodname,detector);
1006     sprintf(pngname,"pics/%s/fnotidTPC%s.png",shortprodname,detector);
1007   }
1008   else{
1009     sprintf(epsname,"pics/%s/fnotidITS%s.eps",shortprodname,detector);
1010     sprintf(pngname,"pics/%s/fnotidITS%s.png",shortprodname,detector);
1011   }
1012
1013   c->SaveAs(epsname);
1014   c->SaveAs(pngname);
1015   delete c;
1016   delete tex;
1017   PHOS->SetName(name);
1018   return PHOS;
1019 }
1020
1021 Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, bool ispp, bool forSim){
1022   if(!forSim){
1023     if(ispp){
1024       return 0.996;
1025     }
1026     else{
1027       return 0.976;
1028     }
1029   }
1030   gStyle->SetOptTitle(0);
1031   gStyle->SetOptStat(0);
1032   gStyle->SetOptFit(0);
1033   TCanvas *c = new TCanvas("c","c",500,400);
1034   c->SetTopMargin(0.04);
1035   c->SetRightMargin(0.04);
1036   c->SetBorderSize(0);
1037   c->SetFillColor(0);
1038   c->SetFillColor(0);
1039   c->SetBorderMode(0);
1040   c->SetFrameFillColor(0);
1041   c->SetFrameBorderMode(0);
1042
1043   TH1D *PHOS = GetHistoCorrNotID(ptcut,name,TPC,false,ispp,forSim);
1044   PHOS->SetMarkerColor(2);
1045   PHOS->SetLineColor(2);
1046   PHOS->SetMaximum(1.01);
1047   PHOS->SetMinimum(0.98);
1048   TF1 *func = new TF1("func","[0]",-etacut,etacut);
1049   PHOS->Fit(func,"","",-etacut,etacut);
1050   PHOS->SetMarkerStyle(20);
1051   PHOS->Draw();
1052   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1053   tex->SetTextSize(0.0537634);
1054   tex->Draw();
1055   char *detector = detectorEMCAL;
1056   if(etacut<0.2) detector = detectorPHOS;
1057   if(TPC){
1058     sprintf(epsname,"pics/%s/fnotidConstTPC%s.eps",shortprodname,detector);
1059     sprintf(pngname,"pics/%s/fnotidConstTPC%s.png",shortprodname,detector);
1060   }
1061   else{
1062     sprintf(epsname,"pics/%s/fnotidConstITS%s.eps",shortprodname,detector);
1063     sprintf(pngname,"pics/%s/fnotidConstITS%s.png",shortprodname,detector);
1064   }
1065
1066   c->SaveAs(epsname);
1067   c->SaveAs(pngname);
1068   delete c;
1069   delete PHOS;
1070   float value = func->GetParameter(0);
1071   delete func;
1072   delete tex;
1073   return value;
1074 }
1075
1076 //==================================CorrNoID=================================================
1077 TH1D *GetHistoNoID(float etacut, char *name, bool eta, bool TPC, bool ispp, bool forSim){
1078   file->cd();
1079   char *myname = mynameITS;
1080   if(TPC) myname = mynameTPC;
1081   TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadronAssumingPion",myname)))->Clone("notid");
1082   TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sChargedHadron",myname)))->Clone("nNotid");
1083
1084   TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadron",myname)))->Clone("id");
1085   int lowbin = id->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accidentally get the wrong bin
1086   int highbin = id->GetYaxis()->FindBin(etacut-.001);
1087
1088   TH1D *nNotidProj;
1089   TH1D *denominator;
1090   TH1D *numerator;
1091   if(eta){
1092     int lowbin = notid->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accv0entally get the wrong bin
1093     int highbin = notid->GetYaxis()->FindBin(etacut-.001);
1094     //cout<<"Projecting from "<<notid->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<notid->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
1095     denominator = id->ProjectionX("name",lowbin,highbin);
1096     numerator = notid->ProjectionX("numerator",lowbin,highbin);
1097     nNotidProj = nNotid->ProjectionX("nNotidProj",lowbin,highbin);
1098   }
1099   else{
1100     //cout<<"Getting eta dependence"<<endl;
1101     int lowbin = id->GetXaxis()->FindBin(etacut);//make sure we don't accidentally get the wrong bin
1102     int highbin;
1103     if(etacut<0.15){//then we actually have ITS standalone tracks and we only want this to run from 0.1 to 0.15 because this will be used only for ITS standalone tracks
1104       highbin = id->GetXaxis()->FindBin(0.15);
1105     }
1106     else{
1107       highbin = id->GetXaxis()->GetNbins();
1108     }
1109     //cout<<"Projecting from "<<id->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<id->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
1110     numerator = notid->ProjectionY("name",lowbin,highbin);
1111     denominator = id->ProjectionY("denominator",lowbin,highbin);
1112     nNotidProj = nNotid->ProjectionY("nNotidProj",lowbin,highbin);
1113   }
1114   if(denominator) numerator->Divide(denominator);
1115
1116   if(numerator->GetNbinsX() != nNotidProj->GetNbinsX()){
1117     cerr<<"Uh-oh!  Can't rescale errors! "<<numerator->GetNbinsX()<<"!="<<nNotidProj->GetNbinsX()<<endl;
1118     return numerator;
1119   }
1120   //fixing the errors
1121   for(int i=1;i<=numerator->GetNbinsX();i++){
1122     Float_t value = numerator->GetBinContent(i);
1123     Float_t valueerr = numerator->GetBinError(i);
1124     Float_t n = nNotidProj->GetBinContent(i);
1125     Float_t err=0.;
1126     if(n>0.0){
1127       err = value/TMath::Power(n,0.5);
1128     }
1129     numerator->SetBinError(i,err);;
1130   }
1131   numerator->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
1132   numerator->GetYaxis()->SetTitleOffset(1.2);
1133   delete denominator;
1134   delete nNotidProj;
1135   delete notid;
1136   delete nNotid;
1137   delete id;
1138   numerator->SetName(name);
1139   return numerator;
1140
1141 }
1142
1143 TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, bool ispp, bool forSim){
1144   gStyle->SetOptTitle(0);
1145   gStyle->SetOptStat(0);
1146   gStyle->SetOptFit(0);
1147   TCanvas *c = new TCanvas("c","c",500,400);
1148   c->SetTopMargin(0.04);
1149   c->SetRightMargin(0.04);
1150   c->SetBorderSize(0);
1151   c->SetFillColor(0);
1152   c->SetFillColor(0);
1153   c->SetBorderMode(0);
1154   c->SetFrameFillColor(0);
1155   c->SetFrameBorderMode(0);
1156
1157   TH1D *PHOS = GetHistoNoID(etacut,name,true,true,ispp,forSim);
1158   PHOS->SetMarkerColor(2);
1159   PHOS->SetLineColor(2);
1160   PHOS->SetAxisRange(0.0,4);
1161   PHOS->SetMaximum(1.1);
1162   PHOS->SetMinimum(0.85);
1163   PHOS->SetMarkerStyle(20);;
1164   PHOS->Draw();
1165   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1166   tex->SetTextSize(0.0537634);
1167   tex->Draw();
1168
1169
1170   char *detector = detectorEMCAL;
1171   if(etacut<0.2) detector = detectorPHOS;
1172   sprintf(epsname,"pics/%s/fnoid%s.eps",shortprodname,detector);
1173   sprintf(pngname,"pics/%s/fnoid%s.png",shortprodname,detector);
1174
1175   c->SaveAs(epsname);
1176   c->SaveAs(pngname);
1177   delete c;
1178   delete tex;
1179   PHOS->SetName(name);
1180   return PHOS;
1181
1182 }
1183
1184 Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, bool ispp, bool forSim){
1185   gStyle->SetOptTitle(0);
1186   gStyle->SetOptStat(0);
1187   gStyle->SetOptFit(0);
1188   TCanvas *c = new TCanvas("c","c",500,400);
1189   c->SetTopMargin(0.04);
1190   c->SetRightMargin(0.04);
1191   c->SetBorderSize(0);
1192   c->SetFillColor(0);
1193   c->SetFillColor(0);
1194   c->SetBorderMode(0);
1195   c->SetFrameFillColor(0);
1196   c->SetFrameBorderMode(0);
1197
1198   bool TPC = true;
1199   if(ptcut<.15) TPC = false;
1200   TH1D *PHOS = GetHistoNoID(ptcut,name,false,TPC,ispp,forSim);
1201   TF1 *func = new TF1("func","[0]",-etacut,etacut);
1202   PHOS->Fit(func,"","",-etacut,etacut);
1203   PHOS->SetMarkerColor(2);
1204   PHOS->SetLineColor(2);
1205   PHOS->SetMaximum(1.1);
1206   PHOS->SetMinimum(0.85);
1207   PHOS->SetMarkerStyle(20);;
1208   PHOS->Draw();
1209   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1210   tex->SetTextSize(0.0537634);
1211   tex->Draw();
1212
1213
1214   char *detector = detectorEMCAL;
1215   if(etacut<0.2) detector = detectorPHOS;
1216   if(TPC){
1217     sprintf(epsname,"pics/%s/fnoid%sTPC.eps",shortprodname,detector);
1218     sprintf(pngname,"pics/%s/fnoid%sTPC.png",shortprodname,detector);
1219   }
1220   else{
1221     sprintf(epsname,"pics/%s/fnoid%sITS.eps",shortprodname,detector);
1222     sprintf(pngname,"pics/%s/fnoid%sITS.png",shortprodname,detector);
1223   }
1224
1225   c->SaveAs(epsname);
1226   c->SaveAs(pngname);
1227   delete c;
1228   delete PHOS;
1229   float value = func->GetParameter(0);
1230   delete func;
1231   delete tex;
1232   return value;
1233
1234 }
1235 //==================================Efficiency=================================================
1236 TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name)
1237 {
1238     if(!numerator){
1239       cerr<<"Error:  numerator does not exist!"<<endl;
1240       return NULL;
1241     }
1242     if(!denominator){
1243       cerr<<"Error:  denominator does not exist!"<<endl;
1244       return NULL;
1245     }
1246     TH1D* result = (TH1D*)numerator->Clone(name);
1247     Int_t nbins = numerator->GetNbinsX();
1248     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
1249       Double_t numeratorVal = numerator->GetBinContent(ibin);
1250       Double_t denominatorVal = denominator->GetBinContent(ibin);
1251       // Check if the errors are right or the thing is scaled
1252       Double_t numeratorValErr = numerator->GetBinError(ibin);
1253       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
1254         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
1255         numeratorVal /= rescale;
1256       }
1257       Double_t denominatorValErr = denominator->GetBinError(ibin);
1258       if (!(denominatorValErr==0. || denominatorVal==0. )) {
1259         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
1260         denominatorVal /= rescale;
1261       }
1262       Double_t quotient = 0.;
1263       if (denominatorVal!=0.) {
1264         quotient = numeratorVal/denominatorVal;
1265       }
1266       Double_t quotientErr=0;
1267       if(numeratorVal>0.0 && denominatorVal>0.0 && (denominatorVal+2.0)>0.0 && (denominatorVal+3.0) >0.0){
1268         double argument = (numeratorVal+1.0)/(denominatorVal+2.0)*
1269           ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0));
1270         double test = TMath::Power(TMath::Abs(argument),0.5);
1271         quotientErr = TMath::Sqrt( TMath::Abs(
1272                                   (numeratorVal+1.0)/(denominatorVal+2.0)*
1273                                   ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))));
1274       }
1275       result->SetBinContent(ibin,quotient);
1276       result->SetBinError(ibin,quotientErr);
1277       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
1278     }
1279     return result;
1280 }
1281
1282
1283
1284 TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC,bool ITS, int cb, int cblast){
1285   bool eta = true;
1286   file->cd();
1287   char *myname = mynameITS;
1288   if(TPC&&!ITS) myname = mynameTPC;
1289   if(TPC&&ITS) myname = mynameTPCITS;
1290   TH2F *numeratorParent;
1291   switch(mycase){
1292   case 0:
1293     if(cblast != -1){//add more centrality bins
1294       for(int i=cb;i<=cblast;i++){
1295         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)))->Clone("RecoHadron");
1296         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)));}
1297         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiMinus",i)));
1298         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KMinus",i)));
1299         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)));
1300         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)));
1301         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"AntiProton",i)));
1302       }
1303     }
1304     else{
1305       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoHadron");
1306       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
1307       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
1308       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")));
1309       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")));
1310       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
1311     }
1312     break;
1313   case 1://pion
1314     if(cblast != -1){//add more centrality bins
1315       for(int i=cb;i<=cblast;i++){
1316         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)))->Clone("RecoPion");
1317         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)));}
1318         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiMinus",i)));
1319       }
1320     }
1321     else{
1322       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoPion");
1323       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
1324     }
1325     break;
1326   case 2://kaon
1327     if(cblast != -1){//add more centrality bins
1328       for(int i=cb;i<=cblast;i++){
1329         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)))->Clone("RecoKaon");
1330         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)));}
1331         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KMinus",i)));
1332       }
1333     }
1334     else{
1335       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")))->Clone("RecoKaon");
1336       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
1337     }
1338     break;
1339   case 3://proton
1340     if(cblast != -1){//add more centrality bins
1341       for(int i=cb;i<=cblast;i++){
1342         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)))->Clone("RecoProton");
1343         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)));}
1344         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"AntiProton",i)));
1345       }
1346     }
1347     else{
1348       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")))->Clone("RecoProton");
1349       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
1350     }
1351     break;
1352   case 4://electron
1353     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",myname,"EPlus",cbname)))->Clone("RecoElectron");
1354     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",myname,"EMinus",cbname)));
1355     break;
1356   }
1357   TH2F *denominatorParent;
1358   switch(mycase){
1359   case 0:
1360     if(cblast != -1){//add more centrality bins
1361       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedChargedHadronCB%i",cb)))->Clone("RecoHadron");
1362       for(int i=cb+1;i<=cblast;i++){
1363         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedChargedHadronCB%i",i)));
1364       }
1365     }
1366     else{
1367       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedChargedHadron"))->Clone("RecoHadron");
1368     }
1369     break;
1370   case 1://pion
1371     if(cblast != -1){//add more centrality bins
1372       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedPiPlusCB%i",cb)))->Clone("RecoPion");
1373       denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiMinusCB%i",cb)));
1374       for(int i=cb+1;i<=cblast;i++){
1375         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiPlusCB%i",i)));
1376         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiMinusCB%i",i)));
1377       }
1378     }
1379     else{
1380       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoPion");
1381       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
1382     }
1383     break;
1384   case 2://kaon
1385     if(cblast != -1){//add more centrality bins
1386       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedKPlusCB%i",cb)))->Clone("RecoKaon");
1387       denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKMinusCB%i",cb)));
1388       for(int i=cb+1;i<=cblast;i++){
1389         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKPlusCB%i",i)));
1390         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKMinusCB%i",i)));
1391       }
1392     }
1393     else{
1394       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedKPlus"))->Clone("RecoKaon");
1395       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
1396     }
1397     break;
1398   case 3://proton
1399     if(cblast != -1){//add more centrality bins
1400       for(int i=cb;i<=cblast;i++){
1401         if(cb==i)denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedProtonCB%i",i)))->Clone("RecoProton");
1402         else{denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedProtonCB%i",i)));}
1403         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedAntiProtonCB%i",i)));
1404       }
1405     }
1406     else{
1407       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedProton"))->Clone("RecoProton");
1408       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
1409     }
1410     break;
1411   case 4://electron
1412     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedEPlus"))->Clone("RecoElectron");
1413     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedEMinus"));
1414     break;
1415   }
1416   TH1D *denominator;
1417   TH1D *numerator;
1418   if(eta){
1419     int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
1420     int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
1421     //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
1422     denominator = denominatorParent->ProjectionX(Form("garbage%s",name),lowbin,highbin);
1423     numerator = numeratorParent->ProjectionX(name,lowbin,highbin);
1424   }
1425   else{
1426     int lowbin = denominatorParent->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
1427     int highbin = denominatorParent->GetXaxis()->GetNbins();
1428     //cout<<"Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
1429     numerator = numeratorParent->ProjectionY(name,lowbin,highbin);
1430     denominator = denominatorParent->ProjectionY(Form("denominator%s",name),lowbin,highbin);
1431   }
1432   delete numeratorParent;
1433   delete denominatorParent;
1434   //numerator->Divide(denominator);
1435   TH1D *result = bayneseffdiv((TH1D*) numerator,(TH1D*)denominator,name);
1436   //result->Rebin(2);
1437   //result->Scale(0.5);
1438   result->SetYTitle("Efficiency");
1439   result->GetYaxis()->SetTitleOffset(0.8);
1440   result->GetXaxis()->SetTitleOffset(0.8);
1441   result->GetYaxis()->SetLabelSize(0.05);
1442   result->GetXaxis()->SetLabelSize(0.05);
1443   result->GetYaxis()->SetTitleSize(0.05);
1444   result->GetXaxis()->SetTitleSize(0.05);
1445   result->SetMarkerColor(color);
1446   result->SetLineColor(color);
1447   result->SetMarkerStyle(marker);
1448   //result->SetName(name);
1449   //result->Draw("e");
1450   delete denominator;
1451   delete numerator;
1452   delete numeratorParent;
1453   delete denominatorParent;
1454   result->SetName(name);
1455   return result;
1456
1457 }
1458
1459 void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname){
1460   bool ITS = true;
1461   gStyle->SetOptTitle(0);
1462   gStyle->SetOptStat(0);
1463   gStyle->SetOptFit(0);
1464   TCanvas *c = new TCanvas("c","c",600,400);
1465   c->SetTopMargin(0.02);
1466   c->SetRightMargin(0.02);
1467   c->SetBorderSize(0);
1468   c->SetFillColor(0);
1469   c->SetFillColor(0);
1470   c->SetBorderMode(0);
1471   c->SetFrameFillColor(0);
1472   c->SetFrameBorderMode(0);
1473   //c->SetLogx();
1474
1475   int colortotal = 1;
1476   int colorpi = 2;
1477   int colork = 3;
1478   int colorp = 4;
1479   int phosmarker = 20;
1480   int emcalmarker = 24;
1481   float ptcut1 = 0.05;
1482   float ptcut2 = 0.1;
1483   TH1D *PHOStotal = GetHistoEfficiency(0.12,"PHOStotal",0,colortotal,phosmarker,TPC,ITS);
1484   TH1D *PHOSpi = GetHistoEfficiency(0.12,"PHOSpi",1,colorpi,phosmarker,TPC,ITS);
1485   TH1D *PHOSp = GetHistoEfficiency(0.12,"PHOSp",2,colork,phosmarker,TPC,ITS);
1486   TH1D *PHOSk = GetHistoEfficiency(0.12,"PHOSk",3,colorp,phosmarker,TPC,ITS);
1487   if(!TPC){PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));}
1488   else{PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.15),PHOStotal->GetXaxis()->FindBin(3.0));}
1489   PHOStotal->SetMinimum(0.0);
1490   PHOStotal->SetMaximum(1.0);
1491   PHOStotal->Draw();
1492   PHOSpi->Draw("same");
1493   PHOSp->Draw("same");
1494   PHOSk->Draw("same");
1495   TH1D *EMCALtotal = GetHistoEfficiency(0.7,"EMCALtotal",0,colortotal,emcalmarker,TPC,ITS);
1496   TH1D *EMCALpi = GetHistoEfficiency(0.7,"EMCALpi",1,colorpi,emcalmarker,TPC,ITS);
1497   TH1D *EMCALp = GetHistoEfficiency(0.7,"EMCALp",2,colork,emcalmarker,TPC,ITS);
1498   TH1D *EMCALk = GetHistoEfficiency(0.7,"EMCALk",3,colorp,emcalmarker,TPC,ITS);
1499   EMCALtotal->Draw("same");
1500   EMCALpi->Draw("same");
1501   EMCALp->Draw("same");
1502   EMCALk->Draw("same");
1503
1504
1505   TLegend *leg = new  TLegend(0.22651,0.247312,0.370805,0.438172);
1506   leg->AddEntry(PHOStotal,"#pi,K,p");
1507   leg->AddEntry(PHOSpi,"#pi^{#pm}");
1508   leg->AddEntry(PHOSk,"K^{#pm}");
1509   leg->AddEntry(PHOSp,"p,#bar{p}");
1510   leg->SetFillStyle(0);
1511   leg->SetFillColor(0);
1512   leg->SetBorderSize(0);
1513   leg->SetTextSize(0.06);
1514  leg->Draw();
1515
1516   TLine *line = new TLine(0.2,0.0,0.2,1.0);
1517   line->Draw();
1518   line->SetLineWidth(3.0);
1519   //line->SetLineColor(TColor::kYellow);
1520   line->SetLineStyle(2);
1521   TLatex *tex = new TLatex(0.497269,0.0513196,prodname);
1522   tex->SetTextSize(0.0537634);
1523   tex->Draw();
1524   TLatex *tex3 = new TLatex(1.16186,0.28348,"Closed symbols |#eta|<0.12 (PHOS)");
1525   tex3->SetTextSize(0.0537634);
1526   tex3->Draw();
1527   TLatex *tex4 = new TLatex(1.16186,0.213221,"Open symbols |#eta|<0.70 (EMCal)");
1528   tex4->SetTextSize(0.0537634);
1529   tex4->Draw();
1530   TLatex *tex2 = new TLatex(0.241937,0.448436,"Likely TPC cut-off 200 MeV/c");
1531   tex2->SetTextSize(0.0537634);
1532   tex2->Draw();
1533   if(TPC){
1534     if(ITS){
1535       sprintf(epsname,"pics/%s/CorrEfficiencyITSTPC.eps",shortprodname);
1536       sprintf(pngname,"pics/%s/CorrEfficiencyITSTPC.png",shortprodname);
1537     }
1538     else{
1539       sprintf(epsname,"pics/%s/CorrEfficiencyTPC.eps",shortprodname);
1540       sprintf(pngname,"pics/%s/CorrEfficiencyTPC.png",shortprodname);
1541     }
1542   }
1543   else{
1544     sprintf(epsname,"pics/%s/CorrEfficiencyITS.eps",shortprodname);
1545     sprintf(pngname,"pics/%s/CorrEfficiencyITS.png",shortprodname);
1546   }
1547   delete PHOStotal;
1548   delete PHOSpi;
1549   delete PHOSp;
1550   delete PHOSk;
1551   delete EMCALtotal;
1552   delete EMCALpi;
1553   delete EMCALp;
1554   delete EMCALk;
1555   delete leg;
1556   c->SaveAs(epsname);
1557   c->SaveAs(pngname);
1558   delete c;
1559   delete line;
1560   delete tex;
1561   delete tex3;
1562 }
1563
1564 //==================================CorrBkgd=================================================
1565 TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC,bool ispp,bool forSim){
1566   file->cd();
1567   char *reweightname = reweightedNo;
1568   if(!forSim) reweightname = reweightedYes;
1569   char *myname = mynameITS;
1570   if(TPC) myname = mynameTPC;
1571   TH2F *signal = ((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)))->Clone("signal");
1572   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
1573   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
1574   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
1575   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
1576   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
1577   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentifiedAssumingPion",myname)));
1578
1579   //Et of all unidentified hadrons (plus hadrons identified as pions) calculated assuming their true mass
1580   TH2F *bkgd = ((TH2F*) out2->FindObject(Form("EtReconstructed%sMisidentifiedElectrons",myname)))->Clone("bkgd");
1581   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sLambdaDaughters%s",myname,reweightname)));
1582   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiLambdaDaughters%s",myname,reweightname)));
1583   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sK0SDaughters%s",myname,reweightname)));
1584   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sXiDaughters",myname)));
1585   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiXiDaughters",myname)));
1586   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sOmegaDaughters",myname)));
1587   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiOmegaDaughters",myname)));
1588   int lowbin = bkgd->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accidentally get the wrong bin
1589   int highbin = bkgd->GetYaxis()->FindBin(etacut-.001);
1590   //cout<<"Projecting from "<<bkgd->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<bkgd->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
1591
1592
1593   TH1D *denominator = signal->ProjectionX(name,lowbin,highbin);
1594   TH1D *numerator = bkgd->ProjectionX("numerator",lowbin,highbin);
1595   numerator->Divide(denominator);
1596   numerator->SetYTitle("Ratio of E_{T}^{background}/E_{T}^{real}");
1597   numerator->GetYaxis()->SetTitleOffset(1.2);
1598   delete signal;
1599   delete bkgd;
1600   delete denominator;
1601   numerator->SetName(name);
1602   return numerator;
1603
1604 }
1605
1606 void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC,bool ispp,bool forSim){
1607   gStyle->SetOptTitle(0);
1608   gStyle->SetOptStat(0);
1609   gStyle->SetOptFit(0);
1610   TCanvas *c = new TCanvas("c","c",500,400);
1611   c->SetTopMargin(0.04);
1612   c->SetRightMargin(0.04);
1613   c->SetBorderSize(0);
1614   c->SetFillColor(0);
1615   c->SetFillColor(0);
1616   c->SetBorderMode(0);
1617   c->SetFrameFillColor(0);
1618   c->SetFrameBorderMode(0);
1619
1620   TH1D *PHOS = GetHistoCorrBkgd(0.12,"PHOS2",TPC,ispp,forSim);
1621   TH1D *EMCAL = GetHistoCorrBkgd(0.7,"EMCAL2",TPC,ispp,forSim);
1622   PHOS->SetMarkerColor(2);
1623   EMCAL->SetMarkerColor(4);
1624   PHOS->SetLineColor(2);
1625   EMCAL->SetLineColor(4);
1626   //EMCAL->SetLineWidth(2);
1627   //PHOS->SetAxisRange(0.0,4);
1628   PHOS->SetMaximum(0.2);
1629   PHOS->SetMinimum(0.0);
1630   PHOS->SetMarkerStyle(20);;
1631   EMCAL->SetMarkerStyle(21);
1632   //  TF1 *funcEMCAL = new TF1("funcEMCAL","[0]+0.0*x",0.05,4);
1633 //   funcEMCAL->SetParameter(0,0.95);
1634 //   funcEMCAL->SetParLimits(0,0.9,1.1);
1635   //EMCAL->Fit(funcEMCAL);//,"","",0.05,3.0);
1636 //   TF1 *funcPHOS = new TF1("funcPHOS","[0]+0.0*x",0.05,4);
1637 //   funcPHOS->SetParameter(0,1.0);
1638   //PHOS->Fit(funcPHOS);
1639   if(TPC)  PHOS->GetXaxis()->SetRange(PHOS->GetXaxis()->FindBin(0.0),PHOS->GetXaxis()->FindBin(4.));
1640   else{  PHOS->GetXaxis()->SetRange(PHOS->GetXaxis()->FindBin(0.0),PHOS->GetXaxis()->FindBin(1.));}
1641   PHOS->Draw();
1642   EMCAL->Draw("same");
1643   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1644   tex->SetTextSize(0.0537634);
1645   tex->Draw();
1646   TLegend *leg = new TLegend(0.145161,0.604839,0.40121,0.860215);
1647   leg->AddEntry(PHOS,"|#eta|<0.12");
1648   leg->AddEntry(EMCAL,"|#eta|<0.70");
1649   leg->SetFillStyle(0);
1650   leg->SetFillColor(0);
1651   leg->SetBorderSize(0);
1652   leg->Draw();
1653   if(TPC){
1654     c->SaveAs(Form("pics/%s/bkgdTPC.eps",shortprodname));
1655     c->SaveAs(Form("pics/%s/bkgdTPC.png",shortprodname));
1656   }
1657   else{
1658     c->SaveAs(Form("pics/%s/bkgdITS.eps",shortprodname));
1659     c->SaveAs(Form("pics/%s/bkgdITS.png",shortprodname));
1660   }
1661   delete c;
1662   delete PHOS;
1663   delete EMCAL;
1664   delete tex;
1665   delete leg;
1666
1667 }