PWG4/totEt/macros/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, char *infilename, 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, char *infilename, bool hadronic = false);
  15
  16 Float_t CorrPtCut(float ptcut, char *prodname = "Enter Production Name", char *shortprodname = "EnterProductionName", char *filename="Et.ESD.new.sim.merged.root", bool ispp = true, bool forSim = true);
  17 TH1D *GetHistoCorrPtCut(float ptcut = 0.15, char *name, char *filename, bool ispp = true, bool forSim = true);
  18
  19 TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta, bool ispp = true, bool forSim = true);
  20 TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp = true, bool forSim = true);
  21 Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim);
  22
  23 TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC, bool ispp, bool forSim);
  24 TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim);
  25 Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, char *infilename, 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, char *infilename);
  29 void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname, char *infilename);
  30
  31 TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC, char *infilename,bool ispp,bool forSim);
  32 void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC, char *infilename,bool ispp,bool forSim);
  33
  34 //===========================================================================================
  35
  36 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"){
  37     TStopwatch timer;
  38     timer.Start();
  39     gSystem->Load("libTree.so");
  40     gSystem->Load("libGeom.so");
  41     gSystem->Load("libVMC.so");
  42     gSystem->Load("libXMLIO.so");
  43
  44     gSystem->Load("libSTEERBase.so");
  45     gSystem->Load("libESD.so");
  46     gSystem->Load("libAOD.so");
  47
  48     gSystem->Load("libANALYSIS");
  49     gSystem->Load("libANALYSISalice");
  50
  51     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  52    gROOT->ProcessLine(".L AliAnalysisEtCuts.cxx+g");
  53    gROOT->ProcessLine(".L AliAnalysisHadEtCorrections.cxx+g");
  54
  55    char outfilename[200];
  56    char *sim = "ForData";
  57    if(forSim) sim = "ForSimulations";
  58    char *system = "PbPb";
  59    if(ispp) system = "pp";
  60    sprintf(outfilename,"corrections.%s.%s.%s.root",shortprodname,system,sim);
  61    TFile *outfile = new TFile(outfilename,"RECREATE");
  62    AliAnalysisHadEtCorrections *hadCorrectionEMCAL = new AliAnalysisHadEtCorrections();
  63    hadCorrectionEMCAL->SetName("hadCorrectionEMCAL");
  64    float etacut = 0.7;
  65    hadCorrectionEMCAL->SetEtaCut(etacut);
  66    //float etacut = hadCorrectionEMCAL->GetEtaCut();
  67    //cout<<"eta cut is "<<etacut<<endl;
  68    cout<<"My name is "<<hadCorrectionEMCAL->GetName()<<endl;
  69    hadCorrectionEMCAL->SetAcceptanceCorrectionFull(1.0);
  70    cout<<"Warning:  Acceptance corrections will have to be updated to include real acceptance maps of the EMCAL and the PHOS"<<endl;
  71    hadCorrectionEMCAL->SetAcceptanceCorrectionPHOS(360.0/60.0);
  72    hadCorrectionEMCAL->SetAcceptanceCorrectionEMCAL(360.0/60.0);
  73
  74    float ptcut = 0.1;
  75    float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,false,etacut);
  76    hadCorrectionEMCAL->SetNeutralCorrection(neutralCorr);
  77    cout<<"Warning:  Setting neutral correction error bars to STAR value of +/-2%.  Use for development purposes only!"<<endl;
  78    hadCorrectionEMCAL->SetNeutralCorrectionLowBound(neutralCorr*0.98);
  79    hadCorrectionEMCAL->SetNeutralCorrectionHighBound(neutralCorr*1.02);
  80
  81
  82    float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,true,etacut);
  83    hadCorrectionEMCAL->SetNotHadronicCorrection(hadronicCorr);
  84    cout<<"Warning:  Setting hadronic correction error bars to value of +/-2%.  Use for development purposes only!"<<endl;
  85    hadCorrectionEMCAL->SetNotHadronicCorrectionLowBound(neutralCorr*0.98);
  86    hadCorrectionEMCAL->SetNotHadronicCorrectionHighBound(neutralCorr*1.02);
  87
  88    float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,infilename,ispp,forSim);
  89    hadCorrectionEMCAL->SetpTCutCorrectionITS(ptcutITS);
  90    float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,infilename,ispp,forSim);
  91    hadCorrectionEMCAL->SetpTCutCorrectionTPC(ptcutTPC);
  92    cout<<"Setting ITS pt cut corr to "<<ptcutITS<<endl;
  93    cout<<"Setting TPC pt cut corr to "<<ptcutTPC<<endl;
  94    cout<<"Warning:  Setting pt cut correction error bars to STAR value of +/-3%.  Use for development purposes only!"<<endl;
  95    hadCorrectionEMCAL->SetpTCutCorrectionITSLowBound(ptcutITS*0.97);
  96    hadCorrectionEMCAL->SetpTCutCorrectionTPCLowBound(ptcutTPC*0.97);
  97    hadCorrectionEMCAL->SetpTCutCorrectionITSHighBound(ptcutITS*1.03);
  98    hadCorrectionEMCAL->SetpTCutCorrectionTPCHighBound(ptcutTPC*1.03);
  99
 100    TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDEMCALTPC",prodname,shortprodname,true,infilename,ispp,forSim);
 101    TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDEMCALITS",prodname,shortprodname,false,infilename,ispp,forSim);
 102    hadCorrectionEMCAL->SetNotIDCorrectionTPC(NotIDTPC);
 103    hadCorrectionEMCAL->SetNotIDCorrectionITS(NotIDITS);
 104
 105    Float_t NotIDConstTPC = CorrNotIDConst(0.15,etacut,"CorrNotIDEMCALTPC2",prodname,shortprodname,true,infilename,ispp,forSim);
 106    Float_t NotIDConstITS = CorrNotIDConst(0.10,etacut,"CorrNotIDEMCALTPC2",prodname,shortprodname,true,infilename,ispp,forSim);
 107    hadCorrectionEMCAL->SetNotIDConstCorrectionTPC(1.0/NotIDConstTPC);
 108    hadCorrectionEMCAL->SetNotIDConstCorrectionITS(1.0/NotIDConstITS);
 109    cout<<"Setting constant PID corrections to "<<NotIDConstTPC<<" and "<<NotIDConstITS<<endl;
 110    cout<<"Warning:  Setting systematic errors on constant correction from unidentified particles at 1%!  For testing and development purposes only!"<<endl;
 111    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*.99);
 112    hadCorrectionEMCAL->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*.99);
 113    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*1.01);
 114    hadCorrectionEMCAL->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*1.01);
 115
 116
 117    TH1D *NoID = CorrNoID(etacut,"CorrNoIDEMCAL",prodname,shortprodname,infilename,ispp,forSim);
 118    hadCorrectionEMCAL->SetNotIDCorrectionNoPID(NoID);
 119
 120    Float_t NoIDTPC = CorrNoIDConst(etacut,0.15,"CorrNoIDEMCAL2",prodname,shortprodname,infilename,ispp,forSim);
 121    Float_t NoIDITS = CorrNoIDConst(etacut,0.1,"CorrNoIDEMCAL2",prodname,shortprodname,infilename,ispp,forSim);
 122    cout<<"Setting constant PID corrections with no PID to "<<NoIDTPC<<" and "<<NoIDITS<<endl;
 123    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoID(1./NoIDTPC);
 124    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoID(1./NoIDITS);
 125    cout<<"Warning:  Setting systematic errors on constant correction from unidentified particles at 2%!  For testing and development purposes only!"<<endl;
 126    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoIDLowBound(1./NoIDTPC*.99);
 127    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoIDLowBound(1./NoIDITS*.99);
 128    hadCorrectionEMCAL->SetNotIDConstCorrectionTPCNoIDHighBound(1./NoIDTPC*1.01);
 129    hadCorrectionEMCAL->SetNotIDConstCorrectionITSNoIDHighBound(1./NoIDITS*1.01);
 130
 131    TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,true,infilename);
 132    hadCorrectionEMCAL->SetEfficiencyPionTPC(efficiencyPionTPC);
 133    if(!efficiencyPionTPC){cerr<<"NOOOOOOOOOOOOOOOOOO!!  We have failed you, Christine!"<<endl;}
 134 //    else{
 135 //      hadCorrectionEMCAL->GetEfficiencyPionTPC()->Draw();
 136 //      cout<< "My name "<<hadCorrectionEMCAL->GetEfficiencyPionTPC()->GetName() <<endl;
 137 //      return;
 138 //    }
 139
 140    TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,true,infilename);
 141    if(!efficiencyKaonTPC){cerr<<"NOOOOOOOOOOOOOOOOOO!!  We have failed you, Christine!"<<endl;}
 142    hadCorrectionEMCAL->SetEfficiencyKaonTPC(efficiencyKaonTPC);
 143    TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,true,infilename);
 144    hadCorrectionEMCAL->SetEfficiencyProtonTPC(efficiencyProtonTPC);
 145    TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,true,infilename);
 146    hadCorrectionEMCAL->SetEfficiencyHadronTPC(efficiencyHadronTPC);
 147    TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,true,infilename);
 148    hadCorrectionEMCAL->SetEfficiencyPionITS(efficiencyPionITS);
 149    TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,true,infilename);
 150    hadCorrectionEMCAL->SetEfficiencyKaonITS(efficiencyKaonITS);
 151    TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,true,infilename);
 152    hadCorrectionEMCAL->SetEfficiencyProtonITS(efficiencyProtonITS);
 153    TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,true,infilename);
 154    hadCorrectionEMCAL->SetEfficiencyHadronITS(efficiencyHadronITS);
 155
 156    //CorrEfficiencyPlots(true,prodname,shortprodname,infilename);
 157    //CorrEfficiencyPlots(false,prodname,shortprodname,infilename);
 158
 159    hadCorrectionEMCAL->GetEfficiencyHadronTPC()->Draw();
 160    TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,infilename,ispp,forSim);
 161    TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,infilename,ispp,forSim);
 162    hadCorrectionEMCAL->SetBackgroundCorrectionTPC(backgroundTPC);
 163    hadCorrectionEMCAL->SetBackgroundCorrectionITS(backgroundITS);
 164    CorrBkgdPlots(prodname,shortprodname,true,infilename,ispp,forSim);
 165    CorrBkgdPlots(prodname,shortprodname,false,infilename,ispp,forSim);
 166
 167    outfile->cd();
 168    hadCorrectionEMCAL->Write();
 169    outfile->Write();
 170    delete hadCorrectionEMCAL;
 171
 172    AliAnalysisHadEtCorrections *hadCorrectionPHOS = new AliAnalysisHadEtCorrections();
 173    hadCorrectionPHOS->SetName("hadCorrectionPHOS");
 174    float etacut = 0.12;
 175    hadCorrectionPHOS->SetEtaCut(etacut);
 176    //float etacut = hadCorrectionPHOS->GetEtaCut();
 177    //cout<<"eta cut is "<<etacut<<endl;
 178    cout<<"My name is "<<hadCorrectionPHOS->GetName()<<endl;
 179    hadCorrectionPHOS->SetAcceptanceCorrectionFull(1.0);
 180    cout<<"Warning:  Acceptance corrections will have to be updated to include real acceptance maps of the PHOS and the PHOS"<<endl;
 181    hadCorrectionPHOS->SetAcceptanceCorrectionPHOS(360.0/60.0);
 182    hadCorrectionPHOS->SetAcceptanceCorrectionEMCAL(360.0/60.0);
 183
 184    float ptcut = 0.1;
 185    float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,false,etacut);
 186    hadCorrectionPHOS->SetNeutralCorrection(neutralCorr);
 187    cout<<"Warning:  Setting neutral correction error bars to STAR value of +/-2%.  Use for development purposes only!"<<endl;
 188    hadCorrectionPHOS->SetNeutralCorrectionLowBound(neutralCorr*0.98);
 189    hadCorrectionPHOS->SetNeutralCorrectionHighBound(neutralCorr*1.02);
 190
 191
 192    float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,true,etacut);
 193    hadCorrectionPHOS->SetNotHadronicCorrection(hadronicCorr);
 194    cout<<"Warning:  Setting hadronic correction error bars to value of +/-2%.  Use for development purposes only!"<<endl;
 195    hadCorrectionPHOS->SetNotHadronicCorrectionLowBound(neutralCorr*0.98);
 196    hadCorrectionPHOS->SetNotHadronicCorrectionHighBound(neutralCorr*1.02);
 197
 198
 199    float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,infilename,ispp,forSim);
 200    hadCorrectionPHOS->SetpTCutCorrectionITS(ptcutITS);
 201    float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,infilename,ispp,forSim);
 202    hadCorrectionPHOS->SetpTCutCorrectionTPC(ptcutTPC);
 203    cout<<"Warning:  Setting pt cut correction error bars to STAR value of +/-3%.  Use for development purposes only!"<<endl;
 204    hadCorrectionPHOS->SetpTCutCorrectionITSLowBound(ptcutITS*0.97);
 205    hadCorrectionPHOS->SetpTCutCorrectionTPCLowBound(ptcutTPC*0.97);
 206    hadCorrectionPHOS->SetpTCutCorrectionITSHighBound(ptcutITS*1.03);
 207    hadCorrectionPHOS->SetpTCutCorrectionTPCHighBound(ptcutTPC*1.03);
 208
 209    TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDPHOSTPC",prodname,shortprodname,true,infilename,ispp,forSim);
 210    TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDPHOSITS",prodname,shortprodname,false,infilename,ispp,forSim);
 211    hadCorrectionPHOS->SetNotIDCorrectionTPC(NotIDTPC);
 212    hadCorrectionPHOS->SetNotIDCorrectionITS(NotIDITS);
 213
 214    Float_t NotIDConstTPC = CorrNotIDConst(0.15,etacut,"CorrNotIDPHOSTPC2",prodname,shortprodname,true,infilename,ispp,forSim);
 215    Float_t NotIDConstITS = CorrNotIDConst(0.10,etacut,"CorrNotIDPHOSTPC2",prodname,shortprodname,true,infilename,ispp,forSim);
 216    hadCorrectionPHOS->SetNotIDConstCorrectionTPC(1./NotIDConstTPC);
 217    hadCorrectionPHOS->SetNotIDConstCorrectionITS(1./NotIDConstITS);
 218    cout<<"Warning:  Setting systematic errors on constant correction from unidentified particles at 1%!  For testing and development purposes only!"<<endl;
 219    hadCorrectionPHOS->SetNotIDConstCorrectionTPCLowBound(1./NotIDConstTPC*.99);
 220    hadCorrectionPHOS->SetNotIDConstCorrectionITSLowBound(1./NotIDConstITS*.99);
 221    hadCorrectionPHOS->SetNotIDConstCorrectionTPCHighBound(1./NotIDConstTPC*1.01);
 222    hadCorrectionPHOS->SetNotIDConstCorrectionITSHighBound(1./NotIDConstITS*1.01);
 223
 224
 225    TH1D *NoID = CorrNoID(etacut,"CorrNoIDPHOS",prodname,shortprodname,infilename,ispp,forSim);
 226    hadCorrectionPHOS->SetNotIDCorrectionNoPID(NoID);
 227
 228
 229    Float_t NoIDTPC = CorrNoIDConst(etacut,0.15,"CorrNoIDPHOS2",prodname,shortprodname,infilename,ispp,forSim);
 230    Float_t NoIDITS = CorrNoIDConst(etacut,0.1,"CorrNoIDPHOS2",prodname,shortprodname,infilename,ispp,forSim);
 231    cout<<"Setting constant PID corrections with no PID to "<<NoIDTPC<<" and "<<NoIDITS<<endl;
 232    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoID(1./NoIDTPC);
 233    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoID(1./NoIDITS);
 234    cout<<"Warning:  Setting systematic errors on constant correction from unidentified particles at 2%!  For testing and development purposes only!"<<endl;
 235    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoIDLowBound(1./NoIDTPC*.99);
 236    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoIDLowBound(1./NoIDITS*.99);
 237    hadCorrectionPHOS->SetNotIDConstCorrectionTPCNoIDHighBound(1./NoIDTPC*1.01);
 238    hadCorrectionPHOS->SetNotIDConstCorrectionITSNoIDHighBound(1./NoIDITS*1.01);
 239
 240    TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,true,infilename);
 241    TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,true,infilename);
 242    TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,true,infilename);
 243    TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,true,infilename);
 244    TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,true,infilename);
 245    TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,true,infilename);
 246    TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,true,infilename);
 247    TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,true,infilename);
 248    //CorrEfficiencyPlots(true,prodname,shortprodname,infilename);
 249    //CorrEfficiencyPlots(false,prodname,shortprodname,infilename);
 250    hadCorrectionPHOS->SetEfficiencyPionTPC(efficiencyPionTPC);
 251    hadCorrectionPHOS->SetEfficiencyKaonTPC(efficiencyKaonTPC);
 252    hadCorrectionPHOS->SetEfficiencyProtonTPC(efficiencyProtonTPC);
 253    hadCorrectionPHOS->SetEfficiencyHadronTPC(efficiencyHadronTPC);
 254    hadCorrectionPHOS->SetEfficiencyPionITS(efficiencyPionITS);
 255    hadCorrectionPHOS->SetEfficiencyKaonITS(efficiencyKaonITS);
 256    hadCorrectionPHOS->SetEfficiencyProtonITS(efficiencyProtonITS);
 257    hadCorrectionPHOS->SetEfficiencyHadronITS(efficiencyHadronITS);
 258
 259    TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,infilename,ispp,forSim);
 260    TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,infilename,ispp,forSim);
 261    hadCorrectionPHOS->SetBackgroundCorrectionTPC(backgroundTPC);
 262    hadCorrectionPHOS->SetBackgroundCorrectionITS(backgroundITS);
 263    CorrBkgdPlots(prodname,shortprodname,true,infilename,ispp,forSim);
 264    CorrBkgdPlots(prodname,shortprodname,false,infilename,ispp,forSim);
 265
 266    //Write the output
 267    outfile->cd();
 268    hadCorrectionPHOS->Write();
 269    outfile->Write();
 270    outfile->Close();
 271
 272    TFile *junk = new TFile("junk.root","RECREATE");
 273    efficiencyPionTPC->Write();
 274    junk->Write();
 275    junk->Close();
 276
 277   timer.Stop();
 278   timer.Print();
 279 }
 280
 281 //==================================CorrNeutral==============================================
 282 Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool ispp, bool forSim, bool TPC, char *infilename, bool hadronic, float etacut){
 283   gStyle->SetOptTitle(0);
 284   gStyle->SetOptStat(0);
 285   gStyle->SetOptFit(0);
 286   TCanvas *c = new TCanvas("c","c",400,400);
 287   c->SetTopMargin(0.0);
 288   c->SetRightMargin(0.0);
 289   c->SetBorderSize(0);
 290   c->SetFillColor(0);
 291   c->SetFillColor(0);
 292   c->SetBorderMode(0);
 293   c->SetFrameFillColor(0);
 294   c->SetFrameBorderMode(0);
 295
 296   TPad *ptpad = c->cd(1);
 297   ptpad->SetTopMargin(0.04);
 298   ptpad->SetRightMargin(0.04);
 299   ptpad->SetLeftMargin(0.149288);
 300   ptpad->SetBorderSize(0);
 301   ptpad->SetFillColor(0);
 302   ptpad->SetFillColor(0);
 303   ptpad->SetBorderMode(0);
 304   ptpad->SetFrameFillColor(0);
 305   ptpad->SetFrameBorderMode(0);
 306
 307   int phosmarker = 20;
 308
 309   char prefix[100];
 310   sprintf(prefix,"%s%2.1f",shortprodname,ptcut);
 311
 312   char histoname[100];
 313   sprintf(histoname,"%stotal",histoname);
 314   int colortotal = 1;
 315   int casetotal = 4;
 316   if(hadronic) casetotal = 8;
 317   TH1D *total = GetHistoCorrNeutral(ptcut,histoname,ispp,forSim,casetotal,false,colortotal,phosmarker,infilename,hadronic);
 318
 319   int colorallneutral = 2;
 320   TH1D *allneutral = GetHistoCorrNeutral(ptcut,"allneutral",ispp,forSim,3,false,colorallneutral,phosmarker,infilename,hadronic);
 321
 322   int colorchargedsecondary = TColor::kViolet-3;
 323   TH1D *chargedsecondary = GetHistoCorrNeutral(ptcut,"chargedsecondary",ispp,forSim,2,false,colorchargedsecondary,phosmarker,infilename,hadronic);
 324
 325   int colorneutralUndet = 4;
 326   TH1D *neutralUndet = GetHistoCorrNeutral(ptcut,"neutralUndet",ispp,forSim,1,false,colorneutralUndet,phosmarker,infilename,hadronic);
 327
 328   int colorv0 = TColor::kGreen+2;
 329   TH1D *v0 = GetHistoCorrNeutral(ptcut,"v0",ispp,forSim,0,false,colorv0,phosmarker,infilename,hadronic);
 330
 331   int colorem = TColor::kCyan;
 332   TH1D *em = GetHistoCorrNeutral(ptcut,"em",ispp,forSim,9,false,colorem,phosmarker,infilename,hadronic);
 333
 334   TF1 *func = new TF1("func","[0]",-.7,.7);
 335   func->SetParameter(0,0.2);
 336   total->Fit(func,"","",-etacut,etacut);
 337
 338   //total->SetAxisRange(0.0,4);
 339   total->GetXaxis()->SetLabelSize(0.05);
 340   total->GetYaxis()->SetLabelSize(0.045);
 341   total->GetXaxis()->SetTitleSize(0.05);
 342   total->GetYaxis()->SetTitleSize(0.06);
 343   if(hadronic){
 344     total->SetMaximum(0.6);
 345     total->SetMinimum(0.0);
 346   }
 347   else{
 348     total->SetMaximum(0.3);
 349     total->SetMinimum(0.0);
 350   }
 351   total->Draw();
 352   allneutral->Draw("same");
 353   chargedsecondary->Draw("same");
 354   neutralUndet->Draw("same");
 355   v0->Draw("same");
 356   if(hadronic) em->Draw("same");
 357
 358   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 359   tex->SetTextSize(0.0537634);
 360   tex->Draw();
 361   TLegend *leg2 = new TLegend(0.518321,0.746873,0.774812,0.955343);
 362   leg2->AddEntry(total,"Total");
 363   leg2->AddEntry(allneutral,"#Lambda,#bar{#Lambda},K^{0}_{S},K^{0}_{L},n,#bar{n}");
 364   leg2->AddEntry(neutralUndet,"K^{0}_{L},n,#bar{n}");
 365   if(hadronic) leg2->AddEntry(em,"e^{#pm},#gamma,#eta,#pi^{0},#omega");
 366   leg2->SetFillStyle(0);
 367   leg2->SetFillColor(0);
 368   leg2->SetBorderSize(0);
 369   leg2->SetTextSize(0.0548607);
 370   leg2->Draw();
 371   char epsname[100];
 372   char pngname[100];
 373   if(hadronic){
 374     sprintf(epsname,"pics/%s/fhadronic.eps",shortprodname);
 375     sprintf(pngname,"pics/%s/fhadronic.png",shortprodname);
 376   }
 377   else{
 378     sprintf(epsname,"pics/%s/fneutral.eps",shortprodname);
 379     sprintf(pngname,"pics/%s/fneutral.png",shortprodname);
 380   }
 381   c->SaveAs(epsname);
 382   c->SaveAs(pngname);
 383
 384
 385   delete total;
 386   delete allneutral;
 387   delete chargedsecondary;
 388   delete neutralUndet;
 389   delete v0;
 390   delete em;
 391   delete c;
 392
 393   float corr = func->GetParameter(0);
 394   //cout<<"Neutral correction: "<<1.0/(1.0-corr)<<endl;
 395   delete func;
 396   return 1.0/(1.0-corr);
 397
 398 }
 399 TH1D *GetHistoCorrNeutral(float cut, char *name, bool ispp, bool forSim, int mycase, bool eta, int color, int marker, char *infilename, bool hadronic){
 400   TFile *file = new TFile(infilename);
 401   TList *list = file->FindObject("out2");
 402   char *reweightname = "";
 403   if(!forSim) reweightname = "Reweighted";
 404   TH2F *numeratorParent;
 405   switch(mycase){
 406   case 0:
 407     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("v0");
 408     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 409     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 410     break;
 411   case 1:
 412     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)))->Clone("Knnbar");
 413     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 414     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 415     break;
 416   case 2:
 417     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedOmega"))->Clone("ch2ndary");
 418     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 419     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 420     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 421     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 422     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 423     break;
 424   case 3:
 425     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 426     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 427     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 428     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 429     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 430     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 431     break;
 432   case 4:
 433     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 434     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 435     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 436     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 437     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 438     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 439     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 440     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 441     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 442     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 443     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 444     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 445     break;
 446   case 5:
 447     numeratorParent= (TH2F*)((TH2F*) out2->FindObject("EtSimulatedXi"))->Clone("allxi");
 448     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 449     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 450     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 451     break;
 452   case 6:
 453     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedOmega"))->Clone("allomega");
 454     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 455     break;
 456   case 7:
 457     numeratorParent = (TH2F*)((TH2F*) out2->FindObject("EtSimulatedSigma"))->Clone("allsigma");
 458     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiSigma"));
 459     break;
 460   case 8:
 461     numeratorParent= (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)))->Clone("allneutral");
 462     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 463     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 464     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 465     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 466     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 467     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 468     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 469     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 470     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 471     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 472     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 473     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedGamma"));
 474     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 475     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 476     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 477     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 478     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 479     break;
 480   case 9:
 481     numeratorParent= (TH2F*)((TH2F*) out2->FindObject("EtSimulatedGamma"))->Clone("allem");
 482     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 483     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 484     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 485     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 486     numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 487     break;
 488   }
 489
 490   TH2F *allhad;
 491   //allhad=(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedAllHadron"))->Clone("id");
 492   allhad=(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedPiPlus"))->Clone("id");
 493   allhad->Add((TH2F*) out2->FindObject("EtSimulatedPiMinus"));
 494   allhad->Add((TH2F*) out2->FindObject("EtSimulatedKMinus"));
 495   allhad->Add((TH2F*) out2->FindObject("EtSimulatedKPlus"));
 496   allhad->Add((TH2F*) out2->FindObject("EtSimulatedProton"));
 497   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiProton"));
 498   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedLambda%s",reweightname)));
 499   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedAntiLambda%s",reweightname)));
 500   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0S%s",reweightname)));
 501   allhad->Add((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)));
 502   allhad->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
 503   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
 504   allhad->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 505   allhad->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 506   allhad->Add((TH2F*) out2->FindObject("EtSimulatedOmega"));
 507   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiOmega"));
 508   allhad->Add((TH2F*) out2->FindObject("EtSimulatedXi"));
 509   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi"));
 510   allhad->Add((TH2F*) out2->FindObject("EtSimulatedSigma"));
 511   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiSigma"));
 512   allhad->Add((TH2F*) out2->FindObject("EtSimulatedXi0"));
 513   allhad->Add((TH2F*) out2->FindObject("EtSimulatedAntiXi0"));
 514
 515   if(hadronic){//if we are getting the correction for the hadronic only case...
 516     allhad->Add((TH2F*) out2->FindObject("EtSimulatedGamma"));
 517     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEta"));
 518     allhad->Add((TH2F*) out2->FindObject("EtSimulatedPi0"));
 519     allhad->Add((TH2F*) out2->FindObject("EtSimulatedOmega0"));
 520     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEPlus"));
 521     allhad->Add((TH2F*) out2->FindObject("EtSimulatedEMinus"));
 522   }
 523
 524   TH1D *denominator;
 525   TH1D *numerator;
 526   if(eta){
 527     int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
 528     int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
 529     //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 530     denominator = allhad->ProjectionX("name",lowbin,highbin);
 531     numerator = numeratorParent->ProjectionX("numerator",lowbin,highbin);
 532   }
 533   else{
 534     int lowbin = allhad->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
 535     int highbin = allhad->GetXaxis()->GetNbins();
 536     //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 537     numerator = numeratorParent->ProjectionY("name",lowbin,highbin);
 538     denominator = allhad->ProjectionY("denominator",lowbin,highbin);
 539   }
 540   numerator->Divide(denominator);
 541   //numerator->Rebin(2);
 542   //numerator->Scale(0.5);
 543   numerator->SetYTitle("E_{T}^{had,sample}/E_{T}^{had,total}");
 544   numerator->GetYaxis()->SetTitleOffset(1.2);
 545   numerator->SetMarkerColor(color);
 546   numerator->SetLineColor(color);
 547   numerator->SetMarkerStyle(marker);
 548   delete denominator;
 549   delete numeratorParent;
 550   delete allhad;
 551   //file->Close();
 552   return numerator;
 553
 554 }
 555
 556 //===============================CorrPtCut=========================================
 557 TH1D *GetHistoCorrPtCut(float ptcut, char *name, char *filename, bool ispp, bool forSim){
 558   TFile *file = new TFile(filename);
 559   TList *list = file->FindObject("out2");
 560   TH2F *allhad = ((TH2F*) out2->FindObject("EtSimulatedAllHadron"))->Clone("allhad");
 561
 562   int lowbin = allhad->GetXaxis()->FindBin(0.0);//make sure we don't accidentally get the wrong bin
 563   int highbin = allhad->GetXaxis()->FindBin(ptcut);
 564   int nbins = allhad->GetXaxis()->GetNbins();
 565   //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 566   //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(nbins)<<endl;
 567
 568   //allhad->Sumw2();
 569
 570   TH1D *numerator = allhad->ProjectionY("name",lowbin,highbin);
 571   TH1D *denominator = allhad->ProjectionY("denominator",lowbin,nbins);
 572   numerator->Divide(denominator);
 573   numerator->SetYTitle("E_{T}^{had, p_{T}<cut-off}/E_{T}^{had, all p_{T}}");
 574   numerator->GetYaxis()->SetTitleOffset(1.);
 575   numerator->GetYaxis()->SetTitleSize(0.08);
 576   numerator->GetYaxis()->SetLabelSize(0.05);
 577   numerator->GetXaxis()->SetTitleSize(0.08);
 578   numerator->GetXaxis()->SetLabelSize(0.05);
 579   numerator->GetXaxis()->SetTitleOffset(.6);
 580   //numerator->Rebin(2);
 581   //numerator->Scale(0.5);
 582   //numerator->Draw("e");
 583   delete allhad;
 584   delete denominator;
 585
 586   return numerator;
 587
 588 }
 589
 590 Float_t CorrPtCut(float ptcut, char *prodname, char *shortprodname, char *filename, bool ispp, bool forSim){
 591
 592   gStyle->SetOptTitle(0);
 593   gStyle->SetOptStat(0);
 594   gStyle->SetOptFit(0);
 595   TCanvas *c = new TCanvas("c","c",500,400);
 596   c->SetTopMargin(0.04);
 597   c->SetRightMargin(0.04);
 598   c->SetLeftMargin(0.181452);
 599   c->SetBottomMargin(0.134409);
 600   c->SetBorderSize(0);
 601   c->SetFillColor(0);
 602   c->SetFillColor(0);
 603   c->SetBorderMode(0);
 604   c->SetFrameFillColor(0);
 605   c->SetFrameBorderMode(0);
 606
 607
 608
 609   TH1D *High = GetHistoCorrPtCut(0.15-.001,"High",filename);
 610   TH1D *Low = GetHistoCorrPtCut(0.1-.001,"Low",filename);
 611   TH1D *Lowest = GetHistoCorrPtCut(0.05-.001,"Lowest",filename);
 612
 613   TF1 *func = new TF1("func","[0]",-.7,.7);
 614   func->SetParameter(0,0.2);
 615   if(ptcut<.125){//its cuts
 616     Low->Fit(func);
 617   }
 618   else{//tpc cuts
 619     High->Fit(func);
 620   }
 621
 622   High->SetMaximum(0.04);
 623   High->SetMinimum(0.0);
 624   High->SetMarkerColor(2);
 625   Low->SetMarkerColor(4);
 626   High->SetLineColor(2);
 627   Low->SetLineColor(4);
 628   High->SetMinimum(0.0);
 629   High->SetMarkerStyle(20);
 630   Low->SetMarkerStyle(21);
 631   Lowest->SetMarkerStyle(22);
 632   High->Draw();
 633   Low->Draw("same");
 634   Lowest->Draw("same");
 635
 636
 637
 638   TLatex *tex = new TLatex(-0.723444,0.0373593,prodname);
 639   tex->SetTextSize(0.0537634);
 640   tex->Draw();
 641   TLegend *leg = new TLegend(0.217742,0.66129,0.477823,0.873656);
 642   leg->AddEntry(High,"p_{T} cut-off = 0.15 GeV/c");
 643   leg->AddEntry(Low,"p_{T} cut-off = 0.1 GeV/c");
 644   leg->AddEntry(Lowest,"p_{T} cut-off = 0.05 GeV/c");
 645   leg->SetFillStyle(0);
 646   leg->SetFillColor(0);
 647   leg->SetBorderSize(0);
 648   leg->SetTextSize(0.0537634);
 649   leg->Draw();
 650
 651   if(ptcut<.125){//its cuts
 652     c->SaveAs(Form("pics/%s/fptcutITS.eps",shortprodname));
 653     c->SaveAs(Form("pics/%s/fptcutITS.png",shortprodname));
 654   }
 655   else{//tpc cuts
 656     c->SaveAs(Form("pics/%s/fptcutTPC.eps",shortprodname));
 657     c->SaveAs(Form("pics/%s/fptcutTPC.png",shortprodname));
 658   }
 659
 660   float corr = func->GetParameter(0);
 661   //cout<<"Pt cut correction: "<<1.0/(1.0-corr)<<endl;
 662   delete High;
 663   delete Low;
 664   delete Lowest;
 665   delete func;
 666   delete c;
 667   return 1.0/(1.0-corr);
 668 }
 669
 670
 671
 672 //==================================CorrNotID=================================================
 673 TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta, bool ispp, bool forSim){
 674   TFile *file = new TFile(infilename);
 675   TList *list = file->FindObject("out2");
 676   char *myname = "ITS";
 677   if(TPC) myname = "TPC";
 678   TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentifiedAssumingPion",myname)))->Clone("notid");
 679   TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sUnidentified",myname)))->Clone("nNotid");
 680   if(!eta){
 681     cout<<"Correction determined for all charged hadrons"<<endl;
 682     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)));
 683     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
 684     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
 685     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
 686     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
 687     notid->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
 688     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sPiPlus",myname)));
 689     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sPiMinus",myname)));
 690     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sKPlus",myname)));
 691     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sKMinus",myname)));
 692     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sProton",myname)));
 693     nNotid->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%sAntiProton",myname)));
 694   }
 695
 696   TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentified",myname)))->Clone("id");
 697   if(!eta){
 698     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)));
 699     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
 700     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
 701     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
 702     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
 703     id->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
 704   }
 705
 706   TH1D *nNotidProj;
 707   TH1D *denominator;
 708   TH1D *numerator;
 709   if(eta){
 710     int lowbin = notid->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accv0entally get the wrong bin
 711     int highbin = notid->GetYaxis()->FindBin(etacut-.001);
 712     cout<<"Projecting from "<<notid->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<notid->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 713     denominator = id->ProjectionX("name",lowbin,highbin);
 714     numerator = notid->ProjectionX("numerator",lowbin,highbin);
 715     nNotidProj = nNotid->ProjectionX("nNotidProj",lowbin,highbin);
 716   }
 717   else{
 718     cout<<"Getting eta dependence"<<endl;
 719     int lowbin = id->GetXaxis()->FindBin(etacut);//make sure we don't accidentally get the wrong bin
 720     int highbin;
 721     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
 722       highbin = id->GetXaxis()->FindBin(0.15);
 723     }
 724     else{
 725       highbin = id->GetXaxis()->GetNbins();
 726     }
 727     cout<<"Projecting from "<<id->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<id->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 728     numerator = notid->ProjectionY("name",lowbin,highbin);
 729     denominator = id->ProjectionY("denominator",lowbin,highbin);
 730     nNotidProj = nNotid->ProjectionY("nNotidProj",lowbin,highbin);
 731   }
 732   TH1D *result = numerator;
 733   if(!denominator){
 734     cerr<<"Uh-oh!  Can't find denominator!!";
 735     return numerator;
 736   }
 737   else{result->Divide(denominator);}
 738   if(result->GetNbinsX() != nNotidProj->GetNbinsX()){
 739     cerr<<"Uh-oh!  Can't rescale errors! "<<result->GetNbinsX()<<"!="<<nNotidProj->GetNbinsX()<<endl;
 740     return result;
 741   }
 742   if(!nNotidProj){
 743     cerr<<"Uh-oh!  Can't find histogram!!";
 744     return numerator;
 745   }
 746   //fixing the errors
 747   for(int i=1;i<=result->GetNbinsX();i++){
 748     Float_t value = result->GetBinContent(i);
 749     Float_t valueerr = result->GetBinError(i);
 750     Float_t n = nNotidProj->GetBinContent(i);
 751     Float_t err;
 752     if(n<=0) err = 0.0;
 753     else{err= value/TMath::Power(n,0.5);}
 754     result->SetBinError(i,err);
 755     //cout<<"Was "<<valueerr<<", setting to "<<err<<endl;
 756   }
 757   result->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
 758   result->GetYaxis()->SetTitleOffset(1.2);
 759   delete denominator;
 760   delete nNotidProj;
 761   return result;
 762
 763 }
 764
 765 TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim){
 766   gStyle->SetOptTitle(0);
 767   gStyle->SetOptStat(0);
 768   gStyle->SetOptFit(0);
 769   TCanvas *c = new TCanvas("c","c",500,400);
 770   c->SetTopMargin(0.04);
 771   c->SetRightMargin(0.04);
 772   c->SetBorderSize(0);
 773   c->SetFillColor(0);
 774   c->SetFillColor(0);
 775   c->SetBorderMode(0);
 776   c->SetFrameFillColor(0);
 777   c->SetFrameBorderMode(0);
 778
 779   TH1D *PHOS = GetHistoCorrNotID(etacut,name,TPC,infilename,true,ispp,forSim);
 780   PHOS->SetMarkerColor(2);
 781   PHOS->SetLineColor(2);
 782   PHOS->SetAxisRange(0.0,4);
 783   if(TPC){
 784     PHOS->SetMaximum(1.1);
 785     PHOS->SetMinimum(0.85);
 786   }
 787   else{
 788     PHOS->SetAxisRange(0.0,0.5);
 789   }
 790   PHOS->SetMarkerStyle(20);
 791   PHOS->Draw();
 792   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 793   tex->SetTextSize(0.0537634);
 794   tex->Draw();
 795   char epsname[100];
 796   char pngname[100];
 797   char *detector = "EMCAL";
 798   if(etacut<0.2) detector = "PHOS";
 799   if(TPC){
 800     sprintf(epsname,"pics/%s/fnotidTPC%s.eps",shortprodname,detector);
 801     sprintf(pngname,"pics/%s/fnotidTPC%s.png",shortprodname,detector);
 802   }
 803   else{
 804     sprintf(epsname,"pics/%s/fnotidITS%s.eps",shortprodname,detector);
 805     sprintf(pngname,"pics/%s/fnotidITS%s.png",shortprodname,detector);
 806   }
 807
 808   c->SaveAs(epsname);
 809   c->SaveAs(pngname);
 810   delete c;
 811   return PHOS;
 812 }
 813
 814 Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim){
 815   gStyle->SetOptTitle(0);
 816   gStyle->SetOptStat(0);
 817   gStyle->SetOptFit(0);
 818   TCanvas *c = new TCanvas("c","c",500,400);
 819   c->SetTopMargin(0.04);
 820   c->SetRightMargin(0.04);
 821   c->SetBorderSize(0);
 822   c->SetFillColor(0);
 823   c->SetFillColor(0);
 824   c->SetBorderMode(0);
 825   c->SetFrameFillColor(0);
 826   c->SetFrameBorderMode(0);
 827
 828   TH1D *PHOS = GetHistoCorrNotID(ptcut,name,TPC,infilename,false,ispp,forSim);
 829   PHOS->SetMarkerColor(2);
 830   PHOS->SetLineColor(2);
 831   PHOS->SetMaximum(1.01);
 832   PHOS->SetMinimum(0.98);
 833   TF1 *func = new TF1("func","[0]",-etacut,etacut);
 834   PHOS->Fit(func,"","",-etacut,etacut);
 835   PHOS->SetMarkerStyle(20);
 836   PHOS->Draw();
 837   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 838   tex->SetTextSize(0.0537634);
 839   tex->Draw();
 840   char epsname[100];
 841   char pngname[100];
 842   char *detector = "EMCAL";
 843   if(etacut<0.2) detector = "PHOS";
 844   if(TPC){
 845     sprintf(epsname,"pics/%s/fnotidConstTPC%s.eps",shortprodname,detector);
 846     sprintf(pngname,"pics/%s/fnotidConstTPC%s.png",shortprodname,detector);
 847   }
 848   else{
 849     sprintf(epsname,"pics/%s/fnotidConstITS%s.eps",shortprodname,detector);
 850     sprintf(pngname,"pics/%s/fnotidConstITS%s.png",shortprodname,detector);
 851   }
 852
 853   c->SaveAs(epsname);
 854   c->SaveAs(pngname);
 855   delete c;
 856   return func->GetParameter(0);
 857 }
 858
 859 //==================================CorrNoID=================================================
 860 TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC, bool ispp, bool forSim){
 861   TFile *file = new TFile(infilename);
 862   char *myname = "ITS";
 863   if(TPC) myname = "TPC";
 864   TList *list = file->FindObject("out2");
 865   TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadronAssumingPion",myname)))->Clone("notid");
 866   TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sChargedHadron",myname)))->Clone("nNotid");
 867
 868   TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadron",myname)))->Clone("id");
 869   int lowbin = id->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accidentally get the wrong bin
 870   int highbin = id->GetYaxis()->FindBin(etacut-.001);
 871
 872   TH1D *nNotidProj;
 873   TH1D *denominator;
 874   TH1D *numerator;
 875   if(eta){
 876     int lowbin = notid->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accv0entally get the wrong bin
 877     int highbin = notid->GetYaxis()->FindBin(etacut-.001);
 878     cout<<"Projecting from "<<notid->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<notid->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 879     denominator = id->ProjectionX("name",lowbin,highbin);
 880     numerator = notid->ProjectionX("numerator",lowbin,highbin);
 881     nNotidProj = nNotid->ProjectionX("nNotidProj",lowbin,highbin);
 882   }
 883   else{
 884     cout<<"Getting eta dependence"<<endl;
 885     int lowbin = id->GetXaxis()->FindBin(etacut);//make sure we don't accidentally get the wrong bin
 886     int highbin;
 887     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
 888       highbin = id->GetXaxis()->FindBin(0.15);
 889     }
 890     else{
 891       highbin = id->GetXaxis()->GetNbins();
 892     }
 893     cout<<"Projecting from "<<id->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<id->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 894     numerator = notid->ProjectionY("name",lowbin,highbin);
 895     denominator = id->ProjectionY("denominator",lowbin,highbin);
 896     nNotidProj = nNotid->ProjectionY("nNotidProj",lowbin,highbin);
 897   }
 898   if(denominator) numerator->Divide(denominator);
 899
 900   if(numerator->GetNbinsX() != nNotidProj->GetNbinsX()){
 901     cerr<<"Uh-oh!  Can't rescale errors! "<<numerator->GetNbinsX()<<"!="<<nNotidProj->GetNbinsX()<<endl;
 902     return numerator;
 903   }
 904   //fixing the errors
 905   for(int i=1;i<=numerator->GetNbinsX();i++){
 906     Float_t value = numerator->GetBinContent(i);
 907     Float_t valueerr = numerator->GetBinError(i);
 908     Float_t n = nNotidProj->GetBinContent(i);
 909     Float_t err=0.;
 910     if(n>0.0){
 911       err = value/TMath::Power(n,0.5);
 912     }
 913     numerator->SetBinError(i,err);;
 914   }
 915   numerator->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
 916   numerator->GetYaxis()->SetTitleOffset(1.2);
 917   delete denominator;
 918   delete nNotidProj;
 919   return numerator;
 920
 921 }
 922
 923 TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim){
 924   gStyle->SetOptTitle(0);
 925   gStyle->SetOptStat(0);
 926   gStyle->SetOptFit(0);
 927   TCanvas *c = new TCanvas("c","c",500,400);
 928   c->SetTopMargin(0.04);
 929   c->SetRightMargin(0.04);
 930   c->SetBorderSize(0);
 931   c->SetFillColor(0);
 932   c->SetFillColor(0);
 933   c->SetBorderMode(0);
 934   c->SetFrameFillColor(0);
 935   c->SetFrameBorderMode(0);
 936
 937   TH1D *PHOS = GetHistoNoID(etacut,name,infilename,true,true,ispp,forSim);
 938   PHOS->SetMarkerColor(2);
 939   PHOS->SetLineColor(2);
 940   PHOS->SetAxisRange(0.0,4);
 941   PHOS->SetMaximum(1.1);
 942   PHOS->SetMinimum(0.85);
 943   PHOS->SetMarkerStyle(20);;
 944   PHOS->Draw();
 945   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 946   tex->SetTextSize(0.0537634);
 947   tex->Draw();
 948
 949
 950   char epsname[100];
 951   char pngname[100];
 952   char *detector = "EMCAL";
 953   if(etacut<0.2) detector = "PHOS";
 954   sprintf(epsname,"pics/%s/fnoid%s.eps",shortprodname,detector);
 955   sprintf(pngname,"pics/%s/fnoid%s.png",shortprodname,detector);
 956
 957   c->SaveAs(epsname);
 958   c->SaveAs(pngname);
 959   delete c;
 960   return PHOS;
 961
 962 }
 963
 964 Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim){
 965   gStyle->SetOptTitle(0);
 966   gStyle->SetOptStat(0);
 967   gStyle->SetOptFit(0);
 968   TCanvas *c = new TCanvas("c","c",500,400);
 969   c->SetTopMargin(0.04);
 970   c->SetRightMargin(0.04);
 971   c->SetBorderSize(0);
 972   c->SetFillColor(0);
 973   c->SetFillColor(0);
 974   c->SetBorderMode(0);
 975   c->SetFrameFillColor(0);
 976   c->SetFrameBorderMode(0);
 977
 978   bool TPC = true;
 979   if(ptcut<.15) TPC = false;
 980   TH1D *PHOS = GetHistoNoID(ptcut,name,infilename,false,TPC,ispp,forSim);
 981   TF1 *func = new TF1("func","[0]",-etacut,etacut);
 982   PHOS->Fit(func,"","",-etacut,etacut);
 983   PHOS->SetMarkerColor(2);
 984   PHOS->SetLineColor(2);
 985   PHOS->SetMaximum(1.1);
 986   PHOS->SetMinimum(0.85);
 987   PHOS->SetMarkerStyle(20);;
 988   PHOS->Draw();
 989   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
 990   tex->SetTextSize(0.0537634);
 991   tex->Draw();
 992
 993
 994   char epsname[100];
 995   char pngname[100];
 996   char *detector = "EMCAL";
 997   if(etacut<0.2) detector = "PHOS";
 998   if(TPC){
 999     sprintf(epsname,"pics/%s/fnoid%sTPC.eps",shortprodname,detector);
1000     sprintf(pngname,"pics/%s/fnoid%sTPC.png",shortprodname,detector);
1001   }
1002   else{
1003     sprintf(epsname,"pics/%s/fnoid%sITS.eps",shortprodname,detector);
1004     sprintf(pngname,"pics/%s/fnoid%sITS.png",shortprodname,detector);
1005   }
1006
1007   c->SaveAs(epsname);
1008   c->SaveAs(pngname);
1009   delete c;
1010   return func->GetParameter(0);
1011
1012 }
1013 //==================================Efficiency=================================================
1014 TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name)
1015 {
1016     if(!numerator){
1017       cerr<<"Error:  numerator does not exist!"<<endl;
1018       return NULL;
1019     }
1020     if(!denominator){
1021       cerr<<"Error:  denominator does not exist!"<<endl;
1022       return NULL;
1023     }
1024     TH1D* result = (TH1D*)numerator->Clone(name);
1025     Int_t nbins = numerator->GetNbinsX();
1026     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
1027       Double_t numeratorVal = numerator->GetBinContent(ibin);
1028       Double_t denominatorVal = denominator->GetBinContent(ibin);
1029       // Check if the errors are right or the thing is scaled
1030       Double_t numeratorValErr = numerator->GetBinError(ibin);
1031       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
1032         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
1033         numeratorVal /= rescale;
1034       }
1035       Double_t denominatorValErr = denominator->GetBinError(ibin);
1036       if (!(denominatorValErr==0. || denominatorVal==0. )) {
1037         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
1038         denominatorVal /= rescale;
1039       }
1040       Double_t quotient = 0.;
1041       if (denominatorVal!=0.) {
1042         quotient = numeratorVal/denominatorVal;
1043       }
1044       Double_t quotientErr=0;
1045       if(numeratorVal>0.0 && denominatorVal>0.0 && (denominatorVal+2.0)>0.0 && (denominatorVal+3.0) >0.0){
1046         double argument = (numeratorVal+1.0)/(denominatorVal+2.0)*
1047           ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0));
1048         double test = TMath::Power(TMath::Abs(argument),0.5);
1049         quotientErr = TMath::Sqrt( TMath::Abs(
1050                                   (numeratorVal+1.0)/(denominatorVal+2.0)*
1051                                   ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))));
1052       }
1053       result->SetBinContent(ibin,quotient);
1054       result->SetBinError(ibin,quotientErr);
1055       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
1056     }
1057     return result;
1058 }
1059
1060
1061
1062 TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC,bool ITS, char *infilename){
1063   bool eta = true;
1064   TFile *file = new TFile(infilename);
1065   TList *list = file->FindObject("out2");
1066   char *myname = "ITS";
1067   if(TPC&&!ITS) myname = "TPC";
1068   if(TPC&&ITS) myname = "TPCITS";
1069   cout<<"Using tracks from "<<myname<<" for efficiency"<<endl;
1070   TH2F *numeratorParent;
1071   switch(mycase){
1072   case 0:
1073     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoHadron");
1074     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
1075     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
1076     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")));
1077     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")));
1078     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
1079     //numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Unidentified")));
1080     break;
1081   case 1://pion
1082     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoPion");
1083     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
1084     break;
1085   case 2://kaon
1086     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")))->Clone("RecoKaon");
1087     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
1088     break;
1089   case 3://proton
1090     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")))->Clone("RecoProton");
1091     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
1092     break;
1093   case 4://electron
1094     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EPlus")))->Clone("RecoElectron");
1095     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"EMinus")));
1096     break;
1097   }
1098   TH2F *denominatorParent;
1099   switch(mycase){
1100   case 0:
1101     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedChargedHadron"))->Clone("RecoHadron");
1102     break;
1103   case 1://pion
1104     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoPion");
1105     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
1106     break;
1107   case 2://kaon
1108     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedKPlus"))->Clone("RecoKaon");
1109     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
1110     break;
1111   case 3://proton
1112     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedProton"))->Clone("RecoProton");
1113     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
1114     break;
1115   case 4://electron
1116     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedEPlus"))->Clone("RecoElectron");
1117     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedEMinus"));
1118     break;
1119   }
1120   TH1D *denominator;
1121   TH1D *numerator;
1122   if(eta){
1123     int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
1124     int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
1125     //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
1126     denominator = denominatorParent->ProjectionX(Form("garbage%s",name),lowbin,highbin);
1127     numerator = numeratorParent->ProjectionX(name,lowbin,highbin);
1128   }
1129   else{
1130     int lowbin = denominatorParent->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
1131     int highbin = denominatorParent->GetXaxis()->GetNbins();
1132     //cout<<"Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
1133     numerator = numeratorParent->ProjectionY(name,lowbin,highbin);
1134     denominator = denominatorParent->ProjectionY(Form("denominator%s",name),lowbin,highbin);
1135   }
1136   delete numeratorParent;
1137   delete denominatorParent;
1138   //numerator->Divide(denominator);
1139   TH1D *result = bayneseffdiv((TH1D*) numerator,(TH1D*)denominator,name);
1140   //result->Rebin(2);
1141   //result->Scale(0.5);
1142   result->SetYTitle("Efficiency");
1143   result->GetYaxis()->SetTitleOffset(0.8);
1144   result->GetXaxis()->SetTitleOffset(0.8);
1145   result->GetYaxis()->SetLabelSize(0.05);
1146   result->GetXaxis()->SetLabelSize(0.05);
1147   result->GetYaxis()->SetTitleSize(0.05);
1148   result->GetXaxis()->SetTitleSize(0.05);
1149   result->SetMarkerColor(color);
1150   result->SetLineColor(color);
1151   result->SetMarkerStyle(marker);
1152   //result->SetName(name);
1153   //result->Draw("e");
1154   delete denominator;
1155   delete numerator;
1156   return result;
1157
1158 }
1159
1160 void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname, char *infilename){
1161   bool ITS = true;
1162   gStyle->SetOptTitle(0);
1163   gStyle->SetOptStat(0);
1164   gStyle->SetOptFit(0);
1165   TCanvas *c = new TCanvas("c","c",600,400);
1166   c->SetTopMargin(0.02);
1167   c->SetRightMargin(0.02);
1168   c->SetBorderSize(0);
1169   c->SetFillColor(0);
1170   c->SetFillColor(0);
1171   c->SetBorderMode(0);
1172   c->SetFrameFillColor(0);
1173   c->SetFrameBorderMode(0);
1174   //c->SetLogx();
1175
1176   int colortotal = 1;
1177   int colorpi = 2;
1178   int colork = 3;
1179   int colorp = 4;
1180   int phosmarker = 20;
1181   int emcalmarker = 24;
1182   float ptcut1 = 0.05;
1183   float ptcut2 = 0.1;
1184   TH1D *PHOStotal = GetHistoEfficiency(0.12,"PHOStotal",0,colortotal,phosmarker,TPC,ITS,infilename);
1185   TH1D *PHOSpi = GetHistoEfficiency(0.12,"PHOSpi",1,colorpi,phosmarker,TPC,ITS,infilename);
1186   TH1D *PHOSp = GetHistoEfficiency(0.12,"PHOSp",2,colork,phosmarker,TPC,ITS,infilename);
1187   TH1D *PHOSk = GetHistoEfficiency(0.12,"PHOSk",3,colorp,phosmarker,TPC,ITS,infilename);
1188   if(!TPC){PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));}
1189   else{PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.15),PHOStotal->GetXaxis()->FindBin(3.0));}
1190   PHOStotal->SetMinimum(0.0);
1191   PHOStotal->SetMaximum(1.0);
1192   PHOStotal->Draw();
1193   PHOSpi->Draw("same");
1194   PHOSp->Draw("same");
1195   PHOSk->Draw("same");
1196   TH1D *EMCALtotal = GetHistoEfficiency(0.7,"EMCALtotal",0,colortotal,emcalmarker,TPC,ITS,infilename);
1197   TH1D *EMCALpi = GetHistoEfficiency(0.7,"EMCALpi",1,colorpi,emcalmarker,TPC,ITS,infilename);
1198   TH1D *EMCALp = GetHistoEfficiency(0.7,"EMCALp",2,colork,emcalmarker,TPC,ITS,infilename);
1199   TH1D *EMCALk = GetHistoEfficiency(0.7,"EMCALk",3,colorp,emcalmarker,TPC,ITS,infilename);
1200   EMCALtotal->Draw("same");
1201   EMCALpi->Draw("same");
1202   EMCALp->Draw("same");
1203   EMCALk->Draw("same");
1204
1205
1206   TLegend *leg = new  TLegend(0.22651,0.247312,0.370805,0.438172);
1207   leg->AddEntry(PHOStotal,"#pi,K,p");
1208   leg->AddEntry(PHOSpi,"#pi^{#pm}");
1209   leg->AddEntry(PHOSk,"K^{#pm}");
1210   leg->AddEntry(PHOSp,"p,#bar{p}");
1211   leg->SetFillStyle(0);
1212   leg->SetFillColor(0);
1213   leg->SetBorderSize(0);
1214   leg->SetTextSize(0.06);
1215  leg->Draw();
1216
1217   TLine *line = new TLine(0.2,0.0,0.2,1.0);
1218   line->Draw();
1219   line->SetLineWidth(3.0);
1220   //line->SetLineColor(TColor::kYellow);
1221   line->SetLineStyle(2);
1222   TLatex *tex = new TLatex(0.497269,0.0513196,prodname);
1223   tex->SetTextSize(0.0537634);
1224   tex->Draw();
1225   TLatex *tex3 = new TLatex(1.16186,0.28348,"Closed symbols |#eta|<0.12 (PHOS)");
1226   tex3->SetTextSize(0.0537634);
1227   tex3->Draw();
1228   TLatex *tex4 = new TLatex(1.16186,0.213221,"Open symbols |#eta|<0.70 (EMCal)");
1229   tex4->SetTextSize(0.0537634);
1230   tex4->Draw();
1231   TLatex *tex2 = new TLatex(0.241937,0.448436,"Likely TPC cut-off 200 MeV/c");
1232   tex2->SetTextSize(0.0537634);
1233   tex2->Draw();
1234   char epsname[100];
1235   char pngname[100];
1236   if(TPC){
1237     if(ITS){
1238       sprintf(epsname,"pics/%s/CorrEfficiencyITSTPC.eps",shortprodname);
1239       sprintf(pngname,"pics/%s/CorrEfficiencyITSTPC.png",shortprodname);
1240     }
1241     else{
1242       sprintf(epsname,"pics/%s/CorrEfficiencyTPC.eps",shortprodname);
1243       sprintf(pngname,"pics/%s/CorrEfficiencyTPC.png",shortprodname);
1244     }
1245   }
1246   else{
1247     sprintf(epsname,"pics/%s/CorrEfficiencyITS.eps",shortprodname);
1248     sprintf(pngname,"pics/%s/CorrEfficiencyITS.png",shortprodname);
1249   }
1250   delete PHOStotal;
1251   delete PHOSpi;
1252   delete PHOSp;
1253   delete PHOSk;
1254   delete EMCALtotal;
1255   delete EMCALpi;
1256   delete EMCALp;
1257   delete EMCALk;
1258   delete leg;
1259   c->SaveAs(epsname);
1260   c->SaveAs(pngname);
1261   delete c;
1262 }
1263
1264 //==================================CorrBkgd=================================================
1265 TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC, char *infilename,bool ispp,bool forSim){
1266   TFile *file = new TFile(infilename);
1267   TList *list = file->FindObject("out2");
1268   char *reweightname = "";
1269   if(!forSim) reweightname = "Reweighted";
1270   char *myname = "ITS";
1271   if(TPC) myname = "TPC";
1272   TH2F *signal = ((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiPlus",myname)))->Clone("signal");
1273   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedPiMinus",myname)));
1274   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKMinus",myname)));
1275   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedKPlus",myname)));
1276   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedProton",myname)));
1277   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sIdentifiedAntiProton",myname)));
1278   signal->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentifiedAssumingPion",myname)));
1279
1280   //Et of all unidentified hadrons (plus hadrons identified as pions) calculated assuming their true mass
1281   TH2F *bkgd = ((TH2F*) out2->FindObject(Form("EtReconstructed%sMisidentifiedElectrons",myname)))->Clone("bkgd");
1282   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sLambdaDaughters%s",myname,reweightname)));
1283   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiLambdaDaughters%s",myname,reweightname)));
1284   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sK0SDaughters%s",myname,reweightname)));
1285   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sXiDaughters",myname)));
1286   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiXiDaughters",myname)));
1287   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sOmegaDaughters",myname)));
1288   bkgd->Add((TH2F*) out2->FindObject(Form("EtReconstructed%sAntiOmegaDaughters",myname)));
1289   int lowbin = bkgd->GetYaxis()->FindBin(-etacut+.001);//make sure we don't accidentally get the wrong bin
1290   int highbin = bkgd->GetYaxis()->FindBin(etacut-.001);
1291   //cout<<"Projecting from "<<bkgd->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<bkgd->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
1292
1293
1294   TH1D *denominator = signal->ProjectionX(name,lowbin,highbin);
1295   TH1D *numerator = bkgd->ProjectionX("numerator",lowbin,highbin);
1296   numerator->Divide(denominator);
1297   numerator->SetYTitle("Ratio of E_{T}^{background}/E_{T}^{real}");
1298   numerator->GetYaxis()->SetTitleOffset(1.2);
1299   delete signal;
1300   delete bkgd;
1301   delete denominator;
1302   return numerator;
1303
1304 }
1305
1306 void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC, char *infilename,bool ispp,bool forSim){
1307   gStyle->SetOptTitle(0);
1308   gStyle->SetOptStat(0);
1309   gStyle->SetOptFit(0);
1310   TCanvas *c = new TCanvas("c","c",500,400);
1311   c->SetTopMargin(0.04);
1312   c->SetRightMargin(0.04);
1313   c->SetBorderSize(0);
1314   c->SetFillColor(0);
1315   c->SetFillColor(0);
1316   c->SetBorderMode(0);
1317   c->SetFrameFillColor(0);
1318   c->SetFrameBorderMode(0);
1319
1320   TH1D *PHOS = GetHistoCorrBkgd(0.12,"PHOS2",TPC,infilename,ispp,forSim);
1321   TH1D *EMCAL = GetHistoCorrBkgd(0.7,"EMCAL2",TPC,infilename,ispp,forSim);
1322   PHOS->SetMarkerColor(2);
1323   EMCAL->SetMarkerColor(4);
1324   PHOS->SetLineColor(2);
1325   EMCAL->SetLineColor(4);
1326   //EMCAL->SetLineWidth(2);
1327   //PHOS->SetAxisRange(0.0,4);
1328   PHOS->SetMaximum(0.2);
1329   PHOS->SetMinimum(0.0);
1330   PHOS->SetMarkerStyle(20);;
1331   EMCAL->SetMarkerStyle(21);
1332   //  TF1 *funcEMCAL = new TF1("funcEMCAL","[0]+0.0*x",0.05,4);
1333 //   funcEMCAL->SetParameter(0,0.95);
1334 //   funcEMCAL->SetParLimits(0,0.9,1.1);
1335   //EMCAL->Fit(funcEMCAL);//,"","",0.05,3.0);
1336 //   TF1 *funcPHOS = new TF1("funcPHOS","[0]+0.0*x",0.05,4);
1337 //   funcPHOS->SetParameter(0,1.0);
1338   //PHOS->Fit(funcPHOS);
1339   if(TPC)  PHOS->GetXaxis()->SetRange(PHOS->GetXaxis()->FindBin(0.0),PHOS->GetXaxis()->FindBin(4.));
1340   else{  PHOS->GetXaxis()->SetRange(PHOS->GetXaxis()->FindBin(0.0),PHOS->GetXaxis()->FindBin(1.));}
1341   PHOS->Draw();
1342   EMCAL->Draw("same");
1343   TLatex *tex = new TLatex(0.161478,1.0835,prodname);
1344   tex->SetTextSize(0.0537634);
1345   tex->Draw();
1346   TLegend *leg = new TLegend(0.145161,0.604839,0.40121,0.860215);
1347   leg->AddEntry(PHOS,"|#eta|<0.12");
1348   leg->AddEntry(EMCAL,"|#eta|<0.70");
1349   leg->SetFillStyle(0);
1350   leg->SetFillColor(0);
1351   leg->SetBorderSize(0);
1352   leg->Draw();
1353   if(TPC){
1354     c->SaveAs(Form("pics/%s/bkgdTPC.eps",shortprodname));
1355     c->SaveAs(Form("pics/%s/bkgdTPC.png",shortprodname));
1356   }
1357   else{
1358     c->SaveAs(Form("pics/%s/bkgdITS.eps",shortprodname));
1359     c->SaveAs(Form("pics/%s/bkgdITS.png",shortprodname));
1360   }
1361   delete c;
1362   delete PHOS;
1363   delete EMCAL;
1364
1365 }