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