Removing PWG3TrackExtrapInterface from PWG3base library (Gines)
[u/mrichter/AliRoot.git] / PWG3 / AliD0toKpiPlots.C
1 void AliD0toKpiPlots(const Char_t *inName="AliD0toKpi.root",
2                      const Char_t *outName="D0histograms.root") {
3   //--------------------------------------------------------------------------
4   // This macro histograms many variables of D0->Kpi candidates
5   //
6   //     Andrea Dainese, andrea.dainese@lnl.infn.it
7   //--------------------------------------------------------------------------
8
9   gSystem->Load("libAOD.so");
10   gSystem->Load("libPWG3base.so");
11
12   // set of cuts
13   Double_t D0Cuts[9] = {0.1,         // mass [GeV]
14                         1000000.,          // dca [micron]
15                         1.1,           // cosThetaStar
16                         0.,           // pT K [GeV/c]
17                         0.,           // pT Pi [GeV/c]    
18                         100000.,       // d0K upper [micron]
19                         100000.,       // d0Pi upper [micron]
20                         10000000000.,            // d0d0 [micron^2]
21                         -1.1};          // cosThetaPointing
22
23   // number of events (for normalization)
24   Bool_t normalize = kFALSE;
25   Double_t  events = 1.;
26
27
28
29   // define histograms
30   TH1F *hptK = new TH1F("hptK","\"K\" p_{t} distribution",50,0,10);
31   hptK->SetXTitle("p_{t} [GeV]");
32
33   TH1F *hptPi = new TH1F("hptPi","\"#pi\" p_{t} distribution",50,0,10);
34   hptPi->SetXTitle("p_{t} [GeV]");
35
36   TH1F *hDCA = new TH1F("hDCA","DCA",50,0,1000);
37   hDCA->SetXTitle("dca [#mu m]");
38
39   TH1F *hptD0 = new TH1F("hptD0","D^{0} p_{t} distribution",40,0,40);
40   hptD0->SetXTitle("p_{t} [GeV]");
41
42   TH1F *hyD0 = new TH1F("hyD0","D^{0} rapidity distribution",50,-2,2);
43   hyD0->SetXTitle("y");
44
45   TH1F *hCPtaD0 = new TH1F("hCPtaD0","cosine of pointing angle distribution",100,-1,1);
46   hCPtaD0->SetXTitle("cos #theta_{point}");
47
48   TH1F *hCPtaXY = new TH1F("hCPtaXY","cosine of pointing angle in (x,y) plane",100,-1,1);
49   hCPtaXY->SetXTitle("cos #theta_{point}");
50
51   TH1F *hCts = new TH1F("hCts","cosine of decay angle",50,-1.2,1.2);
52   hCts->SetXTitle("cos #theta^{*}");
53
54   TH2F *hCtsVsPtK = new TH2F("hCtsVsPtK","cosine of decay angle VS \"K\" p_{t}",50,0,5,50,-1,1);
55   hCtsVsPtK->SetYTitle("cos #theta^{*}");
56   hCtsVsPtK->SetXTitle("p_{t} [GeV]");
57
58   TH1F *hd0d0 = new TH1F("hd0d0","Product of the impact parameters",100,-100000,100000);
59   hd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
60
61   TH1F *hd0K = new TH1F("hd0K","Impact parameter of \"K\"",100,-5000,5000);
62   hd0K->SetXTitle("d_{0}^{K} [#mu m]");
63
64   TH1F *hd0Pi = new TH1F("hd0Pi","Impact parameter of \"#pi\"",100,-5000,5000);
65   hd0Pi->SetXTitle("d_{0}^{#pi} [#mu m]");
66
67   TH2F *hCPtaVsd0d0 = new TH2F("hCPtaVsd0d0","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,100000,100,-1,1);
68   hCPtaVsd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
69   hCPtaVsd0d0->SetYTitle("cos #theta_{point}");
70
71   TH2F *hCPtaVsd0d0zoom = new TH2F("hCPtaVsd0d0zoom","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,0,100,.9,1);
72   hCPtaVsd0d0zoom->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
73   hCPtaVsd0d0zoom->SetYTitle("cos #theta_{point}");
74
75   TH2F *hd0d0VSptD0 = new TH2F("hd0d0VSptD0","d_{0}^{K} #times d_{0}^{#pi} VS D^{0} p_{t}",50,0,25,100,-120000,120000);
76   hd0d0VSptD0->SetYTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
77   hd0d0VSptD0->SetXTitle("D^{0} p_{t} [GeV]");
78
79   TH1F *hMass = new TH1F("hMass","Invariant mass distribution",50,1.765,1.965);
80   hMass->SetXTitle("M[K,#pi] [GeV]");
81
82   TH2F *hArm = new TH2F("hArm","Armenteros plot",50,-2,2,50,0,1);
83   hArm->SetXTitle("#alpha");
84   hArm->SetYTitle("q_{t}");
85
86   // open input file and get tree
87   TFile *inFile = TFile::Open(inName);
88
89   TTree *treeD0 = (TTree*)inFile->Get("TreeD0");
90   AliD0toKpi *D = 0; 
91   treeD0->SetBranchAddress("D0toKpi",&D);
92   Int_t entries = (Int_t)treeD0->GetEntries();
93
94   printf("+++\n+++ Number of D0 in tree:  %d\n+++\n",entries);
95
96   Double_t MD0,MD0bar,ctsD0,ctsD0bar,ctsPiD0,ctsPiD0bar;
97   Double_t WgtD0,WgtD0bar;
98   Double_t sampleABC=0.;
99   Int_t okD0=0,okD0bar=0;
100   Int_t nSel = 0;
101   Int_t ptbin;
102
103   // loop on D0
104   for(Int_t i=0; i<entries; i++) {
105     if(i%10000==0) printf(" candidate %d of %d\n",i,entries);
106
107     // get event from tree
108     treeD0->GetEvent(i);
109     //--- select the PID strategy & compute weights
110     //    D->ApplyPID("TOFparam_PbPb");
111     //    D->ComputeWgts();
112     // get weights for the three samples A+B+C 
113     //    D->GetWgts(WgtD0,WgtD0bar,"ABC");
114     WgtD0 = 1.; WgtD0bar = 1.;
115
116     // normalize to 1 event
117     if(normalize) { WgtD0 /= events; WgtD0bar /= events; }
118
119     // check if candidate passes selection (as D0 or D0bar)
120     D->Select(D0Cuts,okD0,okD0bar);
121
122     // set weights to 0 if the candidate doesn't pass selection
123     if(!okD0)    WgtD0=0.; 
124     if(!okD0bar) WgtD0bar=0.;
125     if(okD0 || okD0bar) nSel++;
126
127     // count selected candidates
128     sampleABC += WgtD0 + WgtD0bar;
129
130     // inv mass and cosThetaStar
131     D->InvMass(MD0,MD0bar);
132     D->CosThetaStar(ctsD0,ctsD0bar);
133     
134     // fill histograms
135     hptK->Fill(D->PtChild(1),WgtD0);
136     hptK->Fill(D->PtChild(0),WgtD0bar);
137     hptPi->Fill(D->PtChild(0),WgtD0);
138     hptPi->Fill(D->PtChild(1),WgtD0bar);
139     hd0K->Fill(D->Getd0Child(1),WgtD0);
140     hd0K->Fill(D->Getd0Child(0),WgtD0bar);
141     hd0Pi->Fill(D->Getd0Child(0),WgtD0);
142     hd0Pi->Fill(D->Getd0Child(1),WgtD0bar);   
143     hMass->Fill(MD0,WgtD0);
144     hMass->Fill(MD0bar,WgtD0bar);
145     hCts->Fill(ctsD0,WgtD0);
146     hCts->Fill(ctsD0bar,WgtD0bar);
147     hCtsVsPtK->Fill(D->PtChild(1),ctsD0,WgtD0);
148     hCtsVsPtK->Fill(D->PtChild(0),ctsD0bar,WgtD0bar);
149     hDCA->Fill(D->GetDCA(),WgtD0+WgtD0bar);
150     hptD0->Fill(D->Pt(),WgtD0+WgtD0bar);
151     hyD0->Fill(D->Rapidity(),WgtD0+WgtD0bar);
152     hd0d0->Fill(D->ProdImpParams(),WgtD0+WgtD0bar);
153     hCPtaD0->Fill(D->CosPointing(),WgtD0+WgtD0bar);
154     hCPtaXY->Fill(D->CosPointingXY(),WgtD0+WgtD0bar);
155     hCPtaVsd0d0->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar);
156     hd0d0VSptD0->Fill(D->Pt(),D->ProdImpParams(),WgtD0+WgtD0bar);
157     hCPtaVsd0d0zoom->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar);
158     hArm->Fill(D->Alpha(),D->Qt(),WgtD0+WgtD0bar);
159     
160    
161   } // end loop on D0 candidates
162
163   inFile->Close();
164
165   printf("\n\n --- Total number of candidates passing selection: %d\n\n --- Sum of weights sample A+B+C: %f\n\n",nSel,sampleABC); 
166   
167   // draw histograms
168   TCanvas *c1 = new TCanvas("c1","pt K & pi",0,0,700,700);
169   c1->SetLogy(); 
170   hptK->Draw();
171   hptPi->Draw("same");
172
173   TCanvas *c2 = new TCanvas("c2","pt D0",0,0,700,700);
174   c2->SetLogy(); 
175   hptD0->Draw();
176
177   TCanvas *c3 = new TCanvas("c3","rapidity D0",0,0,700,700);
178   hyD0->Draw(); 
179
180   TCanvas *c4 = new TCanvas("c4","pointing angle",0,0,700,700);
181   hCPtaD0->Draw();
182
183   TCanvas *c5 = new TCanvas("c5","d0 x d0",0,0,700,700);
184   c5->SetLogy();
185   hd0d0->Draw();
186
187   TCanvas *c6 = new TCanvas("c6","pointing angle VS d0d0",0,0,700,700);
188   c6->SetLogz();
189   hCPtaVsd0d0->Draw("box");
190
191   TCanvas *c7 = new TCanvas("c7","mass",0,0,700,700);
192   hMass->Draw();
193
194   TCanvas *c8 = new TCanvas("c8","armenteros",0,0,700,700);
195   hArm->Draw("box");
196
197   TCanvas *c9 = new TCanvas("c9","decay angle",0,0,700,700); 
198   hCts->Draw();
199
200   TCanvas *c10 = new TCanvas("c10","dca",0,0,700,700);
201   c10->SetLogy();
202   hDCA->Draw();
203
204   TCanvas *c11 = new TCanvas("c11","d0 K & pi",0,0,700,700);
205   c11->SetLogy();
206   hd0K->Draw();
207   hd0Pi->Draw("same");
208
209   // write all histograms to file
210   TFile *outFile = new TFile(outName,"recreate");
211   hMass->Write();
212   hDCA->Write(); 
213   hCts->Write();
214   hCtsVsPtK->Write();
215   hArm->Write();
216   hCPtaVsd0d0->Write();
217   hd0d0VSptD0->Write();
218   hCPtaVsd0d0zoom->Write();
219   hd0d0->Write();
220   hCPtaD0->Write();
221   hCPtaXY->Write();
222   hptK->Write();
223   hptPi->Write();
224   hptD0->Write();
225   hyD0->Write();
226   hd0K->Write();
227   hd0Pi->Write();
228   outFile->Close();
229   
230   return;
231 }
232
233
234
235