]> 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
Updated for changes
[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
0e866ddc 13Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool TPC, char *infilename, bool hadronic = false, float etacut = 0.7);
b64de20c 14TH1D *GetHistoCorrNeutral(float cut, char *name, int mycase, bool eta, int color, int marker, char *infilename, bool hadronic = false);
f427cbed 15
16Float_t CorrPtCut(float ptcut, char *prodname = "Enter Production Name", char *shortprodname = "EnterProductionName", char *filename="Et.ESD.new.sim.merged.root");
17TH1D *GetHistoCorrPtCut(float ptcut = 0.15, char *name, char *filename);
18
0e866ddc 19TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta);
b64de20c 20TH1D *CorrNotID(float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename);
0e866ddc 21Float_t CorrNotIDConst(float ptcut, float etacut,char *name, char *prodname, char *shortprodname, bool TPC, char *infilename);
f427cbed 22
0e866ddc 23TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC);
a02dfa56 24TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename);
0e866ddc 25Float_t CorrNoIDConst(float etacut, float ptcut,char *name, char *prodname, char *shortprodname, char *infilename);
b64de20c 26
27TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name);
28TH1D *GetHistoEfficiency(float cut, char *name, int mycase, int color, int marker,bool TPC, char *infilename);
29void CorrEfficiencyPlots(bool TPC, char *prodname, char *shortprodname, char *infilename);
30
31TH1D *GetHistoCorrBkgd(float etacut,char *name, bool TPC, char *infilename);
32void CorrBkgdPlots(char *prodname, char *shortprodname, bool TPC, char *infilename);
33
34//===========================================================================================
f427cbed 35
a02dfa56 36void GetCorrections(char *prodname = "Enter Production Name", char *shortprodname = "EnterProductionName", 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];
56 sprintf(outfilename,"corrections.%s.root",shortprodname);
f427cbed 57 TFile *outfile = new TFile(outfilename,"RECREATE");
58 AliAnalysisHadEtCorrections *hadCorrectionEMCAL = new AliAnalysisHadEtCorrections();
59 hadCorrectionEMCAL->SetName("hadCorrectionEMCAL");
b64de20c 60 float etacut = 0.7;
61 hadCorrectionEMCAL->SetEtaCut(etacut);
62 //float etacut = hadCorrectionEMCAL->GetEtaCut();
63 //cout<<"eta cut is "<<etacut<<endl;
f427cbed 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);
b64de20c 69
f427cbed 70 float ptcut = 0.1;
0e866ddc 71 float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,TPC,infilename,false,etacut);
f427cbed 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
b64de20c 77
0e866ddc 78 float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,TPC,infilename,true,etacut);
b64de20c 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);
f427cbed 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);
0e866ddc 88 cout<<"Setting ITS pt cut corr to "<<ptcutITS<<endl;
89 cout<<"Setting TPC pt cut corr to "<<ptcutTPC<<endl;
f427cbed 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
b64de20c 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
0e866ddc 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
a02dfa56 113 TH1D *NoID = CorrNoID(etacut,"CorrNoIDEMCAL",prodname,shortprodname,infilename);
b64de20c 114 hadCorrectionEMCAL->SetNotIDCorrectionNoPID(NoID);
115
0e866ddc 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
b64de20c 127 TH1D *efficiencyPionTPC = GetHistoEfficiency(etacut,"hEfficiencyPionTPC",1,1,20,true,infilename);
a02dfa56 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
b64de20c 136 TH1D *efficiencyKaonTPC = GetHistoEfficiency(etacut,"hEfficiencyKaonTPC",2,1,20,true,infilename);
a02dfa56 137 if(!efficiencyKaonTPC){cerr<<"NOOOOOOOOOOOOOOOOOO!! We have failed you, Christine!"<<endl;}
138 hadCorrectionEMCAL->SetEfficiencyKaonTPC(efficiencyKaonTPC);
b64de20c 139 TH1D *efficiencyProtonTPC = GetHistoEfficiency(etacut,"hEfficiencyProtonTPC",3,1,20,true,infilename);
a02dfa56 140 hadCorrectionEMCAL->SetEfficiencyProtonTPC(efficiencyProtonTPC);
b64de20c 141 TH1D *efficiencyHadronTPC = GetHistoEfficiency(etacut,"hEfficiencyHadronTPC",0,1,20,true,infilename);
a02dfa56 142 hadCorrectionEMCAL->SetEfficiencyHadronTPC(efficiencyHadronTPC);
b64de20c 143 TH1D *efficiencyPionITS = GetHistoEfficiency(etacut,"hEfficiencyPionITS",1,1,20,false,infilename);
a02dfa56 144 hadCorrectionEMCAL->SetEfficiencyPionITS(efficiencyPionITS);
b64de20c 145 TH1D *efficiencyKaonITS = GetHistoEfficiency(etacut,"hEfficiencyKaonITS",2,1,20,false,infilename);
a02dfa56 146 hadCorrectionEMCAL->SetEfficiencyKaonITS(efficiencyKaonITS);
b64de20c 147 TH1D *efficiencyProtonITS = GetHistoEfficiency(etacut,"hEfficiencyProtonITS",3,1,20,false,infilename);
a02dfa56 148 hadCorrectionEMCAL->SetEfficiencyProtonITS(efficiencyProtonITS);
b64de20c 149 TH1D *efficiencyHadronITS = GetHistoEfficiency(etacut,"hEfficiencyHadronITS",0,1,20,false,infilename);
a02dfa56 150 hadCorrectionEMCAL->SetEfficiencyHadronITS(efficiencyHadronITS);
151
152 //CorrEfficiencyPlots(true,prodname,shortprodname,infilename);
153 //CorrEfficiencyPlots(false,prodname,shortprodname,infilename);
b64de20c 154
a02dfa56 155 hadCorrectionEMCAL->GetEfficiencyHadronTPC()->Draw();
b64de20c 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);
f427cbed 162
f427cbed 163 outfile->cd();
164 hadCorrectionEMCAL->Write();
165 outfile->Write();
a02dfa56 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;
0e866ddc 181 float neutralCorr = CorrNeutral(ptcut,prodname,shortprodname,TPC,infilename,false,etacut);
a02dfa56 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
0e866ddc 188 float hadronicCorr = CorrNeutral(ptcut,prodname,shortprodname,TPC,infilename,true,etacut);
a02dfa56 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
0e866ddc 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
a02dfa56 221 TH1D *NoID = CorrNoID(etacut,"CorrNoIDPHOS",prodname,shortprodname,infilename);
222 hadCorrectionPHOS->SetNotIDCorrectionNoPID(NoID);
223
0e866ddc 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
a02dfa56 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();
f427cbed 266 outfile->Close();
267
268 timer.Stop();
269 timer.Print();
270}
271
b64de20c 272//==================================CorrNeutral==============================================
0e866ddc 273Float_t CorrNeutral(float ptcut, char *prodname, char *shortprodname, bool TPC, char *infilename, bool hadronic, float etacut){
f427cbed 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;
b64de20c 306 int casetotal = 4;
307 if(hadronic) casetotal = 8;
308 TH1D *total = GetHistoCorrNeutral(ptcut,histoname,casetotal,false,colortotal,phosmarker,infilename,hadronic);
f427cbed 309
310 int colorallneutral = 2;
b64de20c 311 TH1D *allneutral = GetHistoCorrNeutral(ptcut,"allneutral",3,false,colorallneutral,phosmarker,infilename,hadronic);
f427cbed 312
313 int colorchargedsecondary = TColor::kViolet-3;
b64de20c 314 TH1D *chargedsecondary = GetHistoCorrNeutral(ptcut,"chargedsecondary",2,false,colorchargedsecondary,phosmarker,infilename,hadronic);
f427cbed 315
316 int colorneutralUndet = 4;
b64de20c 317 TH1D *neutralUndet = GetHistoCorrNeutral(ptcut,"neutralUndet",1,false,colorneutralUndet,phosmarker,infilename,hadronic);
f427cbed 318
319 int colorv0 = TColor::kGreen+2;
b64de20c 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);
f427cbed 324
325 TF1 *func = new TF1("func","[0]",-.7,.7);
326 func->SetParameter(0,0.2);
0e866ddc 327 total->Fit(func,"","",-etacut,etacut);
f427cbed 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);
b64de20c 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 }
f427cbed 342 total->Draw();
343 allneutral->Draw("same");
344 chargedsecondary->Draw("same");
345 neutralUndet->Draw("same");
346 v0->Draw("same");
b64de20c 347 if(hadronic) em->Draw("same");
f427cbed 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}");
b64de20c 356 if(hadronic) leg2->AddEntry(em,"e^{#pm},#gamma,#eta,#pi^{0},#omega");
f427cbed 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];
b64de20c 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 }
f427cbed 372 c->SaveAs(epsname);
373 c->SaveAs(pngname);
374
a02dfa56 375
376 delete total;
377 delete allneutral;
378 delete chargedsecondary;
379 delete neutralUndet;
380 delete v0;
381 delete em;
382 delete c;
383
f427cbed 384 float corr = func->GetParameter(0);
a02dfa56 385 //cout<<"Neutral correction: "<<1.0/(1.0-corr)<<endl;
386 delete func;
f427cbed 387 return 1.0/(1.0-corr);
388
389}
b64de20c 390TH1D *GetHistoCorrNeutral(float cut, char *name, int mycase, bool eta, int color, int marker, char *infilename, bool hadronic){
f427cbed 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;
b64de20c 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;
f427cbed 477 }
478
479 TH2F *allhad;
480 allhad=(TH2F*) ((TH2F*) out2->FindObject("EtSimulatedAllHadron"))->Clone("id");
b64de20c 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 }
f427cbed 489
a02dfa56 490 //numeratorParent->Sumw2();
491 //allhad->Sumw2();
f427cbed 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);
a02dfa56 497 //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
f427cbed 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();
a02dfa56 504 //cout<<"Projecting from "<<allhad->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<allhad->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
f427cbed 505 numerator = numeratorParent->ProjectionY("name",lowbin,highbin);
506 denominator = allhad->ProjectionY("denominator",lowbin,highbin);
507 }
508 numerator->Divide(denominator);
0e866ddc 509 //numerator->Rebin(2);
510 //numerator->Scale(0.5);
f427cbed 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);
a02dfa56 516 delete denominator;
517 delete numeratorParent;
518 delete allhad;
519 //file->Close();
f427cbed 520 return numerator;
521
522}
523
b64de20c 524//===============================CorrPtCut=========================================
f427cbed 525TH1D *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();
a02dfa56 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;
f427cbed 535
a02dfa56 536 //allhad->Sumw2();
f427cbed 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");
a02dfa56 551 delete allhad;
552 delete denominator;
553
f427cbed 554 return numerator;
555
556}
557
558Float_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);
a02dfa56 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;
f427cbed 635 return 1.0/(1.0-corr);
636}
b64de20c 637
638
639
640//==================================CorrNotID=================================================
0e866ddc 641TH1D *GetHistoCorrNotID(float etacut,char *name, bool TPC, char *infilename, bool eta){
b64de20c 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");
0e866ddc 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 }
b64de20c 663
b64de20c 664 TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sUnidentified",myname)))->Clone("id");
0e866ddc 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 }
b64de20c 673
0e866ddc 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 }
b64de20c 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 }
b64de20c 710 if(!nNotidProj){
711 cerr<<"Uh-oh! Can't find histogram!!";
712 return numerator;
713 }
0e866ddc 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;
b64de20c 724 }
725 result->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
726 result->GetYaxis()->SetTitleOffset(1.2);
a02dfa56 727 delete denominator;
728 delete nNotidProj;
b64de20c 729 return result;
730
731}
732
733TH1D *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
0e866ddc 747 TH1D *PHOS = GetHistoCorrNotID(etacut,name,TPC,infilename,true);
b64de20c 748 PHOS->SetMarkerColor(2);
749 PHOS->SetLineColor(2);
b64de20c 750 PHOS->SetAxisRange(0.0,4);
751 if(TPC){
752 PHOS->SetMaximum(1.1);
753 PHOS->SetMinimum(0.85);
754 }
755 else{
b64de20c 756 PHOS->SetAxisRange(0.0,0.5);
757 }
a02dfa56 758 PHOS->SetMarkerStyle(20);
b64de20c 759 PHOS->Draw();
760 TLatex *tex = new TLatex(0.161478,1.0835,prodname);
761 tex->SetTextSize(0.0537634);
762 tex->Draw();
b64de20c 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);
a02dfa56 778 delete c;
b64de20c 779 return PHOS;
780}
781
0e866ddc 782Float_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
b64de20c 827//==================================CorrNoID=================================================
0e866ddc 828TH1D *GetHistoNoID(float etacut, char *name, char *infilename, bool eta, bool TPC){
b64de20c 829 TFile *file = new TFile(infilename);
0e866ddc 830 char *myname = "ITS";
831 if(TPC) myname = "TPC";
b64de20c 832 TList *list = file->FindObject("out2");
0e866ddc 833 TH2F *notid = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadronAssumingPion",myname)))->Clone("notid");
834 TH2F *nNotid = ((TH2F*) out2->FindObject(Form("EtNReconstructed%sChargedHadron",myname)))->Clone("nNotid");
b64de20c 835
0e866ddc 836 TH2F *id = ((TH2F*) out2->FindObject(Form("EtReconstructed%sChargedHadron",myname)))->Clone("id");
b64de20c 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);
b64de20c 839
0e866ddc 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 }
a02dfa56 866 if(denominator) numerator->Divide(denominator);
b64de20c 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);
a02dfa56 877 Float_t err=0.;
878 if(n>0.0){
879 err = value/TMath::Power(n,0.5);
880 }
881 numerator->SetBinError(i,err);;
b64de20c 882 }
883 numerator->SetYTitle("Ratio of E_{T}^{assuming pion}/E_{T}^{real}");
884 numerator->GetYaxis()->SetTitleOffset(1.2);
a02dfa56 885 delete denominator;
886 delete nNotidProj;
b64de20c 887 return numerator;
888
889}
890
a02dfa56 891TH1D *CorrNoID(float etacut,char *name, char *prodname, char *shortprodname, char *infilename){
b64de20c 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
0e866ddc 905 TH1D *PHOS = GetHistoNoID(etacut,name,infilename,true,true);
b64de20c 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);
a02dfa56 927 delete c;
b64de20c 928 return PHOS;
929
0e866ddc 930}
931
932Float_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
b64de20c 980}
981//==================================Efficiency=================================================
982TH1D* 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
1030TH1D *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")));
a02dfa56 1045 //numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Unidentified")));
b64de20c 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 }
b64de20c 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);
a02dfa56 1091 //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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();
a02dfa56 1098 //cout<<"Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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);
a02dfa56 1118 //result->SetName(name);
b64de20c 1119 //result->Draw("e");
a02dfa56 1120 delete denominator;
1121 delete numerator;
b64de20c 1122 return result;
1123
1124}
1125
1126void 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 }
a02dfa56 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;
b64de20c 1219 c->SaveAs(epsname);
1220 c->SaveAs(pngname);
a02dfa56 1221 delete c;
b64de20c 1222}
1223
1224//==================================CorrBkgd=================================================
1225TH1D *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);
a02dfa56 1249 //cout<<"Projecting from "<<bkgd->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<bkgd->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
b64de20c 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);
a02dfa56 1257 delete signal;
1258 delete bkgd;
1259 delete denominator;
b64de20c 1260 return numerator;
1261
1262}
1263
1264void 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
a02dfa56 1278 TH1D *PHOS = GetHistoCorrBkgd(0.12,"PHOS2",TPC,infilename);
1279 TH1D *EMCAL = GetHistoCorrBkgd(0.7,"EMCAL2",TPC,infilename);
b64de20c 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 }
a02dfa56 1319 delete c;
1320 delete PHOS;
1321 delete EMCAL;
b64de20c 1322
1323}