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