]>
Commit | Line | Data |
---|---|---|
9f3689ba | 1 | //Oystein Djuvsland |
2 | //macro for plotting correction factors for EM et for use with AliAnalysisEt and writes them to a file | |
943dc129 | 3 | #include "TTree.h" |
4 | #include "TFile.h" | |
5 | #include <TList.h> | |
6 | #include <Rtypes.h> | |
7 | #include "TCanvas.h" | |
8 | #include <TH2I.h> | |
9 | #include <TH2F.h> | |
10 | #include <TProfile.h> | |
11 | #include <TFitResultPtr.h> | |
12 | #include <TFitResult.h> | |
13 | #include "TF1.h" | |
14 | #include <iostream> | |
15 | #include <AliAnalysisEtTrackMatchCorrections.h> | |
16 | #include <AliAnalysisEtRecEffCorrection.h> | |
17 | ||
18 | TCanvas *c1 = 0; | |
19 | ||
47151f26 | 20 | TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ; |
21 | TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ; | |
9f3689ba | 22 | |
23 | //creates an empty set of correction factors for debugging purposes | |
943dc129 | 24 | TF1 * generateRecEffFunction(Double_t p0, Double_t p1); |
9f3689ba | 25 | int createDummy(Double_t p0 = 1.0, Double_t p1 = 0.0, char *det = "Phos") |
02c62614 | 26 | { |
27 | TFile *outfile = TFile::Open("calocorrections.root","RECREATE"); | |
28 | TF1 fitneutral("fitneutral","0", 0, 100); | |
29 | TF1 fitcharged("fitcharged","0", 0, 100); | |
30 | TF1 fitgamma("fitgamma","0", 0, 100); | |
31 | TF1 fitsecondary("fitsecondary","0", 0, 100); | |
9f3689ba | 32 | AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(Form("TmCorrections%s",det), fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0); |
02c62614 | 33 | TF1 *func = generateRecEffFunction(p0, p1); |
9f3689ba | 34 | AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(Form("ReCorrections%s",det), *func, 1000); |
02c62614 | 35 | |
36 | cor->Write(); | |
37 | recor->Write(); | |
38 | outfile->Close(); | |
39 | ||
943dc129 | 40 | |
02c62614 | 41 | } |
9f3689ba | 42 | //p0 = correction factors for efficiency from track matching |
43 | //p1 = correction factors for efficiency from track matching | |
44 | //determined from fit of efficiency from track matching, 0.366 from simulations, fit as a function of energy | |
45 | //p0 is the constant | |
46 | //p1 is linear term | |
9793e78d | 47 | int calculateCorrections(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", Double_t p0 = 0.366, Double_t p1 = 0.0, char *det = "Emcal", Bool_t isSim = kFALSE) |
943dc129 | 48 | { |
9f3689ba | 49 | TString detector = det; |
9793e78d | 50 | TString emcal = "Emcal"; |
943dc129 | 51 | c1 = new TCanvas; |
52 | ||
53 | TFile *f = TFile::Open(filename, "READ"); | |
54 | ||
55 | TList *l = (TList*)f->Get("out1"); | |
56 | ||
9f3689ba | 57 | TTree *primaryTree = (TTree*)l->FindObject(("fPrimaryTree"+detector+"MC").Data()); |
58 | TTree *recTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"Rec").Data()); | |
59 | TTree *mcTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"MC").Data()); | |
943dc129 | 60 | std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl; |
61 | ||
62 | Int_t clusterMult = 0; | |
63 | Int_t nChargedNotRemoved = 0; | |
64 | Int_t nNeutralNotRemoved = 0; | |
65 | Int_t nGammaRemoved = 0; | |
66 | Int_t nSecNotRemoved = 0; | |
67 | ||
68 | Double_t emEtRec = 0; | |
69 | Double_t emEtMc = 0; | |
70 | ||
71 | recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult); | |
72 | mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved); | |
73 | mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved); | |
74 | mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved); | |
75 | mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved); | |
76 | ||
9793e78d | 77 | Float_t maxMult = 99.5; |
78 | Int_t nbins = 100; | |
d45a33b8 | 79 | if(detector==emcal){ |
9793e78d | 80 | // //100 seems to be sufficient... |
d45a33b8 | 81 | maxMult = 249.5; |
82 | nbins = 250; | |
83 | } | |
9793e78d | 84 | TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult); |
85 | TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult); | |
86 | TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult); | |
87 | TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult); | |
88 | ||
943dc129 | 89 | Int_t nEvents = mcTree->GetEntriesFast(); |
90 | for(Int_t i = 0; i < nEvents; i++) | |
91 | { | |
92 | mcTree->GetEvent(i); | |
93 | recTree->GetEvent(i); | |
94 | hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved); | |
95 | hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved); | |
96 | hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved); | |
97 | hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved); | |
98 | } | |
99 | ||
100 | c1->Divide(2,2); | |
101 | c1->cd(1); | |
102 | TString title = "Charged particles not removed"; | |
103 | TString xtitle = "Cluster Multiplicity"; | |
104 | TString ytitle = "N_{ch}"; | |
105 | hChargedVsClusterMult->SetTitle(title); | |
106 | hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
107 | hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
108 | hChargedVsClusterMult->SetStats(0); | |
109 | hChargedVsClusterMult->Draw(); | |
110 | TProfile *chProf = hChargedVsClusterMult->ProfileX(); | |
111 | chProf->SetStats(0); | |
112 | chProf->Draw("SAME"); | |
113 | //TF1 fitcharged("fitcharged","([0] + [1]*x)*(0.48/([2] + [3]*[2]))", 0, 100); | |
9793e78d | 114 | TF1 fitcharged("fitcharged","pol2", 0, 100);//fit of number of charged tracks vs detector multiplicity |
9f3689ba | 115 | //if straight line track matching roughly not dependent on centrality |
943dc129 | 116 | // fitcharged.FixParameter(2, p0); |
117 | // fitcharged.FixParameter(3, p1); | |
118 | TFitResultPtr chRes = chProf->Fit(&fitcharged,"S"); | |
119 | TArrayD ch; | |
120 | if(!chRes) | |
121 | { | |
122 | ch = TArrayD(chRes->NPar(),chRes->GetParams()); | |
123 | } | |
124 | else | |
125 | { | |
126 | std::cout << "Could not extract charged contribution params" << std::endl; | |
127 | } | |
128 | c1->cd(2); | |
129 | title = "Neutral particles not removed"; | |
130 | ytitle = "N_{neutral}"; | |
131 | hNeutralVsClusterMult->SetTitle(title); | |
132 | hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
133 | hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
134 | hNeutralVsClusterMult->SetStats(0); | |
135 | hNeutralVsClusterMult->Draw(); | |
136 | TProfile *neuProf = hNeutralVsClusterMult->ProfileX(); | |
137 | neuProf->SetStats(0); | |
138 | neuProf->Draw("SAME"); | |
9793e78d | 139 | TF1 fitneutral("fitneutral","pol2", 0, 100);//fit of number of neutral particles in calo that we call background in calo vs detector multiplicity |
9f3689ba | 140 | //may include K0S |
943dc129 | 141 | TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S"); |
142 | TArrayD neu; | |
143 | if(!neuRes) | |
144 | { | |
145 | neu = TArrayD(neuRes->NPar(),neuRes->GetParams()); | |
146 | } | |
147 | else | |
148 | { | |
149 | std::cout << "Could not extract charged contribution params" << std::endl; | |
150 | } | |
151 | c1->cd(3); | |
152 | ||
153 | title = "Gammas removed"; | |
154 | ytitle = "N_{#gamma}"; | |
155 | hGammaVsClusterMult->SetTitle(title); | |
156 | hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
157 | hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
158 | hGammaVsClusterMult->SetStats(0); | |
159 | hGammaVsClusterMult->Draw(); | |
160 | TProfile *gamProf = hGammaVsClusterMult->ProfileX(); | |
161 | gamProf->SetStats(0); | |
162 | gamProf->Draw("SAME"); | |
9793e78d | 163 | TF1 fitgamma("fitgamma","pol2", 0, 100);//fit of number of gammas removed erroneously vs detector multiplicity |
943dc129 | 164 | TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S"); |
165 | TArrayD gamma; | |
166 | if(!gammaRes) | |
167 | { | |
168 | gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams()); | |
169 | } | |
170 | else | |
171 | { | |
172 | std::cout << "Could not extract charged contribution params" << std::endl; | |
173 | } | |
174 | ||
175 | c1->cd(4); | |
176 | ||
177 | title = "Secondaries not removed"; | |
178 | ytitle = "N_{sec}"; | |
179 | hSecVsClusterMult->SetTitle(title); | |
180 | hSecVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
181 | hSecVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
182 | hSecVsClusterMult->SetStats(0); | |
183 | hSecVsClusterMult->Draw(); | |
184 | TProfile *secProf = hSecVsClusterMult->ProfileX(); | |
185 | secProf->SetStats(0); | |
186 | secProf->Draw("SAME"); | |
9793e78d | 187 | TF1 fitsecondary("fitsecondary","pol2", 0, 100);//fit of number of secondary particles that leave deposits in calo vs detector multiplicity |
943dc129 | 188 | TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S"); |
189 | TArrayD sec; | |
190 | if(!secRes) | |
191 | { | |
192 | sec = TArrayD(secRes->NPar(),secRes->GetParams()); | |
193 | } | |
194 | else | |
195 | { | |
196 | std::cout << "Could not extract charged contribution params" << std::endl; | |
197 | } | |
9f3689ba | 198 | //ugly hack for getting the energy |
199 | //changing number of particles to energy of particles by multiplying by the <energy> and dividing by the efficiency for the average energy | |
200 | //average energy from each of these: in excel file | |
943dc129 | 201 | Double_t meanCharged = 0.48/(p0 + p1*0.48); |
202 | Double_t meanNeutral = 0.53/(p0 + p1*0.53); | |
203 | Double_t meanGamma = 0.51/(p0 + p1*0.51); | |
204 | Double_t meanSecondary = meanGamma; | |
205 | ||
95309578 | 206 | TH2F *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllCent"); |
207 | TH2F *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoCent"); | |
208 | TH2F *eff2D; | |
d45a33b8 | 209 | // if(fHistHadronDepositsRecoMult && fHistHadronDepositsAllMult ){ |
95309578 | 210 | eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D"); |
d45a33b8 | 211 | // } |
212 | // else{ | |
213 | // cerr<<"Warning! Did not calculate reconstruction efficiency!!"<<endl; | |
214 | // eff2D = new TH2F("eff2D", "eff2D",200, 0, 10,20,-0.5,19.5); | |
215 | // } | |
47151f26 | 216 | //cor->SetReconstructionEfficiency(eff2D); |
217 | ||
95309578 | 218 | TH2F *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent"); |
219 | TH2F *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent"); | |
220 | TH2F *gammaEff2D; | |
d45a33b8 | 221 | //if(fHistGammasGeneratedMult && fHistGammasFoundMult){ |
95309578 | 222 | gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D"); |
d45a33b8 | 223 | // } |
224 | // else{ | |
225 | // cerr<<"Warning! Did not calculate reconstruction efficiency!!"<<endl; | |
226 | // gammaEff2D = new TH2F("gammaEff2D", "gammaEff2D",200, 0, 10,20,-0.5,19.5); | |
227 | // } | |
47151f26 | 228 | |
229 | AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D, | |
943dc129 | 230 | meanCharged, meanNeutral, meanGamma, meanSecondary ); |
231 | ||
47151f26 | 232 | |
943dc129 | 233 | TF1 *func = generateRecEffFunction(p0, p1); |
47151f26 | 234 | AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000); |
943dc129 | 235 | TFile *outfile = TFile::Open("calocorrections.root","RECREATE"); |
236 | cor->Write(); | |
237 | recor->Write(); | |
238 | return 0; | |
239 | } | |
240 | ||
241 | ||
242 | TF1* generateRecEffFunction(Double_t p0, Double_t p1) | |
243 | { | |
244 | TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200); | |
245 | ||
246 | Double_t params[2] = {p0, p1}; | |
247 | f->SetParameters(params); | |
248 | return f; | |
249 | } | |
47151f26 | 250 | |
251 | TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) | |
252 | { | |
253 | if(!numerator){ | |
254 | cerr<<"Error: numerator does not exist!"<<endl; | |
255 | return NULL; | |
256 | } | |
257 | if(!denominator){ | |
258 | cerr<<"Error: denominator does not exist!"<<endl; | |
259 | return NULL; | |
260 | } | |
261 | TH1* result = (TH1*)numerator->Clone(name); | |
262 | Int_t nbins = numerator->GetNbinsX(); | |
263 | for (Int_t ibin=0; ibin<= nbins+1; ++ibin) { | |
264 | Double_t numeratorVal = numerator->GetBinContent(ibin); | |
265 | Double_t denominatorVal = denominator->GetBinContent(ibin); | |
266 | // Check if the errors are right or the thing is scaled | |
267 | Double_t numeratorValErr = numerator->GetBinError(ibin); | |
268 | if (!(numeratorValErr==0. || numeratorVal ==0.) ) { | |
269 | Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal; | |
270 | numeratorVal /= rescale; | |
271 | } | |
272 | Double_t denominatorValErr = denominator->GetBinError(ibin); | |
273 | if (!(denominatorValErr==0. || denominatorVal==0. )) { | |
274 | Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal; | |
275 | denominatorVal /= rescale; | |
276 | } | |
277 | Double_t quotient = 0.; | |
278 | if (denominatorVal!=0.) { | |
279 | quotient = numeratorVal/denominatorVal; | |
280 | } | |
281 | Double_t quotientErr=0; | |
282 | quotientErr = TMath::Sqrt( | |
283 | (numeratorVal+1.0)/(denominatorVal+2.0)* | |
284 | ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))); | |
285 | result->SetBinContent(ibin,quotient); | |
286 | result->SetBinError(ibin,quotientErr); | |
287 | //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl; | |
288 | } | |
289 | return result; | |
290 | } | |
291 | ||
292 | ||
293 | TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) | |
294 | { | |
295 | if(!numerator){ | |
296 | cerr<<"Error: numerator does not exist!"<<endl; | |
297 | return NULL; | |
298 | } | |
299 | if(!denominator){ | |
300 | cerr<<"Error: denominator does not exist!"<<endl; | |
301 | return NULL; | |
302 | } | |
303 | TH2* result = (TH2*)numerator->Clone(name); | |
304 | Int_t nbinsX = numerator->GetNbinsX(); | |
305 | Int_t nbinsY = numerator->GetNbinsY(); | |
306 | for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) { | |
307 | for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) { | |
308 | Double_t numeratorVal = numerator->GetBinContent(ibin,jbin); | |
309 | Double_t denominatorVal = denominator->GetBinContent(ibin,jbin); | |
310 | // Check if the errors are right or the thing is scaled | |
311 | Double_t numeratorValErr = numerator->GetBinError(ibin,jbin); | |
312 | if (!(numeratorValErr==0. || numeratorVal ==0.) ) { | |
313 | Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal; | |
314 | numeratorVal /= rescale; | |
315 | } | |
316 | Double_t denominatorValErr = denominator->GetBinError(ibin,jbin); | |
317 | if (!(denominatorValErr==0. || denominatorVal==0. )) { | |
318 | Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal; | |
319 | denominatorVal /= rescale; | |
320 | } | |
321 | Double_t quotient = 0.; | |
322 | if (denominatorVal!=0.) { | |
323 | quotient = numeratorVal/denominatorVal; | |
324 | } | |
325 | Double_t quotientErr=0; | |
326 | quotientErr = TMath::Sqrt( | |
327 | (numeratorVal+1.0)/(denominatorVal+2.0)* | |
328 | ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))); | |
329 | result->SetBinContent(ibin,jbin,quotient); | |
330 | result->SetBinError(ibin,jbin,quotientErr); | |
331 | //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl; | |
332 | } | |
333 | } | |
334 | return result; | |
335 | } |