AliPhysics  608b256 (608b256)
EstimateDmesonPIDsyst.C
Go to the documentation of this file.
1 #if !defined (__CINT__) || defined (__CLING__)
2 
3 #include <Riostream.h>
4 #include <TFile.h>
5 #include <TDirectoryFile.h>
6 #include <TCanvas.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TStyle.h>
10 #include <TString.h>
11 #include <TList.h>
12 #include <TMath.h>
13 #include <TGaxis.h>
14 
16 
17 #endif
18 
20 // \brief macro for estimation of PID Systematic uncertainty propagated from the single track to the D mesons //
21 // \main function: EstimateDmesonPIDsyst //
22 // \author: A. M. Barbano, anastasia.maria.barbano@cern.ch //
23 // \author: F. Grosa, fabrizio.grosa@cern.ch //
25 
26 using namespace std;
27 
28 //_____________________________________________________
29 //METHOD PROTOTYPES
32  TString infilename="LHC18a4a2_prop.root",
33  TString suffix="ppMB_kINT7",
34  double ptmin=2.,
35  double ptmax=24.);
36 void SetStyle();
37 
38 //_____________________________________________________
39 //METHOD IMPLEMENTATIONS
40 int EstimateDmesonPIDsyst(int ch, int pid, TString infilename, TString suffix, double ptmin, double ptmax) {
41 
42  SetStyle();
43 
44  TString mesonname = "";
45  if(ch == AliAnalysisTaskSEDmesonPIDSysProp::kD0toKpi) mesonname = "D0";
46  else if(ch == AliAnalysisTaskSEDmesonPIDSysProp::kDplustoKpipi) mesonname = "Dplus";
47  else if(ch == AliAnalysisTaskSEDmesonPIDSysProp::kDstartoKpipi) mesonname = "Dstar";
48  else if(ch == AliAnalysisTaskSEDmesonPIDSysProp::kDstoKKpi) mesonname = "Ds";
49 
50  TString PIDname = "";
51  double ymax = 0.05;
52  if(pid == AliAnalysisTaskSEDmesonPIDSysProp::kConservativePID) PIDname = "conservativePID";
54  PIDname = "strongPID";
55  ymax = 0.1;
56  }
57  else if(pid == AliAnalysisTaskSEDmesonPIDSysProp::knSigmaPID) PIDname = "nSigmaPID";
58 
59  TFile* infile = TFile::Open(infilename.Data());
60  if(!infile || !infile->IsOpen()) return 1;
61  TDirectoryFile* indir = (TDirectoryFile*)infile->Get(Form("PWGHF_D2H_PIDeffsyst_%s_%s_%s",mesonname.Data(),PIDname.Data(),suffix.Data()));
62  cout << Form("PWGHF_D2H_PIDeffsyst_%s_%s_%s",mesonname.Data(),PIDname.Data(),suffix.Data()) << endl;
63  if(!indir) {
64  cerr << "TDirectoryFile not found, please check the name. Exit" << endl;
65  return 2;
66  }
67  TList* inlist = (TList*)indir->Get(Form("systUnc_%s_%s_%s",mesonname.Data(),PIDname.Data(),suffix.Data()));
68  if(!inlist) {
69  cerr << "TList not found, please check the name. Exit" << endl;
70  return 3;
71  }
72  TH2F* hSystVPt = (TH2F*)inlist->FindObject("fHistSystPIDEffD");
73  TH2F* hPtCorr = (TH2F*)inlist->FindObject("fHistPtDauVsD");
74  hSystVPt->SetDirectory(0);
75  hPtCorr->SetDirectory(0);
76  infile->Close();
77 
78  TH1F* hSyst = (TH1F*)hSystVPt->ProjectionX("hSyst");
79  hSyst->SetLineWidth(2);
80  hSyst->SetLineColor(kRed);
81  hSyst->SetMarkerColor(kRed);
82  hSyst->SetMarkerStyle(kFullCircle);
83  for(int iPt=0; iPt<hSyst->GetNbinsX(); iPt++) {
84  TH1F* hTmp = (TH1F*)hSystVPt->ProjectionY("hTmp",iPt+1,iPt+1);
85  hSyst->SetBinContent(iPt+1,hTmp->GetMean());
86  hSyst->SetBinError(iPt+1,hTmp->GetMeanError());
87  }
88 
89  hSystVPt->GetYaxis()->SetTitleOffset(1.4);
90  hSystVPt->GetXaxis()->SetTitleOffset(1.2);
91  hSystVPt->GetYaxis()->SetTitleSize(0.045);
92  hSystVPt->GetXaxis()->SetTitleSize(0.045);
93  hSystVPt->GetYaxis()->SetLabelSize(0.04);
94  hSystVPt->GetXaxis()->SetLabelSize(0.04);
95  hPtCorr->GetYaxis()->SetTitleOffset(1.4);
96  hPtCorr->GetXaxis()->SetTitleOffset(1.2);
97  hPtCorr->GetYaxis()->SetTitleSize(0.045);
98  hPtCorr->GetXaxis()->SetTitleSize(0.045);
99  hPtCorr->GetYaxis()->SetLabelSize(0.04);
100  hPtCorr->GetXaxis()->SetLabelSize(0.04);
101 
102  TCanvas* cSyst = new TCanvas("cSyst","",1000,500);
103  cSyst->Divide(2,1);
104  cSyst->cd(1)->DrawFrame(ptmin,0.,ptmax,ptmax,";#it{p}_{T}^{D} (GeV/#it{c});#it{p}_{T}^{daugh} (GeV/#it{c})");
105  cSyst->cd(1)->SetLogz();
106  hPtCorr->Draw("same colz");
107  cSyst->cd(2)->DrawFrame(ptmin,0.,ptmax,ymax,";#it{p}_{T}^{D} (GeV/#it{c});relative systematic uncertainty");
108  cSyst->cd(2)->SetLogz();
109  hSystVPt->Draw("colz same");
110  hSyst->Draw("same");
111 
112  TFile outfile(Form("%s_%s_PIDsyst.root",mesonname.Data(),PIDname.Data()),"recreate");
113  hPtCorr->Write();
114  hSystVPt->Write();
115  hSyst->Write();
116  outfile.Close();
117 
118  cSyst->SaveAs(Form("%s_%s_PIDsyst.pdf",mesonname.Data(),PIDname.Data()));
119  return 0;
120 }
121 
122 //_____________________________________________________
123 void SetStyle() {
124  gStyle->SetPadTopMargin(0.05);
125  gStyle->SetPadRightMargin(0.12);
126  gStyle->SetPadLeftMargin(0.14);
127  gStyle->SetPadBottomMargin(0.12);
128  gStyle->SetTitleOffset(1.4,"y");
129  gStyle->SetTitleOffset(1.2,"x");
130  gStyle->SetTitleSize(0.045,"xyz");
131  gStyle->SetLabelSize(0.04,"xyz");
132  gStyle->SetOptTitle(0);
133 
134  gStyle->SetPadTickX(1);
135  gStyle->SetPadTickY(1);
136  gStyle->SetOptStat(0);
137 
138  TGaxis::SetMaxDigits(3);
139 }
const Double_t ymax
Definition: AddTaskCFDStar.C:7
Definition: External.C:236
int EstimateDmesonPIDsyst(int ch=AliAnalysisTaskSEDmesonPIDSysProp::kDstoKKpi, int pid=AliAnalysisTaskSEDmesonPIDSysProp::kStrongPID, TString infilename="LHC18a4a2_prop.root", TString suffix="ppMB_kINT7", double ptmin=2., double ptmax=24.)
const Double_t ptmax
const Double_t ptmin
void SetStyle()