Access to the ESD objects stored in a tree. Coding conventions (T.Kuhr)
[u/mrichter/AliRoot.git] / ANALYSIS / AliD0toKpiAnalysis.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //----------------------------------------------------------------------------
17 //    Implementation of the D0toKpi reconstruction and analysis class
18 // Note: the two decay tracks are labelled: 0 (positive track)
19 //                                          1 (negative track)
20 // An example of usage can be found in the macro AliD0toKpiTest.C
21 //            Origin: A. Dainese    andrea.dainese@pd.infn.it            
22 //----------------------------------------------------------------------------
23 #include <TKey.h>
24 #include <TFile.h>
25 #include <TTree.h>
26 #include <TString.h>
27 #include <TSystem.h>
28 #include <TParticle.h>
29 #include "AliESD.h"
30 #include "AliStack.h"
31 #include "AliRunLoader.h"
32 #include "AliITStrackV2.h"
33 #include "AliITSVertexerTracks.h"
34 #include "AliESDVertex.h"
35 #include "AliV0vertexer.h"
36 #include "AliV0vertex.h"
37 #include "AliD0toKpi.h"
38 #include "AliD0toKpiAnalysis.h"
39
40 typedef struct {
41   Int_t lab;
42   Int_t pdg;
43   Int_t mumlab;
44   Int_t mumpdg;
45   Float_t Vx,Vy,Vz;
46   Float_t Px,Py,Pz;
47 } REFTRACK;
48
49 ClassImp(AliD0toKpiAnalysis)
50
51 //----------------------------------------------------------------------------
52 AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
53   // Default constructor
54
55   SetBz();
56   SetPtCut();
57   Setd0Cut();
58   SetMassCut();
59   SetD0Cuts();
60   SetVertex1();
61   SetPID();
62   fVertexOnTheFly = kFALSE;
63   fSim = kFALSE;
64   fOnlySignal = kFALSE;
65   fDebug = kFALSE;
66 }
67 //----------------------------------------------------------------------------
68 AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
69 //----------------------------------------------------------------------------
70 void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
71   // select candidates that pass fD0Cuts and write them to a new file
72
73   TFile *inFile = TFile::Open(inName);
74
75   TTree *treeD0in=(TTree*)inFile->Get("TreeD0");
76   AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis");
77   printf("+++\n+++  I N P U T   S T A T U S:\n+++\n");
78   inAnalysis->PrintStatus();
79
80
81   AliD0toKpi *d = 0; 
82   treeD0in->SetBranchAddress("D0toKpi",&d);
83   Int_t entries = (Int_t)treeD0in->GetEntries();
84
85   printf("+++\n+++ Number of D0 in input tree:  %d\n+++\n",entries);
86
87   TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
88   treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
89
90
91   Int_t okD0=0,okD0bar=0;
92   Int_t nSel = 0;
93
94   for(Int_t i=0; i<entries; i++) {
95     // get event from tree
96     treeD0in->GetEvent(i);
97
98     if(fSim && fOnlySignal && !d->IsSignal()) continue; 
99
100     // check if candidate passes selection (as D0 or D0bar)
101     if(d->Select(fD0Cuts,okD0,okD0bar)) {
102       nSel++;
103       treeD0out->Fill();
104     }
105
106   }
107
108   AliD0toKpiAnalysis *outAnalysis = (AliD0toKpiAnalysis*)inAnalysis->Clone("D0toKpiAnalysis");
109   outAnalysis->SetD0Cuts(fD0Cuts);
110   printf("------------------------------------------\n");
111   printf("+++\n+++  O U T P U T   S T A T U S:\n+++\n");
112   outAnalysis->PrintStatus();
113
114   printf("+++\n+++ Number of D0 in output tree:  %d\n+++\n",nSel);
115
116   TFile* outFile = new TFile(outName,"recreate");
117   treeD0out->Write();
118   outAnalysis->Write();
119   outFile->Close();
120
121   return;
122 }
123 //----------------------------------------------------------------------------
124 Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
125                                               Double_t time) const {
126   // calculated the mass from momentum, track length from vertex to TOF
127   // and time measured by the TOF
128   if(length==0.) return -1000.;
129   Double_t a = time*time/length/length;
130   if(a > 1.) {
131     a = TMath::Sqrt(a-1.);
132   } else {
133     a = -TMath::Sqrt(1.-a);
134   }
135
136   return mom*a;
137 }
138 //----------------------------------------------------------------------------
139 void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
140                                         const Char_t *outName) {
141   // Find D0 candidates and calculate parameters
142
143   if(fBz<-9000.) {
144     printf("AliD0toKpiAnalysis::FindCandidates():  Set B!\n");
145     return;
146   }
147   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
148
149   TString trkName("AliITStracksV2.root");
150   if(gSystem->AccessPathName(trkName.Data(),kFileExists)) {
151     printf("AliD0toKpiAnalysis::FindCandidates(): No tracks file found!\n"); 
152     return;
153   }
154
155   TString outName1=outName;
156   TString outName2("nTotEvents.dat");
157
158   Int_t    ev;
159   Int_t    nTotEv=0,nD0rec=0,nD0rec1ev=0;
160   Double_t dca;
161   Double_t v2[3],mom[6],d0[2];
162   Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
163   Int_t    iTrkP,iTrkN,trkEntries;
164   Int_t    nTrksP=0,nTrksN=0;
165   Int_t    trkNum[2];
166   Int_t    okD0=0,okD0bar=0;
167   Char_t   trkTreeName[100],vtxName[100];
168   AliITStrackV2 *postrack = 0;
169   AliITStrackV2 *negtrack = 0;
170
171   // create the AliITSVertexerTracks object
172   // (it will be used only if fVertexOnTheFly=kTrue)
173   AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
174   vertexer1->SetMinTracks(2);
175   vertexer1->SetDebug(0);
176   Int_t  skipped[2];
177   Bool_t goodVtx1;
178   
179
180   // define the cuts for vertexing
181   Double_t vtxcuts[]={50., // max. allowed chi2
182                       0.0, // min. allowed negative daughter's impact param 
183                       0.0, // min. allowed positive daughter's impact param 
184                       1.0, // max. allowed DCA between the daughter tracks
185                      -1.0, // min. allowed cosine of pointing angle
186                       0.0, // min. radius of the fiducial volume
187                       2.9};// max. radius of the fiducial volume
188   
189   // create the AliV0vertexer object
190   AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
191
192   // create tree for reconstructed D0s
193   AliD0toKpi *ioD0toKpi=0;
194   TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
195   treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
196
197   // open file with tracks
198   TFile *trkFile = TFile::Open(trkName.Data());
199
200   // loop on events in file
201   for(ev=evFirst; ev<=evLast; ev++) {
202     printf(" --- Looking for D0->Kpi in event  %d ---\n",ev);
203
204     // retrieve primary vertex from file
205     sprintf(vtxName,"VertexTracks_%d",ev);
206     AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
207     vertex1stored->GetXYZ(fV1);
208     delete vertex1stored;
209
210     // retrieve tracks from file
211     sprintf(trkTreeName,"TreeT_ITS_%d",ev);
212     TTree *trkTree=(TTree*)trkFile->Get(trkTreeName);
213     if(!trkTree) { 
214       printf("AliD0toKpiAnalysis::FindCandidates():\n  tracks tree not found for evet %d\n",ev); 
215       continue; 
216     }
217     trkEntries = (Int_t)trkTree->GetEntries();
218     printf(" Number of tracks: %d\n",trkEntries);
219
220     // count the total number of events
221     nTotEv++;
222
223     // call function which applies sigle-track selection and
224     // separetes positives and negatives
225     TObjArray trksP(trkEntries/2);
226     Int_t *trkEntryP = new Int_t[trkEntries];
227     TObjArray trksN(trkEntries/2);
228     Int_t *trkEntryN = new Int_t[trkEntries];
229     SelectTracks(*trkTree,trksP,trkEntryP,nTrksP,trksN,trkEntryN,nTrksN);
230       
231     nD0rec1ev = 0;
232
233     // loop on positive tracks
234     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
235       if(iTrkP%10==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
236           
237       // get track from track array
238       postrack = (AliITStrackV2*)trksP.At(iTrkP);
239       trkNum[0] = trkEntryP[iTrkP];      
240
241       // loop on negative tracks 
242       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
243         // get track from tracks array
244         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
245         trkNum[1] = trkEntryN[iTrkN];      
246
247         AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
248
249         //
250         // ----------- DCA MINIMIZATION ------------------
251         //
252         // find the DCA and propagate the tracks to the DCA 
253         dca = vertexer2->PropagateToDCA(pnt,ppt);
254
255         // define the AliV0vertex object
256         AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
257           
258         // get position of the secondary vertex
259         vertex2->GetXYZ(v2[0],v2[1],v2[2]);
260         
261         delete vertex2;
262   
263         // momenta of the tracks at the vertex
264         ptP = 1./TMath::Abs(ppt->Get1Pt());
265         alphaP = ppt->GetAlpha();
266         phiP = alphaP+TMath::ASin(ppt->GetSnp());
267         mom[0] = ptP*TMath::Cos(phiP); 
268         mom[1] = ptP*TMath::Sin(phiP);
269         mom[2] = ptP*ppt->GetTgl();
270           
271         ptN = 1./TMath::Abs(pnt->Get1Pt());
272         alphaN = pnt->GetAlpha();
273         phiN = alphaN+TMath::ASin(pnt->GetSnp());
274         mom[3] = ptN*TMath::Cos(phiN); 
275         mom[4] = ptN*TMath::Sin(phiN);
276         mom[5] = ptN*pnt->GetTgl();
277           
278         goodVtx1 = kTRUE;
279         // no vertexing if DeltaMass > fMassCut 
280         if(fVertexOnTheFly) {
281           goodVtx1 = kFALSE;
282           if(SelectInvMass(mom)) {
283             // primary vertex from *other* tracks in event
284             vertexer1->SetVtxStart(fV1[0],fV1[1]);
285             skipped[0] = trkEntryP[iTrkP];
286             skipped[1] = trkEntryN[iTrkN];
287             vertexer1->SetSkipTracks(2,skipped);
288             AliESDVertex *vertex1onfly = 
289               (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree); 
290             if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
291             vertex1onfly->GetXYZ(fV1);
292             //vertex1onfly->PrintStatus();
293             delete vertex1onfly;        
294           }
295         }         
296
297         // impact parameters of the tracks w.r.t. the primary vertex
298         d0[0] =  10000.*ppt->GetD(fV1[0],fV1[1]);
299         d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
300
301         // create the object AliD0toKpi
302         AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
303
304         // select D0s
305         if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
306               
307           // fill the tree
308           ioD0toKpi=&theD0;
309           treeD0->Fill();
310
311           nD0rec++; nD0rec1ev++;
312           ioD0toKpi=0;  
313         }
314         
315         negtrack = 0;
316       } // loop on negative tracks
317       postrack = 0;
318     }   // loop on positive tracks
319
320     trksP.Delete();
321     trksN.Delete();
322     delete [] trkEntryP;
323     delete [] trkEntryN;
324     delete trkTree;
325
326     printf(" Number of D0 candidates: %d\n",nD0rec1ev);
327   }    // loop on events in file
328
329
330   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
331   printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
332
333   delete vertexer1;
334   delete vertexer2;
335
336   trkFile->Close();
337
338   // create a copy of this class to be written to output file
339   AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
340   copy->PrintStatus();
341
342   // add PDG codes to decay tracks in found candidates (in simulation mode)
343   // and store tree in the output file
344   if(!fSim) {
345     TFile *outroot = new TFile(outName1.Data(),"recreate");
346     treeD0->Write();
347     copy->Write();
348     outroot->Close();
349     delete outroot;
350   } else {
351     printf(" Now adding information from simulation (PDG codes) ...\n");
352     TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
353     MakeTracksRefFile(evFirst,evLast);
354     SimulationInfo(treeD0,treeD0sim);
355     delete treeD0;
356     TFile *outroot = new TFile(outName1.Data(),"recreate");
357     treeD0sim->Write();
358     copy->Write();
359     outroot->Close();
360     delete outroot;
361   }
362
363   // write to a file the total number of events
364   FILE *outdat = fopen(outName2.Data(),"w");
365   fprintf(outdat,"%d\n",nTotEv);
366   fclose(outdat);
367
368   return;
369 }
370 //----------------------------------------------------------------------------
371 void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
372                                            const Char_t *outName) {
373   // Find D0 candidates and calculate parameters
374
375   if(fBz<-9000.) {
376     printf("AliD0toKpiAnalysis::FindCandidatesESD():  Set B!\n");
377     return;
378   }
379   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
380
381   TString esdName("AliESDs.root");
382   if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
383     printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n"); 
384     return;
385   }
386
387   TString outName1=outName;
388   TString outName2("nTotEvents.dat");
389
390   Int_t    nTotEv=0,nD0rec=0,nD0rec1ev=0;
391   Double_t dca;
392   Double_t v2[3],mom[6],d0[2];
393   Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
394   Int_t    iTrkP,iTrkN,trkEntries;
395   Int_t    nTrksP=0,nTrksN=0;
396   Int_t    trkNum[2];
397   Double_t tofmass[2];
398   Int_t    okD0=0,okD0bar=0;
399   AliITStrackV2 *postrack = 0;
400   AliITStrackV2 *negtrack = 0;
401
402   // create the AliITSVertexerTracks object
403   // (it will be used only if fVertexOnTheFly=kTrue)
404   AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
405   vertexer1->SetMinTracks(2);
406   vertexer1->SetDebug(0);
407   Int_t  skipped[2];
408   Bool_t goodVtx1;
409   
410
411   // define the cuts for vertexing
412   Double_t vtxcuts[]={50., // max. allowed chi2
413                       0.0, // min. allowed negative daughter's impact param 
414                       0.0, // min. allowed positive daughter's impact param 
415                       1.0, // max. allowed DCA between the daughter tracks
416                      -1.0, // min. allowed cosine of pointing angle
417                       0.0, // min. radius of the fiducial volume
418                       2.9};// max. radius of the fiducial volume
419   
420   // create the AliV0vertexer object
421   AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
422
423   // create tree for reconstructed D0s
424   AliD0toKpi *ioD0toKpi=0;
425   TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
426   treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
427
428   // open file with tracks
429   TFile *esdFile = TFile::Open(esdName.Data());
430   AliESD* event = new AliESD;
431   TTree* tree = (TTree*) esdFile->Get("esdTree");
432   if (!tree) {
433     Error("FindCandidatesESD", "no ESD tree found");
434     return;
435   }
436   tree->SetBranchAddress("ESD", &event);
437
438   // loop on events in file
439   for (Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
440     if (iEvent > evLast) break;
441     tree->GetEvent(iEvent);
442     Int_t ev = (Int_t)event->GetEventNumber();
443     printf("--- Finding D0 -> Kpi in event %d\n",ev);
444
445     // retrieve primary vertex from file
446     //sprintf(vtxName,"Vertex_%d",ev);
447     //AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
448     //vertex1stored->GetXYZ(fV1);
449     //delete vertex1stored;
450     event->GetVertex()->GetXYZ(fV1);
451
452     trkEntries = event->GetNumberOfTracks();
453     printf(" Number of tracks: %d\n",trkEntries);
454
455     // count the total number of events
456     nTotEv++;
457
458     // call function which applies sigle-track selection and
459     // separetes positives and negatives
460     TObjArray trksP(trkEntries/2);
461     Int_t    *trkEntryP   = new Int_t[trkEntries];
462     TObjArray trksN(trkEntries/2);
463     Int_t    *trkEntryN   = new Int_t[trkEntries];
464     TTree *trkTree = new TTree();
465     if(fVertexOnTheFly) {
466       SelectTracksESDvtx(*event,trkTree,trksP,trkEntryP,nTrksP,
467                                         trksN,trkEntryN,nTrksN);
468     } else {
469       SelectTracksESD(*event,trksP,trkEntryP,nTrksP,
470                              trksN,trkEntryN,nTrksN);
471     }      
472
473     if(fDebug) printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
474
475     nD0rec1ev = 0;
476
477     // loop on positive tracks
478     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
479       if((iTrkP%10==0) || fDebug) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
480           
481       // get track from track array
482       postrack = (AliITStrackV2*)trksP.UncheckedAt(iTrkP);
483       trkNum[0] = trkEntryP[iTrkP];      
484
485       // loop on negative tracks 
486       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
487
488         // get track from tracks array
489         negtrack = (AliITStrackV2*)trksN.UncheckedAt(iTrkN);
490         trkNum[1] = trkEntryN[iTrkN];      
491
492         AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
493
494         //
495         // ----------- DCA MINIMIZATION ------------------
496         //
497         // find the DCA and propagate the tracks to the DCA 
498         dca = vertexer2->PropagateToDCA(pnt,ppt);
499
500         // define the AliV0vertex object
501         AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
502           
503         // get position of the secondary vertex
504         vertex2->GetXYZ(v2[0],v2[1],v2[2]);
505         
506         delete vertex2;
507   
508         // momenta of the tracks at the vertex
509         ptP = 1./TMath::Abs(ppt->Get1Pt());
510         alphaP = ppt->GetAlpha();
511         phiP = alphaP+TMath::ASin(ppt->GetSnp());
512         mom[0] = ptP*TMath::Cos(phiP); 
513         mom[1] = ptP*TMath::Sin(phiP);
514         mom[2] = ptP*ppt->GetTgl();
515           
516         ptN = 1./TMath::Abs(pnt->Get1Pt());
517         alphaN = pnt->GetAlpha();
518         phiN = alphaN+TMath::ASin(pnt->GetSnp());
519         mom[3] = ptN*TMath::Cos(phiN); 
520         mom[4] = ptN*TMath::Sin(phiN);
521         mom[5] = ptN*pnt->GetTgl();
522           
523         goodVtx1 = kTRUE;
524         // no vertexing if DeltaMass > fMassCut 
525         if(fVertexOnTheFly) {
526           goodVtx1 = kFALSE;
527           if(SelectInvMass(mom)) {
528             // primary vertex from *other* tracks in event
529             vertexer1->SetVtxStart(fV1[0],fV1[1]);
530             skipped[0] = trkEntryP[iTrkP];
531             skipped[1] = trkEntryN[iTrkN];
532             vertexer1->SetSkipTracks(2,skipped);
533             AliESDVertex *vertex1onfly = 
534               (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree); 
535             if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
536             vertex1onfly->GetXYZ(fV1);
537             //vertex1onfly->PrintStatus();
538             delete vertex1onfly;        
539           }
540         }         
541
542         // impact parameters of the tracks w.r.t. the primary vertex
543         d0[0] =  10000.*ppt->GetD(fV1[0],fV1[1]);
544         d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
545
546         // create the object AliD0toKpi
547         AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
548
549         // select D0s
550         if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
551               
552           // get PID info from ESD
553           AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
554           Double_t esdpid0[5];
555           t0->GetESDpid(esdpid0);
556           if(t0->GetStatus()&AliESDtrack::kTOFpid) {
557             tofmass[0] = CalculateTOFmass(t0->GetP(),
558                                           t0->GetIntegratedLength(),
559                                           t0->GetTOFsignal());
560           } else {
561             tofmass[0] = -1000.;
562           }
563           AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
564           Double_t esdpid1[5];
565           t1->GetESDpid(esdpid1);
566           if(t1->GetStatus()&AliESDtrack::kTOFpid) {
567             tofmass[1] = CalculateTOFmass(t1->GetP(),
568                                           t1->GetIntegratedLength(),
569                                           t1->GetTOFsignal());
570           } else {
571             tofmass[1] = -1000.;
572           }
573
574           theD0.SetPIDresponse(esdpid0,esdpid1);
575           theD0.SetTOFmasses(tofmass);
576
577           // fill the tree
578           ioD0toKpi=&theD0;
579           treeD0->Fill();
580
581           nD0rec++; nD0rec1ev++;
582           ioD0toKpi=0;  
583         }
584         
585         negtrack = 0;
586       } // loop on negative tracks
587       postrack = 0;
588     }   // loop on positive tracks
589
590     trksP.Delete();
591     trksN.Delete();
592     delete [] trkEntryP;
593     delete [] trkEntryN;
594     delete trkTree;
595
596
597     printf(" Number of D0 candidates: %d\n",nD0rec1ev);
598   }    // loop on events in file
599
600
601   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
602   printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
603
604   delete vertexer1;
605   delete vertexer2;
606
607   esdFile->Close();
608
609   // create a copy of this class to be written to output file
610   AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
611
612   // add PDG codes to decay tracks in found candidates (in simulation mode)
613   // and store tree in the output file
614   if(!fSim) {
615     TFile *outroot = new TFile(outName1.Data(),"recreate");
616     treeD0->Write();
617     copy->Write();
618     outroot->Close();
619     delete outroot;
620   } else {
621     printf(" Now adding information from simulation (PDG codes) ...\n");
622     TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
623     MakeTracksRefFileESD();
624     SimulationInfo(treeD0,treeD0sim);
625     delete treeD0;
626     TFile *outroot = new TFile(outName1.Data(),"recreate");
627     treeD0sim->Write();
628     copy->Write();
629     outroot->Close();
630     delete outroot;
631   }
632
633   // write to a file the total number of events
634   FILE *outdat = fopen(outName2.Data(),"w");
635   fprintf(outdat,"%d\n",nTotEv);
636   fclose(outdat);
637
638   return;
639 }
640 //-----------------------------------------------------------------------------
641 void AliD0toKpiAnalysis::PrintStatus() const {
642   // Print parameters being used
643
644   printf(" fBz      = %f T\n",fBz);
645   printf("Preselections:\n");
646   printf("    fPtCut   = %f GeV\n",fPtCut);
647   printf("    fd0Cut   = %f micron\n",fd0Cut);
648   printf("    fMassCut = %f GeV\n",fMassCut);
649   if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
650   if(fSim) { 
651     printf("Simulation mode\n");
652     if(fOnlySignal) printf("  Only signal goes to file\n");
653   }
654   printf("Cuts on candidates:\n");
655   printf("    |M-MD0| [GeV]    < %f\n",fD0Cuts[0]);
656   printf("    dca    [micron]  < %f\n",fD0Cuts[1]);
657   printf("    cosThetaStar     < %f\n",fD0Cuts[2]);
658   printf("    pTK     [GeV]    > %f\n",fD0Cuts[3]);
659   printf("    pTpi    [GeV]    > %f\n",fD0Cuts[4]);
660   printf("    |d0K|  [micron]  < %f\n",fD0Cuts[5]);
661   printf("    |d0pi| [micron]  < %f\n",fD0Cuts[6]);
662   printf("    d0d0  [micron^2] < %f\n",fD0Cuts[7]);
663   printf("    cosThetaPoint    > %f\n",fD0Cuts[8]);
664
665   return;
666 }
667 //-----------------------------------------------------------------------------
668 Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
669   // Apply preselection in the invariant mass of the pair
670
671   Double_t mD0 = 1.8645;
672   Double_t mPi = 0.13957;
673   Double_t mKa = 0.49368;
674
675   Double_t energy[2];
676   Double_t mom2[2],momTot2;
677
678   mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
679   mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
680
681   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
682             (p[1]+p[4])*(p[1]+p[4])+
683             (p[2]+p[5])*(p[2]+p[5]);
684
685   // D0 -> K- pi+
686   energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
687   energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
688
689   Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
690     
691   // D0bar -> K+ pi-
692   energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
693   energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
694
695   Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
696
697   if(TMath::Abs(minvD0-mD0)    < fMassCut) return kTRUE;
698   if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
699   return kFALSE;
700 }
701 //-----------------------------------------------------------------------------
702 void AliD0toKpiAnalysis::SelectTracks(TTree &trkTree,
703                   TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
704                   TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
705   // Create two TObjArrays with positive and negative tracks and 
706   // apply single-track preselection
707
708   nTrksP=0,nTrksN=0;
709  
710   Int_t entr = (Int_t)trkTree.GetEntries();
711
712   // trasfer tracks from tree to arrays
713   for(Int_t i=0; i<entr; i++) {
714
715     AliITStrackV2 *track = new AliITStrackV2; 
716     trkTree.SetBranchAddress("tracks",&track);
717
718     trkTree.GetEvent(i);
719
720     // single track selection
721     if(!SingleTrkCuts(*track)) { delete track; continue; }
722
723     if(track->Get1Pt()>0.) { // negative track
724       trksN.AddLast(track);
725       trkEntryN[nTrksN] = i;
726       nTrksN++;
727     } else {                 // positive track
728       trksP.AddLast(track);
729       trkEntryP[nTrksP] = i;
730       nTrksP++;
731     }
732
733   }
734
735   return;
736 }
737 //-----------------------------------------------------------------------------
738 void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
739         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
740         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
741   // Create two TObjArrays with positive and negative tracks and 
742   // apply single-track preselection
743
744   nTrksP=0,nTrksN=0;
745
746   Int_t entr = event.GetNumberOfTracks();
747  
748   // transfer ITS tracks from ESD to arrays and to a tree
749   for(Int_t i=0; i<entr; i++) {
750     //if(fDebug) printf(" SelectTracksESD: %d/%d\n",i,entr);
751
752     AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
753     UInt_t status = esdtrack->GetStatus();
754
755     if(!(status&AliESDtrack::kITSrefit)) continue;
756
757     AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack); 
758
759     // single track selection
760     if(!SingleTrkCuts(*itstrack)) 
761       { delete itstrack; continue; }
762
763     if(itstrack->Get1Pt()>0.) { // negative track
764       trksN.AddLast(itstrack);
765       trkEntryN[nTrksN] = i;
766       nTrksN++;
767     } else {                 // positive track
768       trksP.AddLast(itstrack);
769       trkEntryP[nTrksP] = i;
770       nTrksP++;
771     }
772
773   } // loop on ESD tracks
774
775   return;
776 }
777 //-----------------------------------------------------------------------------
778 void AliD0toKpiAnalysis::SelectTracksESDvtx(AliESD &event,TTree *trkTree,
779         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
780         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
781   // Create two TObjArrays with positive and negative tracks and 
782   // apply single-track preselection
783
784   nTrksP=0,nTrksN=0;
785
786   Int_t entr = event.GetNumberOfTracks();
787  
788   AliITStrackV2 *itstrackfortree = 0;
789   trkTree->Branch("tracks","AliITStrackV2",&itstrackfortree,entr,0);
790
791
792   // transfer ITS tracks from ESD to arrays and to a tree
793   for(Int_t i=0; i<entr; i++) {
794
795     AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
796     UInt_t status = esdtrack->GetStatus();
797
798     if(!(status&AliESDtrack::kITSrefit)) continue;
799
800     AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack);
801
802     // store track in the tree to be used for primary vertex finding
803     itstrackfortree = new AliITStrackV2(*esdtrack);
804     trkTree->Fill();
805     itstrackfortree = 0;
806
807     // single track selection
808     if(!SingleTrkCuts(*itstrack)) 
809       { delete itstrack; continue; }
810
811     if(itstrack->Get1Pt()>0.) { // negative track
812       trksN.AddLast(itstrack);
813       trkEntryN[nTrksN] = i;
814       nTrksN++;
815     } else {                 // positive track
816       trksP.AddLast(itstrack);
817       trkEntryP[nTrksP] = i;
818       nTrksP++;
819     }
820
821   } // loop on esd tracks
822
823   delete itstrackfortree;
824
825   return;
826 }
827 //-----------------------------------------------------------------------------
828 void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
829                                    Double_t cut2,Double_t cut3,Double_t cut4,
830                                    Double_t cut5,Double_t cut6,
831                                    Double_t cut7,Double_t cut8) {
832   // Set the cuts for D0 selection
833   fD0Cuts[0] = cut0;
834   fD0Cuts[1] = cut1;
835   fD0Cuts[2] = cut2;
836   fD0Cuts[3] = cut3;
837   fD0Cuts[4] = cut4;
838   fD0Cuts[5] = cut5;
839   fD0Cuts[6] = cut6;
840   fD0Cuts[7] = cut7;
841   fD0Cuts[8] = cut8;
842
843   return;
844 }
845 //-----------------------------------------------------------------------------
846 void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
847   // Set the cuts for D0 selection
848
849   for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
850
851   return;
852 }
853 //-----------------------------------------------------------------------------
854 Bool_t AliD0toKpiAnalysis::SingleTrkCuts(const AliITStrackV2& trk) const {
855   // Check if track passes some kinematical cuts  
856
857   if(TMath::Abs(1./trk.Get1Pt()) < fPtCut) 
858     return kFALSE;
859   if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1])) < fd0Cut) 
860     return kFALSE;
861
862   return kTRUE;
863 }
864 //----------------------------------------------------------------------------
865 void AliD0toKpiAnalysis::MakeTracksRefFile(Int_t evFirst,Int_t evLast) 
866   const {
867   // Create a file with simulation info for the reconstructed tracks
868   
869   TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
870   TFile *trk = TFile::Open("AliITStracksV2.root");
871   AliRunLoader *rl = AliRunLoader::Open("galice.root");
872
873   // load kinematics
874   rl->LoadKinematics();
875   
876   Int_t      label;
877   TParticle *part;  
878   TParticle *mumpart;
879   REFTRACK   reftrk;
880   
881   for(Int_t ev=evFirst; ev<=evLast; ev++){
882     rl->GetEvent(ev);  
883     AliStack *stack = rl->Stack();
884
885     trk->cd();
886
887     // Tree with tracks
888     char tname[100];
889     sprintf(tname,"TreeT_ITS_%d",ev);
890
891     TTree *tracktree=(TTree*)trk->Get(tname);
892     if(!tracktree) continue;
893     AliITStrackV2 *track = new AliITStrackV2; 
894     tracktree->SetBranchAddress("tracks",&track);
895     Int_t nentr=(Int_t)tracktree->GetEntries();
896
897     // Tree for true track parameters
898     char ttname[100];
899     sprintf(ttname,"Tree_Ref_%d",ev);
900     TTree *reftree = new TTree(ttname,"Tree with true track params");
901     reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
902
903     for(Int_t i=0; i<nentr; i++) {
904       tracktree->GetEvent(i);
905       label = TMath::Abs(track->GetLabel());
906
907       part = (TParticle*)stack->Particle(label);
908       reftrk.lab = label;
909       reftrk.pdg = part->GetPdgCode();
910       reftrk.mumlab = part->GetFirstMother();
911       if(part->GetFirstMother()>=0) {
912         mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
913         reftrk.mumpdg = mumpart->GetPdgCode();
914       } else {
915         reftrk.mumpdg=-1;
916       }
917       reftrk.Vx = part->Vx();
918       reftrk.Vy = part->Vy();
919       reftrk.Vz = part->Vz();
920       reftrk.Px = part->Px();
921       reftrk.Py = part->Py();
922       reftrk.Pz = part->Pz();
923       
924       reftree->Fill();
925     } // loop on tracks   
926
927     out->cd();
928     reftree->Write();
929
930     delete track;
931     delete reftree;
932     delete tracktree;
933     delete stack;
934   } // loop on events
935
936   trk->Close();
937   out->Close();
938
939   return;
940 }
941 //----------------------------------------------------------------------------
942 void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
943   // Create a file with simulation info for the reconstructed tracks
944   
945   TFile *outFile = TFile::Open("ITStracksRefFile.root","recreate");
946   TFile *esdFile = TFile::Open("AliESDs.root");
947   AliRunLoader *rl = AliRunLoader::Open("galice.root");
948
949   // load kinematics
950   rl->LoadKinematics();
951   
952   Int_t      label;
953   TParticle *part;  
954   TParticle *mumpart;
955   REFTRACK   reftrk;
956   
957   TKey *key=0;
958   TIter next(esdFile->GetListOfKeys());
959   // loop on events in file
960   while ((key=(TKey*)next())!=0) {
961     AliESD *event=(AliESD*)key->ReadObj();
962     Int_t ev = (Int_t)event->GetEventNumber();
963
964     rl->GetEvent(ev);  
965     AliStack *stack = rl->Stack();
966
967     Int_t nentr=(Int_t)event->GetNumberOfTracks();
968
969     // Tree for true track parameters
970     char ttname[100];
971     sprintf(ttname,"Tree_Ref_%d",ev);
972     TTree *reftree = new TTree(ttname,"Tree with true track params");
973     reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
974
975     for(Int_t i=0; i<nentr; i++) {
976       AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
977       label = TMath::Abs(esdtrack->GetLabel());
978
979       part = (TParticle*)stack->Particle(label);
980       reftrk.lab = label;
981       reftrk.pdg = part->GetPdgCode();
982       reftrk.mumlab = part->GetFirstMother();
983       if(part->GetFirstMother()>=0) {
984         mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
985         reftrk.mumpdg = mumpart->GetPdgCode();
986       } else {
987         reftrk.mumpdg=-1;
988       }
989       reftrk.Vx = part->Vx();
990       reftrk.Vy = part->Vy();
991       reftrk.Vz = part->Vz();
992       reftrk.Px = part->Px();
993       reftrk.Py = part->Py();
994       reftrk.Pz = part->Pz();
995       
996       reftree->Fill();
997
998     } // loop on tracks   
999
1000     outFile->cd();
1001     reftree->Write();
1002
1003     delete reftree;
1004     delete event;
1005     delete stack;
1006   } // loop on events
1007
1008   esdFile->Close();
1009   outFile->Close();
1010
1011   return;
1012 }
1013 //-----------------------------------------------------------------------------
1014 void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
1015   // add pdg codes to candidate decay tracks (for sim)
1016
1017   TString refFileName("ITStracksRefFile.root");
1018   if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) { 
1019     printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n"); 
1020     return;
1021   }
1022   TFile *refFile = TFile::Open(refFileName.Data());
1023
1024   Char_t refTreeName[100];
1025   Int_t  event;
1026   Int_t  pdg[2],mumpdg[2],mumlab[2];
1027   REFTRACK reftrk;
1028
1029   // read-in reference tree for event 0 (the only event for Pb-Pb)
1030   sprintf(refTreeName,"Tree_Ref_%d",0);
1031   TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
1032   refTree0->SetBranchAddress("rectracks",&reftrk);
1033
1034   AliD0toKpi *theD0 = 0; 
1035   treeD0in->SetBranchAddress("D0toKpi",&theD0);
1036   treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
1037
1038   Int_t entries = (Int_t)treeD0in->GetEntries();
1039
1040   for(Int_t i=0; i<entries; i++) {
1041     if(i%100==0) printf("  done %d candidates of %d\n",i,entries);    
1042
1043     treeD0in->GetEvent(i);
1044     event = theD0->EventNo();
1045
1046     if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
1047       refTree0->GetEvent(theD0->GetTrkNum(0));
1048       pdg[0]    = reftrk.pdg;
1049       mumpdg[0] = reftrk.mumpdg;
1050       mumlab[0] = reftrk.mumlab;
1051       refTree0->GetEvent(theD0->GetTrkNum(1));
1052       pdg[1]    = reftrk.pdg;
1053       mumpdg[1] = reftrk.mumpdg;
1054       mumlab[1] = reftrk.mumlab;
1055     } else {
1056       sprintf(refTreeName,"Tree_Ref_%d",event);
1057       TTree *refTree = (TTree*)refFile->Get(refTreeName);
1058       refTree->SetBranchAddress("rectracks",&reftrk);
1059       refTree->GetEvent(theD0->GetTrkNum(0));
1060       pdg[0]    = reftrk.pdg;
1061       mumpdg[0] = reftrk.mumpdg;
1062       mumlab[0] = reftrk.mumlab;
1063       refTree->GetEvent(theD0->GetTrkNum(1));
1064       pdg[1]    = reftrk.pdg;
1065       mumpdg[1] = reftrk.mumpdg;
1066       mumlab[1] = reftrk.mumlab;
1067       delete refTree;
1068     }
1069     
1070     theD0->SetPdgCodes(pdg);
1071     theD0->SetMumPdgCodes(mumpdg);
1072     
1073     if(TMath::Abs(mumpdg[0])==421 && TMath::Abs(mumpdg[1])==421 
1074        && mumlab[0]==mumlab[1]) theD0->SetSignal();
1075     
1076     if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
1077
1078   }
1079
1080   delete refTree0;
1081
1082   refFile->Close();
1083
1084   return;
1085 }
1086
1087
1088
1089
1090
1091