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.so");
  40         gSystem->Load("libANALYSISalice.so");
  41         gSystem->Load("$ALICE_ROOT/CORRFW/libCORRFW.so") ;
  42         gSystem->Load("libPWG3base.so");
  43         gSystem->Load("libPWG3vertexingHF.so");
  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 }