TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / totEt / macros / calculateCorrections.C
1 //Oystein Djuvsland
2 //macro for plotting correction factors for EM et for use with AliAnalysisEt and writes them to a file
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
20 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
21 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
22
23 //creates an empty set of correction factors for debugging purposes
24 TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
25 int createDummy(Double_t p0 = 1.0, Double_t p1 = 0.0, char *det = "Phos")
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);
32   AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(Form("TmCorrections%s",det), fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
33   TF1 *func = generateRecEffFunction(p0, p1);
34   AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(Form("ReCorrections%s",det), *func, 1000);
35
36   cor->Write();
37   recor->Write();
38   outfile->Close();
39
40
41 }
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
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)
48 {
49   TString detector = det;
50   TString emcal = "Emcal";
51   c1 = new TCanvas;
52   
53   TFile *f = TFile::Open(filename, "READ");
54   
55   TList *l = (TList*)f->Get("out1");
56   
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());
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   
77   Float_t maxMult = 99.5;
78   Int_t nbins = 100;
79    if(detector==emcal){
80 //     //100 seems to be sufficient...
81      maxMult = 249.5;
82      nbins = 250;
83   }
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
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);
114   TF1 fitcharged("fitcharged","pol2", 0, 100);//fit of number of charged tracks vs detector multiplicity
115   //if straight line track matching roughly not dependent on centrality
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");
139   TF1 fitneutral("fitneutral","pol2", 0, 100);//fit of number of neutral particles in calo that we call background in calo vs detector multiplicity
140   //may include K0S
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");
163   TF1 fitgamma("fitgamma","pol2", 0, 100);//fit of number of gammas removed erroneously vs detector multiplicity
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");
187   TF1 fitsecondary("fitsecondary","pol2", 0, 100);//fit of number of secondary particles that leave deposits in calo vs detector multiplicity
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   }
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
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   
206   TH2F  *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllCent");
207   TH2F  *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoCent");
208   TH2F *eff2D;
209 //   if(fHistHadronDepositsRecoMult && fHistHadronDepositsAllMult ){
210     eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D");
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 //   }
216   //cor->SetReconstructionEfficiency(eff2D);
217   
218   TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
219   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
220   TH2F *gammaEff2D;
221   //if(fHistGammasGeneratedMult && fHistGammasFoundMult){
222     gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D");
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 //   }
228
229   AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D,
230                                                                                    meanCharged, meanNeutral, meanGamma, meanSecondary );
231   
232
233   TF1 *func = generateRecEffFunction(p0, p1);
234   AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000);
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 }
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 }