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