]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/totEt/macros/GetCorrections.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added mock up of baryon enhancement, added code 20100 as lead collisions from 2010...
[u/mrichter/AliRoot.git] / PWG4 / totEt / macros / GetCorrections.C
CommitLineData
f427cbed 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
f43fc416 13Float_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);
14TH1D *GetHistoCorrNeutral(float cut, char *name, bool ispp, bool forSim, int mycase, bool eta, int color, int marker, char *infilename, bool hadronic = false);
f427cbed 15
f43fc416 16Float_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);
17TH1D *GetHistoCorrPtCut(float ptcut = 0.15, char *name, char *filename, bool ispp = true, bool forSim = true);
f427cbed 18
f43fc416 19TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta, bool ispp = true, bool forSim = true);
20TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp = true, bool forSim = true);
21Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim);
f427cbed 22
f43fc416 23TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC, bool ispp, bool forSim);
24TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim);
25Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim);
b64de20c 26
27TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name);
f43fc416 28TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC, bool ITS, char *infilename);
b64de20c 29void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname, char *infilename);
30
f43fc416 31TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC, char *infilename,bool ispp,bool forSim);
32void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC, char *infilename,bool ispp,bool forSim);
b64de20c 33
34//===========================================================================================
f427cbed 35
f43fc416 36void 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"){
f427cbed 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");
0e866ddc 52 gROOT->ProcessLine(".L AliAnalysisEtCuts.cxx+g");
f427cbed 53 gROOT->ProcessLine(".L AliAnalysisHadEtCorrections.cxx+g");
54
a02dfa56 55 char outfilename[200];
f43fc416 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);
f427cbed 61 TFile *outfile = new TFile(outfilename,"RECREATE");
62 AliAnalysisHadEtCorrections *hadCorrectionEMCAL = new AliAnalysisHadEtCorrections();
63 hadCorrectionEMCAL->SetName("hadCorrectionEMCAL");
b64de20c 64 float etacut = 0.7;
65 hadCorrectionEMCAL->SetEtaCut(etacut);
66 //float etacut = hadCorrectionEMCAL->GetEtaCut();
67 //cout<<"eta cut is "<<etacut<<endl;
f427cbed 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);
b64de20c 73
f427cbed 74 float ptcut = 0.1;
f43fc416 75 float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,false,etacut);
f427cbed 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
b64de20c 81
f43fc416 82 float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,true,etacut);
b64de20c 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);
f427cbed 87
f43fc416 88 float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,infilename,ispp,forSim);
f427cbed 89 hadCorrectionEMCAL->SetpTCutCorrectionITS(ptcutITS);
f43fc416 90 float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,infilename,ispp,forSim);
f427cbed 91 hadCorrectionEMCAL->SetpTCutCorrectionTPC(ptcutTPC);
0e866ddc 92 cout<<"Setting ITS pt cut corr to "<<ptcutITS<<endl;
93 cout<<"Setting TPC pt cut corr to "<<ptcutTPC<<endl;
f427cbed 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
f43fc416 100 TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDEMCALTPC",prodname,shortprodname,true,infilename,ispp,forSim);
101 TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDEMCALITS",prodname,shortprodname,false,infilename,ispp,forSim);
b64de20c 102 hadCorrectionEMCAL->SetNotIDCorrectionTPC(NotIDTPC);
103 hadCorrectionEMCAL->SetNotIDCorrectionITS(NotIDITS);
104
f43fc416 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);
0e866ddc 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
f43fc416 117 TH1D *NoID = CorrNoID(etacut,"CorrNoIDEMCAL",prodname,shortprodname,infilename,ispp,forSim);
b64de20c 118 hadCorrectionEMCAL->SetNotIDCorrectionNoPID(NoID);
119
f43fc416 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);
0e866ddc 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
f43fc416 131 TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,true,infilename);
a02dfa56 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
f43fc416 140 TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,true,infilename);
a02dfa56 141 if(!efficiencyKaonTPC){cerr<<"NOOOOOOOOOOOOOOOOOO!! We have failed you, Christine!"<<endl;}
142 hadCorrectionEMCAL->SetEfficiencyKaonTPC(efficiencyKaonTPC);
f43fc416 143 TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,true,infilename);
a02dfa56 144 hadCorrectionEMCAL->SetEfficiencyProtonTPC(efficiencyProtonTPC);
f43fc416 145 TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,true,infilename);
a02dfa56 146 hadCorrectionEMCAL->SetEfficiencyHadronTPC(efficiencyHadronTPC);
f43fc416 147 TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,true,infilename);
a02dfa56 148 hadCorrectionEMCAL->SetEfficiencyPionITS(efficiencyPionITS);
f43fc416 149 TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,true,infilename);
a02dfa56 150 hadCorrectionEMCAL->SetEfficiencyKaonITS(efficiencyKaonITS);
f43fc416 151 TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,true,infilename);
a02dfa56 152 hadCorrectionEMCAL->SetEfficiencyProtonITS(efficiencyProtonITS);
f43fc416 153 TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,true,infilename);
a02dfa56 154 hadCorrectionEMCAL->SetEfficiencyHadronITS(efficiencyHadronITS);
155
156 //CorrEfficiencyPlots(true,prodname,shortprodname,infilename);
157 //CorrEfficiencyPlots(false,prodname,shortprodname,infilename);
b64de20c 158
a02dfa56 159 hadCorrectionEMCAL->GetEfficiencyHadronTPC()->Draw();
f43fc416 160 TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,infilename,ispp,forSim);
161 TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,infilename,ispp,forSim);
b64de20c 162 hadCorrectionEMCAL->SetBackgroundCorrectionTPC(backgroundTPC);
163 hadCorrectionEMCAL->SetBackgroundCorrectionITS(backgroundITS);
f43fc416 164 CorrBkgdPlots(prodname,shortprodname,true,infilename,ispp,forSim);
165 CorrBkgdPlots(prodname,shortprodname,false,infilename,ispp,forSim);
f427cbed 166
f427cbed 167 outfile->cd();
168 hadCorrectionEMCAL->Write();
169 outfile->Write();
a02dfa56 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;
f43fc416 185 float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,false,etacut);
a02dfa56 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
f43fc416 192 float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,ispp,forSim,TPC,infilename,true,etacut);
a02dfa56 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
f43fc416 199 float ptcutITS = CorrPtCut(0.1,prodname,shortprodname,infilename,ispp,forSim);
a02dfa56 200 hadCorrectionPHOS->SetpTCutCorrectionITS(ptcutITS);
f43fc416 201 float ptcutTPC = CorrPtCut(0.15,prodname,shortprodname,infilename,ispp,forSim);
a02dfa56 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
f43fc416 209 TH1D *NotIDTPC = CorrNotID(etacut,"CorrNotIDPHOSTPC",prodname,shortprodname,true,infilename,ispp,forSim);
210 TH1D *NotIDITS = CorrNotID(etacut,"CorrNotIDPHOSITS",prodname,shortprodname,false,infilename,ispp,forSim);
a02dfa56 211 hadCorrectionPHOS->SetNotIDCorrectionTPC(NotIDTPC);
212 hadCorrectionPHOS->SetNotIDCorrectionITS(NotIDITS);
213
f43fc416 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);
0e866ddc 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
f43fc416 225 TH1D *NoID = CorrNoID(etacut,"CorrNoIDPHOS",prodname,shortprodname,infilename,ispp,forSim);
a02dfa56 226 hadCorrectionPHOS->SetNotIDCorrectionNoPID(NoID);
227
0e866ddc 228
f43fc416 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);
0e866ddc 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
f43fc416 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);
a02dfa56 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
f43fc416 259 TH1D *backgroundTPC = GetHistoCorrBkgd(etacut,"hBackgroundTPC",true,infilename,ispp,forSim);
260 TH1D *backgroundITS = GetHistoCorrBkgd(etacut,"hBackgroundITS",false,infilename,ispp,forSim);
a02dfa56 261 hadCorrectionPHOS->SetBackgroundCorrectionTPC(backgroundTPC);
262 hadCorrectionPHOS->SetBackgroundCorrectionITS(backgroundITS);
f43fc416 263 CorrBkgdPlots(prodname,shortprodname,true,infilename,ispp,forSim);
264 CorrBkgdPlots(prodname,shortprodname,false,infilename,ispp,forSim);
a02dfa56 265
266 //Write the output
267 outfile->cd();
268 hadCorrectionPHOS->Write();
269 outfile->Write();
f427cbed 270 outfile->Close();
271
f43fc416 272 TFile *junk = new TFile("junk.root","RECREATE");
273 efficiencyPionTPC->Write();
274 junk->Write();
275 junk->Close();
276
f427cbed 277 timer.Stop();
278 timer.Print();
279}
280
b64de20c 281//==================================CorrNeutral==============================================
f43fc416 282Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool ispp, bool forSim, bool TPC, char *infilename, bool hadronic, float etacut){
f427cbed 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;
b64de20c 315 int casetotal = 4;
316 if(hadronic) casetotal = 8;
f43fc416 317 TH1D *total = GetHistoCorrNeutral(ptcut,histoname,ispp,forSim,casetotal,false,colortotal,phosmarker,infilename,hadronic);
f427cbed 318
319 int colorallneutral = 2;
f43fc416 320 TH1D *allneutral = GetHistoCorrNeutral(ptcut,"allneutral",ispp,forSim,3,false,colorallneutral,phosmarker,infilename,hadronic);
f427cbed 321
322 int colorchargedsecondary = TColor::kViolet-3;
f43fc416 323 TH1D *chargedsecondary = GetHistoCorrNeutral(ptcut,"chargedsecondary",ispp,forSim,2,false,colorchargedsecondary,phosmarker,infilename,hadronic);
f427cbed 324
325 int colorneutralUndet = 4;
f43fc416 326 TH1D *neutralUndet = GetHistoCorrNeutral(ptcut,"neutralUndet",ispp,forSim,1,false,colorneutralUndet,phosmarker,infilename,hadronic);
f427cbed 327
328 int colorv0 = TColor::kGreen+2;
f43fc416 329 TH1D *v0 = GetHistoCorrNeutral(ptcut,"v0",ispp,forSim,0,false,colorv0,phosmarker,infilename,hadronic);
b64de20c 330
331 int colorem = TColor::kCyan;
f43fc416 332 TH1D *em = GetHistoCorrNeutral(ptcut,"em",ispp,forSim,9,false,colorem,phosmarker,infilename,hadronic);
f427cbed 333
334 TF1 *func = new TF1("func","[0]",-.7,.7);
335 func->SetParameter(0,0.2);
0e866ddc 336 total->Fit(func,"","",-etacut,etacut);
f427cbed 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);
b64de20c 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 }
f427cbed 351 total->Draw();
352 allneutral->Draw("same");
353 chargedsecondary->Draw("same");
354 neutralUndet->Draw("same");
355 v0->Draw("same");
b64de20c 356 if(hadronic) em->Draw("same");
f427cbed 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}");
b64de20c 365 if(hadronic) leg2->AddEntry(em,"e^{#pm},#gamma,#eta,#pi^{0},#omega");
f427cbed 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];
b64de20c 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 }
f427cbed 381 c->SaveAs(epsname);
382 c->SaveAs(pngname);
383
a02dfa56 384
385 delete total;
386 delete allneutral;
387 delete chargedsecondary;
388 delete neutralUndet;
389 delete v0;
390 delete em;
391 delete c;
392
f427cbed 393 float corr = func->GetParameter(0);
a02dfa56 394 //cout<<"Neutral correction: "<<1.0/(1.0-corr)<<endl;
395 delete func;
f427cbed 396 return 1.0/(1.0-corr);
397
398}
f43fc416 399TH1D *GetHistoCorrNeutral(float cut, char *name, bool ispp, bool forSim, int mycase, bool eta, int color, int marker, char *infilename, bool hadronic){
f427cbed 400 TFile *file = new TFile(infilename);
401 TList *list = file->FindObject("out2");
f43fc416 402 char *reweightname = "";
403 if(!forSim) reweightname = "Reweighted";
f427cbed 404 TH2F *numeratorParent;
405 switch(mycase){
406 case 0:
f43fc416 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)));
f427cbed 410 break;
411 case 1:
f43fc416 412 numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtSimulatedK0L%s",reweightname)))->Clone("Knnbar");
f427cbed 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:
f43fc416 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)));
f427cbed 429 numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedNeutron"));
430 numeratorParent->Add((TH2F*) out2->FindObject("EtSimulatedAntiNeutron"));
431 break;
432 case 4:
f43fc416 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)));
f427cbed 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;
b64de20c 460 case 8:
f43fc416 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)));
b64de20c 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;
f427cbed 488 }
489
490 TH2F *allhad;
f43fc416 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
b64de20c 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 }
f427cbed 523
f427cbed 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);
a02dfa56 529 //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
f427cbed 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();
a02dfa56 536 //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
f427cbed 537 numerator = numeratorParent->ProjectionY("name",lowbin,highbin);
538 denominator = allhad->ProjectionY("denominator",lowbin,highbin);
539 }
540 numerator->Divide(denominator);
0e866ddc 541 //numerator->Rebin(2);
542 //numerator->Scale(0.5);
f427cbed 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);
a02dfa56 548 delete denominator;
549 delete numeratorParent;
550 delete allhad;
551 //file->Close();
f427cbed 552 return numerator;
553
554}
555
b64de20c 556//===============================CorrPtCut=========================================
f43fc416 557TH1D *GetHistoCorrPtCut(float ptcut, char *name, char *filename, bool ispp, bool forSim){
f427cbed 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();
a02dfa56 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;
f427cbed 567
a02dfa56 568 //allhad->Sumw2();
f427cbed 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");
a02dfa56 583 delete allhad;
584 delete denominator;
585
f427cbed 586 return numerator;
587
588}
589
f43fc416 590Float_t CorrPtCut(float ptcut, char *prodname, char *shortprodname, char *filename, bool ispp, bool forSim){
f427cbed 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);
a02dfa56 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;
f427cbed 667 return 1.0/(1.0-corr);
668}
b64de20c 669
670
671
672//==================================CorrNotID=================================================
f43fc416 673TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta, bool ispp, bool forSim){
b64de20c 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");
0e866ddc 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 }
b64de20c 695
b64de20c 696 TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentified",myname)))->Clone("id");
0e866ddc 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 }
b64de20c 705
0e866ddc 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 }
b64de20c 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 }
b64de20c 742 if(!nNotidProj){
743 cerr<<"Uh-oh! Can't find histogram!!";
744 return numerator;
745 }
0e866ddc 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;
b64de20c 756 }
757 result->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
758 result->GetYaxis()->SetTitleOffset(1.2);
a02dfa56 759 delete denominator;
760 delete nNotidProj;
b64de20c 761 return result;
762
763}
764
f43fc416 765TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim){
b64de20c 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
f43fc416 779 TH1D *PHOS = GetHistoCorrNotID(etacut,name,TPC,infilename,true,ispp,forSim);
b64de20c 780 PHOS->SetMarkerColor(2);
781 PHOS->SetLineColor(2);
b64de20c 782 PHOS->SetAxisRange(0.0,4);
783 if(TPC){
784 PHOS->SetMaximum(1.1);
785 PHOS->SetMinimum(0.85);
786 }
787 else{
b64de20c 788 PHOS->SetAxisRange(0.0,0.5);
789 }
a02dfa56 790 PHOS->SetMarkerStyle(20);
b64de20c 791 PHOS->Draw();
792 TLatex *tex = new TLatex(0.161478,1.0835,prodname);
793 tex->SetTextSize(0.0537634);
794 tex->Draw();
b64de20c 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);
a02dfa56 810 delete c;
b64de20c 811 return PHOS;
812}
813
f43fc416 814Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename, bool ispp, bool forSim){
0e866ddc 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
f43fc416 828 TH1D *PHOS = GetHistoCorrNotID(ptcut,name,TPC,infilename,false,ispp,forSim);
0e866ddc 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
b64de20c 859//==================================CorrNoID=================================================
f43fc416 860TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC, bool ispp, bool forSim){
b64de20c 861 TFile *file = new TFile(infilename);
0e866ddc 862 char *myname = "ITS";
863 if(TPC) myname = "TPC";
b64de20c 864 TList *list = file->FindObject("out2");
0e866ddc 865 TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadronAssumingPion",myname)))->Clone("notid");
866 TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sChargedHadron",myname)))->Clone("nNotid");
b64de20c 867
0e866ddc 868 TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadron",myname)))->Clone("id");
b64de20c 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);
b64de20c 871
0e866ddc 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 }
a02dfa56 898 if(denominator) numerator->Divide(denominator);
b64de20c 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);
a02dfa56 909 Float_t err=0.;
910 if(n>0.0){
911 err = value/TMath::Power(n,0.5);
912 }
913 numerator->SetBinError(i,err);;
b64de20c 914 }
915 numerator->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
916 numerator->GetYaxis()->SetTitleOffset(1.2);
a02dfa56 917 delete denominator;
918 delete nNotidProj;
b64de20c 919 return numerator;
920
921}
922
f43fc416 923TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim){
b64de20c 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
f43fc416 937 TH1D *PHOS = GetHistoNoID(etacut,name,infilename,true,true,ispp,forSim);
b64de20c 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);
a02dfa56 959 delete c;
b64de20c 960 return PHOS;
961
0e866ddc 962}
963
f43fc416 964Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, char *infilename, bool ispp, bool forSim){
0e866ddc 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;
f43fc416 980 TH1D *PHOS = GetHistoNoID(ptcut,name,infilename,false,TPC,ispp,forSim);
0e866ddc 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
b64de20c 1012}
1013//==================================Efficiency=================================================
1014TH1D* 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
f43fc416 1062TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC,bool ITS, char *infilename){
b64de20c 1063 bool eta = true;
1064 TFile *file = new TFile(infilename);
1065 TList *list = file->FindObject("out2");
1066 char *myname = "ITS";
f43fc416 1067 if(TPC&&!ITS) myname = "TPC";
1068 if(TPC&&ITS) myname = "TPCITS";
1069 cout<<"Using tracks from "<<myname<<" for efficiency"<<endl;
b64de20c 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")));
a02dfa56 1079 //numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Unidentified")));
b64de20c 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 }
b64de20c 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);
a02dfa56 1125 //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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();
a02dfa56 1132 //cout<<"Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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);
a02dfa56 1152 //result->SetName(name);
b64de20c 1153 //result->Draw("e");
a02dfa56 1154 delete denominator;
1155 delete numerator;
b64de20c 1156 return result;
1157
1158}
1159
1160void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname, char *infilename){
f43fc416 1161 bool ITS = true;
b64de20c 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;
f43fc416 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);
b64de20c 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");
f43fc416 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);
b64de20c 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){
f43fc416 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 }
b64de20c 1245 }
1246 else{
1247 sprintf(epsname,"pics/%s/CorrEfficiencyITS.eps",shortprodname);
1248 sprintf(pngname,"pics/%s/CorrEfficiencyITS.png",shortprodname);
1249 }
a02dfa56 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;
b64de20c 1259 c->SaveAs(epsname);
1260 c->SaveAs(pngname);
a02dfa56 1261 delete c;
b64de20c 1262}
1263
1264//==================================CorrBkgd=================================================
f43fc416 1265TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC, char *infilename,bool ispp,bool forSim){
b64de20c 1266 TFile *file = new TFile(infilename);
1267 TList *list = file->FindObject("out2");
f43fc416 1268 char *reweightname = "";
1269 if(!forSim) reweightname = "Reweighted";
b64de20c 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");
f43fc416 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)));
b64de20c 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);
a02dfa56 1291 //cout<<"Projecting from "<<bkgd->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<bkgd->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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);
a02dfa56 1299 delete signal;
1300 delete bkgd;
1301 delete denominator;
b64de20c 1302 return numerator;
1303
1304}
1305
f43fc416 1306void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC, char *infilename,bool ispp,bool forSim){
b64de20c 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
f43fc416 1320 TH1D *PHOS = GetHistoCorrBkgd(0.12,"PHOS2",TPC,infilename,ispp,forSim);
1321 TH1D *EMCAL = GetHistoCorrBkgd(0.7,"EMCAL2",TPC,infilename,ispp,forSim);
b64de20c 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 }
a02dfa56 1361 delete c;
1362 delete PHOS;
1363 delete EMCAL;
b64de20c 1364
1365}