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