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