TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / totEt / macros / calculateCorrections.C
CommitLineData
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
18TCanvas *c1 = 0;
19
47151f26 20TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
21TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
9f3689ba 22
23//creates an empty set of correction factors for debugging purposes
943dc129 24TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
9f3689ba 25int 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 47int 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
242TF1* 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
251TH1* 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
293TH2* 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}