.so cleanup: more gSystem->Load()
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / DrawEfficiency.C
1
2 // macro to calculate the pt RecAnCuts/MCAcc
3 // efficiency, from a CF container 
4 // Author: C. Zampolli 
5
6 // channel could be:
7 // D0: D0 -> Kpi, from old task
8 // D0_New: D0 -> Kpi, from new CF common implementation
9 // Dplus_New: D+ -> Kpipi, from new CF common implementation
10 //
11 // eff = sum of the efficiency indexes to compute
12 // MCAcc_over_MCLimAcc = 0x001
13 // RecPPR_over_MCAcc = 0x002
14 // RecPID_over_MCAcc = 0x004
15 // all = 0x007
16 //
17 // selection = D origin
18 // 0 --> from c only
19 // 1 --> from b only
20 // 2 --> from both c and b
21
22 #include <Riostream.h>
23
24 extern TRandom *gRandom;
25 extern TBenchmark *gBenchmark;
26 extern TSystem *gSystem;
27
28 void DrawEfficiency(const char* channel, Int_t selection = 0, Int_t ieff = 7){
29
30         gROOT->SetStyle("Plain");
31         gStyle->SetPalette(1);
32         gStyle->SetOptStat(0);
33         gStyle->SetPalette(1);
34         gStyle->SetCanvasColor(0);
35         gStyle->SetFrameFillColor(0);
36         gStyle->SetOptTitle(0);
37         
38         gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include  -I$ROOTSYS/include");
39         gSystem->Load("libANALYSIS");
40         gSystem->Load("libANALYSISalice");
41         gSystem->Load("$ALICE_ROOT/CORRFW/libCORRFW") ;
42         gSystem->Load("libPWGHFbase");
43         gSystem->Load("libPWGHFvertexingHF");
44                 
45         Int_t mcAcc_over_mcLimAcc = 0x001;
46         Int_t recPPR_over_mcAcc = 0x002;
47         Int_t recPID_over_mcAcc = 0x004;
48
49         // pt index
50         Int_t ipt =0;
51
52         Int_t stepNum = -1;
53         Int_t stepDen = -1;
54         
55         // Read the  container from file
56         TFile* f = new TFile("AnalysisResults.root");   
57         TString directoryName;
58         TString containerName;
59         TString cutName;
60         TString outfileName;
61
62         if (channel == "D0") {
63                 if (selection == 0){
64                         directoryName = "PWG3_D2H_CFtaskD0toKpi";
65                         containerName = "CFHFccontainer0";
66                         cutName = "Cuts";
67                         outfileName = "fileEff_D0_from_c.root";
68                 }
69                 else if (selection == 1){
70                         directoryName = "PWG3_D2H_CFtaskD0toKpiKeepD0fromBOnly";
71                         containerName = "CFHFccontainer0D0fromB";
72                         cutName = "Cuts";
73                         outfileName = "fileEff_D0_from_b.root";
74                 }
75                 else if (selection == 2){
76                         directoryName = "PWG3_D2H_CFtaskD0toKpiKeepD0fromB";
77                         containerName = "CFHFccontainer0allD0";
78                         cutName = "Cuts";
79                         outfileName = "fileEff_D0_from_c_and_b.root";
80                 }
81                 else {
82                         Printf("not a valid selection, return");
83                         return;
84                 }
85         }
86         else if (channel == "D0_New"){
87                 if (selection == 0){
88                         directoryName = "PWG3_D2H_CFtaskD0toKpi_NEW";
89                         containerName = "CFHFccontainer0_New";
90                         cutName = "Cuts_New";
91                         //directoryName = "PWG3_D2H_CFtaskD0toKpi_CommonFramework";
92                         //containerName = "CFHFccontainer0_CommonFramework";
93                         //cutName = "Cuts_CommonFramework";
94                         outfileName = "fileEff_D0_CommonFramework_from_c.root";
95                 }
96                 else if (selection == 1){
97                         directoryName = "PWG3_D2H_CFtaskD0toKpiKeepDfromBOnly";
98                         containerName = "CFHFccontainer0DfromB_New";
99                         cutName = "Cuts_New";
100                         //directoryName = "PWG3_D2H_CFtaskD0toKpiKeepDfromBOnly_CommonFramework";
101                         //containerName = "CFHFccontainer0DfromB_CommonFramework";      
102                         //cutName = "Cuts_CommonFramework";
103                         outfileName = "fileEff_D0_CommonFramework_from_b.root";
104                 }
105                 else if (selection == 2){
106                         directoryName = "PWG3_D2H_CFtaskD0toKpiKeepDfromB_NEW";
107                         containerName = "CFHFccontainer0allD_New";
108                         cutName = "Cuts_New";
109                         //directoryName = "PWG3_D2H_CFtaskD0toKpiKeepDfromB_CommonFramework";
110                         //containerName = "CFHFccontainer0allD_CommonFramework";
111                         //cutName = "Cuts_CommonFramework";
112                         outfileName = "fileEff_D0_CommonFramework_from_c_and_b.root";
113                 }
114                 else {
115                         Printf("not a valid selection, return");
116                         return;
117                 }
118         }
119         else if (channel == "Dplus_New"){
120                 if (selection == 0){
121                         directoryName = "PWG3_D2H_CFtaskDplustoKpipi_NEW";
122                         containerName = "CFHFccontainer0_New_3Prong";
123                         cutName = "Cuts_3Prong";
124                         //directoryName = "PWG3_D2H_CFtaskDplustoKpipi_CommonFramework";
125                         //containerName = "CFHFccontainer0_3Prong_CommonFramework";
126                         //cutName =  "Cuts_3Prong_CommonFramework";
127                         outfileName = "fileEff_Dplus_CommonFramework_from_c.root";
128                 }
129                 else if (selection == 1){
130                         directoryName = "PWG3_D2H_CFtaskDplustoKpipiKeepDfromBOnly";
131                         containerName = "CFHFccontainer0DfromB_New_3Prong";
132                         cutName = "Cuts_3Prong";
133                         //directoryName = "PWG3_D2H_CFtaskDplustoKpipiKeepDfromBOnly_CommonFramework";
134                         //containerName = "CFHFccontainer0DfromB_3Prong_CommonFramework";
135                         //cutName =  "Cuts_3Prong_CommonFramework";
136                         outfileName = "fileEff_Dplus_CommonFramework_from_b.root";
137                 }
138                 else if (selection == 2){
139                         directoryName = "PWG3_D2H_CFtaskDplustoKpipiKeepDfromB_NEW";
140                         containerName = "CFHFccontainer0allD_New_3Prong";
141                         cutName = "Cuts_3Prong";
142                         //directoryName = "PWG3_D2H_CFtaskDplustoKpipiKeepDfromB_CommonFramework";
143                         //containerName = "CFHFccontainer0allD_3Prong_CommonFramework";
144                         //cutName =  "Cuts_3Prong_CommonFramework";
145                         outfileName = "fileEff_Dplus_CommonFramework_from_c_and_b.root";
146                 }
147                 else {
148                         Printf("not a valid selection, return");
149                         return;
150                 }
151         }
152         else {
153                 Printf("not a valid channel, return");
154                 return;
155         }
156
157         Printf("Opening file Analysisresults.root");
158         Printf("Reading Directory \"%s\"",directoryName.Data());
159         Printf("Getting CF Container \"%s\"",containerName.Data());
160         Printf("Getting Cut Object \"%s\"",cutName.Data());
161
162
163         TDirectoryFile* d = (TDirectoryFile*)f->Get(directoryName.Data());
164         if (!d){
165                 Printf("Directory does not exist! Check directory name (%s) in file AnalysisResults.root - returning...", directoryName.Data());
166                 return;
167         }
168         AliCFContainer *data = (AliCFContainer*) (d->Get(containerName.Data()));
169         AliRDHFCuts *cutsRDHF = (AliRDHFCuts*)(d->Get(cutName.Data()));
170
171         if (!data){
172                 Printf("Container does not exist! Check container name (%s) in directory %s - returning...", containerName.Data(), directoryName.Data());
173                 return;
174         }
175         
176         TFile* fileEff = new TFile(outfileName.Data(), "RECREATE");
177         TString plotDir(Form("EffPlots/%s",channel));
178         gSystem->Exec(Form("mkdir EffPlots"));
179         gSystem->Exec(Form("mkdir %s",plotDir.Data()));
180         
181         //construct the efficiency grid from the data container 
182         AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
183
184         TCanvas *ceffpt = new TCanvas("ceffpt","Efficiency vs pt",50,50,550,550);
185         ceffpt->cd();
186         ceffpt->SetLeftMargin(0.15);
187         ceffpt->SetRightMargin(0.05);
188         TH1D *hpteffCF = 0x0; //the efficiency vs pt
189
190         if (ieff & mcAcc_over_mcLimAcc){
191                 AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
192                 stepDen = (Int_t)(AliCFTaskVertexingHF::kStepGeneratedLimAcc);  
193                 stepNum = (Int_t)(AliCFTaskVertexingHF::kStepAcceptance);       
194                 printf("Calculating efficiency for mcAcc_over_mcLimAcc: stepDen = %d, stepNum = %d\n",stepDen,stepNum); 
195                 eff->CalculateEfficiency(stepNum,stepDen); //eff= step1/step0
196                 
197                 //canvas
198                 ceffpt->cd();
199
200                 //The efficiency along the variables
201                 hpteffCF = eff->Project(ipt); 
202                 SetHistoEff(hpteffCF,8,20,"mcAcc_over_mcLimAcc");
203                 hpteffCF->Draw("hist");
204                 hpteffCF->Draw("err same");
205                 fileEff->cd();
206                 hpteffCF->Write("hpteffCF_mcAcc_over_mcLimAcc");
207                 
208                 // printing png files
209                 ceffpt->Print(Form("%s/effpt_mcAcc_over_mcLimAcc.png", plotDir.Data()));
210                 ceffpt->Print(Form("%s/effpt_mcAcc_over_mcLimAcc.eps", plotDir.Data()));
211                 ceffpt->Print(Form("%s/effpt_mcAcc_over_mcLimAcc.gif", plotDir.Data()));
212                 delete eff;
213                 eff = 0x0;
214         }
215
216         if (ieff & recPPR_over_mcAcc){
217                 AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
218                 stepDen = (Int_t)(AliCFTaskVertexingHF::kStepAcceptance);       
219                 stepNum = (Int_t)(AliCFTaskVertexingHF::kStepRecoPPR);  
220                 printf("Calculating efficiency for RecPPR_over_mcAcc: stepDen = %d, stepNum = %d\n",stepDen,stepNum);   
221                 eff->CalculateEfficiency(stepNum,stepDen); //eff= step1/step0
222                 
223                 //canvas
224                 ceffpt->cd();
225                 
226                 //The efficiency along the variables
227                 hpteffCF = eff->Project(ipt); 
228                 SetHistoEff(hpteffCF,8,20, "recAnCuts_over_mcAcc");
229                 hpteffCF->Draw("hist");
230                 hpteffCF->Draw("err same");
231                 fileEff->cd();
232                 hpteffCF->Write("hpteffCF_RecAnCut_over_mcAcc");
233                 
234                 // printing png files
235                 ceffpt->Print(Form("%s/effpt_RecAnCut_over_mcAcc.png", plotDir.Data()));
236                 ceffpt->Print(Form("%s/effpt_RecAnCut_over_mcAcc.eps", plotDir.Data()));
237                 ceffpt->Print(Form("%s/effpt_RecAnCut_over_mcAcc.gif", plotDir.Data()));
238                 delete eff;
239                 eff = 0x0;
240         }
241
242         if (ieff & recPID_over_mcAcc){
243                 AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
244                 stepDen = (Int_t)(AliCFTaskVertexingHF::kStepAcceptance);       
245                 stepNum = (Int_t)(AliCFTaskVertexingHF::kStepRecoPID);  
246                 printf("Calculating efficiency for RecPID_over_mcAcc: stepDen = %d, stepNum = %d\n",stepDen,stepNum);   
247                 eff->CalculateEfficiency(stepNum,stepDen); //eff= step1/step0
248                 
249                 //canvas
250                 ceffpt->cd();
251                 
252                 //The efficiency along the variables
253                 hpteffCF = eff->Project(ipt); 
254                 SetHistoEff(hpteffCF,8,20,"recPID_over_mcAcc");
255                 hpteffCF->Draw("hist");
256                 hpteffCF->Draw("err same");
257                 fileEff->cd();
258                 hpteffCF->Write("hpteffCF_RecPID_over_mcAcc");
259                 
260                 // printing png files
261                 ceffpt->Print(Form("%s/effpt_RecPID_over_mcAcc.png", plotDir.Data()));
262                 ceffpt->Print(Form("%s/effpt_RecPID_over_mcAcc.eps", plotDir.Data()));
263                 ceffpt->Print(Form("%s/effpt_RecPID_over_mcAcc.gif", plotDir.Data()));
264                 delete eff;
265                 eff = 0x0;
266         }
267         
268         cutsRDHF->Write("Cuts");
269
270         // writing single distributions
271         TH1D *hMCAccpt = data->ShowProjection(ipt, AliCFTaskVertexingHF::kStepAcceptance);
272         SetHistoDistribution(hMCAccpt, 1, 20);
273         hMCAccpt->Draw();
274         TH1D *hMCLimAccpt = data->ShowProjection(ipt, AliCFHeavyFlavourTaskMultiVarMultiStep::kStepGeneratedLimAcc);
275         SetHistoDistribution(hMCLimAccpt, 4, 20);
276         TH1D *hRecoAnCutspt = data->ShowProjection(ipt, AliCFTaskVertexingHF::kStepRecoPPR);
277         SetHistoDistribution(hRecoAnCutspt, 8, 20);
278         TH1D *hRecoPIDpt = data->ShowProjection(ipt, AliCFTaskVertexingHF::kStepRecoPID);
279         SetHistoDistribution(hRecoPIDpt, 6, 20);
280         hMCAccpt->Write("hMCAccpt");
281         hMCLimAccpt->Write("hMCLimAccpt");
282         hRecoAnCutspt->Write("hRecoAnCutspt");
283         hRecoPIDpt->Write("hRecoPIDpt");
284
285         //      fileEff->Close(); // commented out to see the canvas on the screen....
286
287 }
288
289 void SetHistoEff(TH1D* h, Int_t color, Int_t style, const char* effType){
290
291         h->SetLineColor(color);
292         h->SetLineWidth(3);
293         h->SetMarkerStyle(style);
294         h->SetMarkerColor(color);
295         h->SetMarkerSize(1.2);
296         h->GetYaxis()->SetTitleOffset(1.5);
297         h->GetXaxis()->SetTitleOffset(1.2);
298         h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
299         h->GetYaxis()->SetTitle(Form("%s, Efficiency",effType));
300         return;
301 }
302 void SetHistoDistribution(TH1D* h, Int_t color, Int_t style){
303
304         h->SetLineColor(color);
305         h->SetLineWidth(3);
306         h->SetMarkerStyle(style);
307         h->SetMarkerColor(color);
308         h->SetMarkerSize(1.2);
309         h->GetXaxis()->SetTitleOffset(1.2);
310         h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
311         return;
312 }