AliPhysics  7f2a7c4 (7f2a7c4)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ReadAODVertexingHF.C
Go to the documentation of this file.
1 void ReadAODVertexingHF(const char *aodFileName="AliAOD.root",
2  const char *aodHFFileName="AliAOD.VertexingHF.root")
3 {
4  //
5  // Example macro to read D0->Kpi candidates from AOD (having the
6  // standard AOD + a friend heavy-flavour AOD) and apply cuts
7  // Origin: A.Dainese
8  //
9 
10  TStopwatch t;
11  t.Start();
12 
13  Bool_t useParFiles=kFALSE;
14  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
15  LoadLibraries(useParFiles);
16 
17  // create a test histogram
18  TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
19  hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
20  hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
21  TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
22  hMass->SetXTitle("Invariant mass [GeV]");
23  hMass->SetYTitle("Entries");
24  TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
25  hSecVtxZ->SetXTitle("z of decay vertex [cm]");
26  hSecVtxZ->SetYTitle("Entries");
27  TH1F *hDeltaMassDstar = new TH1F("hDeltaMassDstar","D* delta mass plot",100,0,0.3);
28  hDeltaMassDstar->SetXTitle("M(Kpipi)-M(Kpi) [GeV]");
29  hDeltaMassDstar->SetYTitle("Entries");
30 
31 
32  // open input file and get the TTree
33  TFile inFile(aodFileName,"READ");
34  if (!inFile.IsOpen()) return;
35 
36  TTree *aodTree = (TTree*)inFile.Get("aodTree");
37  aodTree->AddFriend("aodTree",aodHFFileName);
38 
39  AliAODEvent *aod = new AliAODEvent();
40 
41  aod->ReadFromTree(aodTree);
42 
43  // load heavy flavour vertices
44  TClonesArray *arrayVerticesHF =
45  (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
46 
47  // load D0->Kpi candidates
48  TClonesArray *arrayD0toKpi =
49  (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
50 
51  // load 3prong candidates
52  TClonesArray *array3Prong =
53  (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
54 
55  // load D* candidates
56  TClonesArray *arrayDstar =
57  (TClonesArray*)aod->GetList()->FindObject("Dstar");
58 
59  // load cascade (V0+track) candidates
60  TClonesArray *arrayCascades =
61  (TClonesArray*)aod->GetList()->FindObject("CascadesHF");
62 
63 
64 
65  Int_t nTotHF=0,nTotD0toKpi=0,nTotDstar=0,nTot3Prong=0,nTotCasc=0;
66  AliAODVertex *vtx1=0;
67 
68  AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi("CutsD0toKpi");
69  cutsD0toKpi->SetStandardCutsPP2010();
70  cutsD0toKpi->SetRemoveDaughtersFromPrim(kFALSE);
71 
72  AliRDHFCutsDStartoKpipi *cutsDStar = new AliRDHFCutsDStartoKpipi("CutsDStartoKpipi");
73  cutsDStar->SetStandardCutsPP2010();
74 
75  AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi("CutsDplustoKpipi");
76  cutsDplustoKpipi->SetStandardCutsPP2010();
77  cutsDplustoKpipi->SetRemoveDaughtersFromPrim(kFALSE);
78 
79  Int_t nTot3ProngSele=0;
80  Int_t nTotD0toKpiSele=0;
81  Int_t nTotDStarSele=0;
82 
83  // loop over events
84  Int_t nEvents = aodTree->GetEntries();
85  for (Int_t nEv = 0; nEv < nEvents; nEv++) {
86  cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
87 
88  // read event
89  aodTree->GetEvent(nEv);
90  //aodTree->BranchRef();
91 
92  //print event info
93  aod->GetHeader()->Print();
94 
95  // primary vertex
96  vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
97  vtx1->Print();
98 
99  /*
100  // make trkIDtoEntry register (temporary)
101  Int_t trkIDtoEntry[100000];
102  for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
103  AliAODTrack *track = aod->GetTrack(it);
104  trkIDtoEntry[track->GetID()]=it;
105  }
106  */
107 
108  // Fix references to daughter tracks
109  //AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
110  //fixer->FixReferences(aod);
111  //delete fixer;
112  //
113  //AliRDHFCutsD0toKpi *mycuts=new AliRDHFCutsD0toKpi();
114  //mycuts->SetFixRefs(kTRUE);
115  //mycuts->IsEventSelected(aod);
116 
117  // loop over D0->Kpi candidates
118  Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
119  nTotD0toKpi += nD0toKpi;
120  cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
121 
122  for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
123  AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
124 
125  d->SetOwnPrimaryVtx(vtx1);
126 
127  /*
128  // get daughter AOD tracks
129  AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
130  AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
131 
132  if(trk0->GetStatus()) printf("ok %d\n",iD0toKpi);
133  */
134  printf("D0 %d\n",iD0toKpi);
135  if(cutsD0toKpi->IsSelected(d,AliRDHFCuts::kAll)) {printf("D0 %d passed\n",iD0toKpi); nTotD0toKpiSele++;}
136 
137  }
138 
139 
140  // loop over D* candidates
141  Int_t nDstar = arrayDstar->GetEntriesFast();
142  nTotDstar += nDstar;
143  cout<<"Number of D*->D0pi: "<<nDstar<<endl;
144 
145 
146 
147  for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
148  AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar);
149  printf("D* %d\n",iDstar);
150  if(cutsDStar->IsSelected(c,AliRDHFCuts::kCandidate)) {printf("D* %d passed\n",iDstar); nTotDStarSele++;}
151 
152  }
153 
154 
155 
156  // count 3prong candidates
157  Int_t n3Prong = array3Prong->GetEntriesFast();
158  nTot3Prong += n3Prong;
159  cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
160 
161  for (Int_t i3p = 0; i3p < n3Prong; i3p++) {
162  AliAODRecoDecayHF3Prong *ccc = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3p);
163  printf("3p %d\n",i3p);
164  //if(cutsDplustoKpipi->IsSelected(ccc,AliRDHFCuts::kCandidate)) {printf("3p %d passed\n",i3p); nTot3ProngSele++;}
165 
166  }
167 
168  /*
169  // loop over HF vertices
170  Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
171  nTotHF += nVtxsHF;
172  cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
173  for (Int_t iVtx = 0; iVtx <0; iVtx++) {
174  AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
175  // print info
176  //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
177  }
178 
179 
180  // count cascade candidates
181  if (arrayCascades){
182  Int_t nCasc = arrayCascades->GetEntriesFast();
183  nTotCasc+=nCasc;
184  cout << "Number of Cascades: "<<nCasc<<endl;
185  }
186  */
187  }
188 
189  printf("\n Total HF vertices: %d\n",nTotHF);
190  printf("\n Total D0->Kpi: %d; selected %d\n",nTotD0toKpi,nTotD0toKpiSele);
191  printf("\n Total D*->D0pi: %d; selected %d\n",nTotDstar,nTotDStarSele);
192  printf("\n Total Charm->3Prong: %d; selected %d\n",nTot3Prong,nTot3ProngSele);
193  if (arrayCascades) printf("\n Total Cascades: %d\n",nTotCasc);
194 
195  /*
196  TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
197  c1->Divide(2,2);
198  c1->cd(1);
199  hCPtaVSd0d0->Draw("colz");
200  c1->cd(2);
201  hMass->SetFillColor(4);
202  hMass->Draw();
203  c1->cd(3);
204  hSecVtxZ->SetFillColor(2);
205  hSecVtxZ->Draw();
206  c1->cd(4);
207  hDeltaMassDstar->SetFillColor(3);
208  hDeltaMassDstar->Draw();
209  */
210 
211  t.Stop();
212  t.Print();
213 
214  return;
215 }
216 //------------------------------------------------------------------------
217 void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
218 {
219  //
220  // Example macro to read D0->Kpi candidates from a stand-alone
221  // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
222  // Origin: A.Dainese
223  //
224  Bool_t useParFiles=kFALSE;
225  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
226  LoadLibraries(useParFiles);
227 
228  // create a test histogram
229  TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
230  hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
231  hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
232  TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
233  hMass->SetXTitle("Invariant mass [GeV]");
234  hMass->SetYTitle("Entries");
235  TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
236  hSecVtxZ->SetXTitle("z of decay vertex [cm]");
237  hSecVtxZ->SetYTitle("Entries");
238 
239  // open input file and get the TTree
240  TFile inFile(aodHFFileName,"READ");
241  if (!inFile.IsOpen()) return;
242 
243  TTree *aodTree = (TTree*)inFile.Get("aodTree");
244 
245  AliAODEvent *aod = new AliAODEvent();
246 
247  aod->ReadFromTree(aodTree);
248 
249  // load heavy flavour vertices
250  TClonesArray *arrayVerticesHF =
251  (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
252 
253  // load D0->Kpi candidates
254  TClonesArray *arrayD0toKpi =
255  (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
256 
257 
258  Double_t cutsD0[9]=
259  // cutsD0[0] = inv. mass half width [GeV]
260  // cutsD0[1] = dca [cm]
261  // cutsD0[2] = cosThetaStar
262  // cutsD0[3] = pTK [GeV/c]
263  // cutsD0[4] = pTPi [GeV/c]
264  // cutsD0[5] = d0K [cm] upper limit!
265  // cutsD0[6] = d0Pi [cm] upper limit!
266  // cutsD0[7] = d0d0 [cm^2]
267  // cutsD0[8] = cosThetaPoint
268  {1000.,
269  100000.,
270  1.1,
271  0.,
272  0.,
273  100000.,
274  100000.,
275  100000000.,
276  -1.1};
277 
278  Int_t nTotHF=0,nTotD0toKpi=0;
279  AliAODVertex *vtx1=0;
280 
281  // loop over events
282  Int_t nEvents = aodTree->GetEntries();
283  for (Int_t nEv = 0; nEv < nEvents; nEv++) {
284  cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
285 
286  // read event
287  aodTree->GetEvent(nEv);
288  //aodTree->BranchRef();
289 
290 
291  // loop over D0->Kpi candidates
292  Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
293  nTotD0toKpi += nD0toKpi;
294  cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
295 
296  for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
297  AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
298  Int_t okD0=0,okD0bar=0;
299  if(d->SelectD0(cutsD0,okD0,okD0bar)) {
300  //cout<<1e8*d->Prodd0d0()<<endl;
301  hMass->Fill(d->InvMassD0(),0.5);
302  hMass->Fill(d->InvMassD0bar(),0.5);
303  hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
304  hSecVtxZ->Fill(d->GetSecVtxZ());
305  //cout<<d->GetSecVtxZ() <<endl;
306 
307  }
308  }
309 
310  }
311 
312  printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
313 
314  TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
315  c->Divide(2,2);
316  c->cd(1);
317  hCPtaVSd0d0->Draw("colz");
318  c->cd(2);
319  hMass->SetFillColor(4);
320  hMass->Draw();
321  c->cd(3);
322  hSecVtxZ->SetFillColor(2);
323  hSecVtxZ->Draw();
324 
325  return;
326 }
double Double_t
Definition: External.C:58
Definition: External.C:236
virtual void SetStandardCutsPP2010()
TCanvas * c
Definition: TestFitELoss.C:172
Class for cuts on AOD reconstructed D+->Kpipi.
int Int_t
Definition: External.C:63
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
Bool_t SelectD0(const Double_t *cuts, Int_t &okD0, Int_t &okD0bar) const
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
void LoadLibraries()
Float_t nEvents[nProd]
void SetRemoveDaughtersFromPrim(Bool_t removeDaughtersPrim)
Definition: AliRDHFCuts.h:216
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
void ReadAODVertexingHF(const char *aodFileName="AliAOD.root", const char *aodHFFileName="AliAOD.VertexingHF.root")