move D0toKpi to PWG3 (Andrea)
[u/mrichter/AliRoot.git] / PWG3 / 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@lnl.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 "AliMC.h"
31 #include "AliRun.h"
32 #include "AliRunLoader.h"
33 #include "AliVertexerTracks.h"
34 #include "AliESDVertex.h"
35 #include "AliESDv0.h"
36 #include "AliD0toKpi.h"
37 #include "AliD0toKpiAnalysis.h"
38 #include "AliLog.h"
39
40 typedef struct {
41   Int_t lab;
42   Int_t pdg;
43   Int_t mumlab;
44   Int_t mumpdg;
45   Int_t mumprongs;
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   SetPtCut();
57   Setd0Cut();
58   SetMassCut();
59   SetD0Cuts();
60   SetVertex1();
61   SetPID();
62   fVertexOnTheFly = kFALSE;
63   fSim = kFALSE;
64   fOnlySignal = kFALSE;
65 }
66 //----------------------------------------------------------------------------
67 AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
68 //----------------------------------------------------------------------------
69 void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
70   // select candidates that pass fD0Cuts and write them to a new file
71
72   TFile *inFile = TFile::Open(inName);
73
74   TTree *treeD0in=(TTree*)inFile->Get("TreeD0");
75   AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis");
76   printf("+++\n+++  I N P U T   S T A T U S:\n+++\n");
77   inAnalysis->PrintStatus();
78
79
80   AliD0toKpi *d = 0; 
81   treeD0in->SetBranchAddress("D0toKpi",&d);
82   Int_t entries = (Int_t)treeD0in->GetEntries();
83
84   printf("+++\n+++ Number of D0 in input tree:  %d\n+++\n",entries);
85
86   TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
87   treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
88
89
90   Int_t okD0=0,okD0bar=0;
91   Int_t nSel = 0;
92
93   for(Int_t i=0; i<entries; i++) {
94     // get event from tree
95     treeD0in->GetEvent(i);
96
97     if(fSim && fOnlySignal && !d->IsSignal()) continue; 
98
99     // check if candidate passes selection (as D0 or D0bar)
100     if(d->Select(fD0Cuts,okD0,okD0bar)) {
101       nSel++;
102       treeD0out->Fill();
103     }
104
105   }
106
107   AliD0toKpiAnalysis *outAnalysis = (AliD0toKpiAnalysis*)inAnalysis->Clone("D0toKpiAnalysis");
108   outAnalysis->SetD0Cuts(fD0Cuts);
109   printf("------------------------------------------\n");
110   printf("+++\n+++  O U T P U T   S T A T U S:\n+++\n");
111   outAnalysis->PrintStatus();
112
113   printf("+++\n+++ Number of D0 in output tree:  %d\n+++\n",nSel);
114
115   TFile* outFile = new TFile(outName,"recreate");
116   treeD0out->Write();
117   outAnalysis->Write();
118   outFile->Close();
119
120   return;
121 }
122 //----------------------------------------------------------------------------
123 Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
124                                               Double_t time) const {
125   // calculated the mass from momentum, track length from vertex to TOF
126   // and time measured by the TOF
127   if(length==0.) return -1000.;
128   Double_t a = time*time/length/length;
129   if(a > 1.) {
130     a = TMath::Sqrt(a-1.);
131   } else {
132     a = -TMath::Sqrt(1.-a);
133   }
134
135   return mom*a;
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
143   TString esdName="AliESDs.root";
144   if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
145     printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n"); 
146     return;
147   }
148
149   TString outName1=outName;
150   TString outName2="nTotEvents.dat";
151
152   Int_t    nTotEv=0,nD0rec=0,nD0rec1ev=0;
153   Double_t dca;
154   Double_t v2[3],mom[6],d0[2];
155   Int_t    iTrkP,iTrkN,trkEntries;
156   Int_t    nTrksP=0,nTrksN=0;
157   Int_t    trkNum[2];
158   Double_t tofmass[2];
159   Int_t    okD0=0,okD0bar=0;
160   AliESDtrack *postrack = 0;
161   AliESDtrack *negtrack = 0;
162
163   // create the AliVertexerTracks object
164   // (it will be used only if fVertexOnTheFly=kTrue)
165   AliVertexerTracks *vertexer1 = new AliVertexerTracks;
166   if(fVertexOnTheFly) {
167     // open the mean vertex
168     TFile *invtx = new TFile("AliESDVertexMean.root");
169     AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
170     invtx->Close();
171     vertexer1->SetVtxStart(initVertex);
172     delete invtx;
173   }
174   Int_t  skipped[2];
175   Bool_t goodVtx1;
176   
177   // create tree for reconstructed D0s
178   AliD0toKpi *ioD0toKpi=0;
179   TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
180   treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
181
182   // open file with tracks
183   TFile *esdFile = TFile::Open(esdName.Data());
184   AliESD* event = new AliESD;
185   TTree* tree = (TTree*) esdFile->Get("esdTree");
186   if(!tree) {
187     Error("FindCandidatesESD", "no ESD tree found");
188     return;
189   }
190   tree->SetBranchAddress("ESD",&event);
191
192   // loop on events in file
193   for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
194     if(iEvent > evLast) break;
195     tree->GetEvent(iEvent);
196     Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
197     printf("--- Finding D0 -> Kpi in event %d\n",ev);
198     // count the total number of events
199     nTotEv++;
200
201     trkEntries = (Int_t)event->GetNumberOfTracks();
202     printf(" Number of tracks: %d\n",trkEntries);
203     if(trkEntries<2) continue;
204
205     // retrieve primary vertex from file
206     if(!event->GetPrimaryVertex()) { 
207       printf(" No vertex\n");
208       continue;
209     }
210     event->GetPrimaryVertex()->GetXYZ(fV1);
211
212     // call function which applies sigle-track selection and
213     // separetes positives and negatives
214     TObjArray trksP(trkEntries/2);
215     Int_t    *trkEntryP   = new Int_t[trkEntries];
216     TObjArray trksN(trkEntries/2);
217     Int_t    *trkEntryN   = new Int_t[trkEntries];
218     TTree *trkTree = new TTree();
219     SelectTracks(event,trksP,trkEntryP,nTrksP,
220                        trksN,trkEntryN,nTrksN);
221
222     printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
223
224
225     nD0rec1ev = 0;
226
227     // LOOP ON  POSITIVE  TRACKS
228     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
229       if(iTrkP%1==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
230           
231       // get track from track array
232       postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
233       trkNum[0] = trkEntryP[iTrkP];      
234
235       // LOOP ON  NEGATIVE  TRACKS
236       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
237
238         // get track from tracks array
239         negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
240         trkNum[1] = trkEntryN[iTrkN];      
241
242         //
243         // ----------- DCA MINIMIZATION ------------------
244         //
245         // find the DCA and propagate the tracks to the DCA
246         Double_t b=event->GetMagneticField(); 
247         AliESDtrack nt(*negtrack), pt(*postrack);
248         dca = nt.PropagateToDCA(&pt,b);
249
250         // define the AliESDv0 object
251         AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
252           
253         // get position of the secondary vertex
254         vertex2.GetXYZ(v2[0],v2[1],v2[2]);
255         vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
256         vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
257         // impact parameters of the tracks w.r.t. the primary vertex
258         d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
259         d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
260         goodVtx1 = kTRUE;
261         
262         // no vertexing if DeltaMass > fMassCut 
263         if(fVertexOnTheFly) {
264           goodVtx1 = kFALSE;
265           if(SelectInvMass(mom)) {
266             // primary vertex from *other* tracks in the event
267             skipped[0] = trkEntryP[iTrkP];
268             skipped[1] = trkEntryN[iTrkN];
269             vertexer1->SetSkipTracks(2,skipped);
270             AliESDVertex *vertex1onfly = 
271               (AliESDVertex*)vertexer1->FindPrimaryVertex(event); 
272             if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
273             vertex1onfly->GetXYZ(fV1);
274             //vertex1onfly->PrintStatus();
275             delete vertex1onfly;        
276           }
277         }
278         
279
280         // create the object AliD0toKpi
281         AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
282         // select D0s
283         if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
284           // get PID info from ESD
285           AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
286           Double_t esdpid0[5];
287           t0->GetESDpid(esdpid0);
288           if(t0->GetStatus()&AliESDtrack::kTOFpid) {
289             tofmass[0] = CalculateTOFmass(t0->GetP(),
290                                           t0->GetIntegratedLength(),
291                                           t0->GetTOFsignal());
292           } else {
293             tofmass[0] = -1000.;
294           }
295           AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
296           Double_t esdpid1[5];
297           t1->GetESDpid(esdpid1);
298           if(t1->GetStatus()&AliESDtrack::kTOFpid) {
299             tofmass[1] = CalculateTOFmass(t1->GetP(),
300                                           t1->GetIntegratedLength(),
301                                           t1->GetTOFsignal());
302           } else {
303             tofmass[1] = -1000.;
304           }
305
306           theD0.SetPIDresponse(esdpid0,esdpid1);
307           theD0.SetTOFmasses(tofmass);
308
309           // fill the tree
310           ioD0toKpi=&theD0;
311           treeD0->Fill();
312
313           nD0rec++; nD0rec1ev++;
314           ioD0toKpi=0;  
315         }
316         
317         negtrack = 0;
318       } // loop on negative tracks
319       postrack = 0;
320     }   // loop on positive tracks
321     
322     delete [] trkEntryP;
323     delete [] trkEntryN;
324     delete trkTree;
325
326     printf(" Number of D0 candidates: %d\n",nD0rec1ev);
327   }    // loop on events in file
328
329
330   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
331   printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
332
333   delete vertexer1;
334
335   esdFile->Close();
336
337   // create a copy of this class to be written to output file
338   AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
339
340   // add PDG codes to decay tracks in found candidates (in simulation mode)
341   // and store tree in the output file
342   if(!fSim) {
343     TFile *outroot = new TFile(outName1.Data(),"recreate");
344     treeD0->Write();
345     copy->Write();
346     outroot->Close();
347     delete outroot;
348   } else {
349     printf(" Now adding information from simulation (PDG codes) ...\n");
350     TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
351     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 void AliD0toKpiAnalysis::PrintStatus() const {
369   // Print parameters being used
370
371   printf("Preselections:\n");
372   printf("    fPtCut   = %f GeV\n",fPtCut);
373   printf("    fd0Cut   = %f micron\n",fd0Cut);
374   printf("    fMassCut = %f GeV\n",fMassCut);
375   if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
376   if(fSim) { 
377     printf("Simulation mode\n");
378     if(fOnlySignal) printf("  Only signal goes to file\n");
379   }
380   printf("Cuts on candidates:\n");
381   printf("    |M-MD0| [GeV]    < %f\n",fD0Cuts[0]);
382   printf("    dca    [micron]  < %f\n",fD0Cuts[1]);
383   printf("    cosThetaStar     < %f\n",fD0Cuts[2]);
384   printf("    pTK     [GeV]    > %f\n",fD0Cuts[3]);
385   printf("    pTpi    [GeV]    > %f\n",fD0Cuts[4]);
386   printf("    |d0K|  [micron]  < %f\n",fD0Cuts[5]);
387   printf("    |d0pi| [micron]  < %f\n",fD0Cuts[6]);
388   printf("    d0d0  [micron^2] < %f\n",fD0Cuts[7]);
389   printf("    cosThetaPoint    > %f\n",fD0Cuts[8]);
390
391   return;
392 }
393 //-----------------------------------------------------------------------------
394 Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
395   // Apply preselection in the invariant mass of the pair
396
397   Double_t mD0 = 1.8645;
398   Double_t mPi = 0.13957;
399   Double_t mKa = 0.49368;
400
401   Double_t energy[2];
402   Double_t mom2[2],momTot2;
403
404   mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
405   mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
406
407   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
408             (p[1]+p[4])*(p[1]+p[4])+
409             (p[2]+p[5])*(p[2]+p[5]);
410
411   // D0 -> K- pi+
412   energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
413   energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
414
415   Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
416     
417   // D0bar -> K+ pi-
418   energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
419   energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
420
421   Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
422
423   if(TMath::Abs(minvD0-mD0)    < fMassCut) return kTRUE;
424   if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
425   return kFALSE;
426 }
427 //-----------------------------------------------------------------------------
428 void AliD0toKpiAnalysis::SelectTracks(AliESD *event,
429         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
430         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
431   // Create two TObjArrays with positive and negative tracks and 
432   // apply single-track preselection
433
434   nTrksP=0,nTrksN=0;
435
436   Int_t entr = event->GetNumberOfTracks();
437  
438   // transfer ITS tracks from ESD to arrays and to a tree
439   for(Int_t i=0; i<entr; i++) {
440
441     AliESDtrack *esdtrack = event->GetTrack(i);
442     UInt_t status = esdtrack->GetStatus();
443
444     if(!(status&AliESDtrack::kITSin)) continue;
445
446     // single track selection
447     if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
448
449     if(esdtrack->GetSign()<0) { // negative track
450       trksN.AddLast(esdtrack);
451       trkEntryN[nTrksN] = i;
452       nTrksN++;
453     } else {                 // positive track
454       trksP.AddLast(esdtrack);
455       trkEntryP[nTrksP] = i;
456       nTrksP++;
457     }
458
459   } // loop on ESD tracks
460
461   return;
462 }
463 //-----------------------------------------------------------------------------
464 void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
465                                    Double_t cut2,Double_t cut3,Double_t cut4,
466                                    Double_t cut5,Double_t cut6,
467                                    Double_t cut7,Double_t cut8) {
468   // Set the cuts for D0 selection
469   fD0Cuts[0] = cut0;
470   fD0Cuts[1] = cut1;
471   fD0Cuts[2] = cut2;
472   fD0Cuts[3] = cut3;
473   fD0Cuts[4] = cut4;
474   fD0Cuts[5] = cut5;
475   fD0Cuts[6] = cut6;
476   fD0Cuts[7] = cut7;
477   fD0Cuts[8] = cut8;
478
479   return;
480 }
481 //-----------------------------------------------------------------------------
482 void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
483   // Set the cuts for D0 selection
484
485   for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
486
487   return;
488 }
489 //-----------------------------------------------------------------------------
490 Bool_t 
491 AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
492   // Check if track passes some kinematical cuts  
493   // Magnetic field "b" (kG)
494
495   if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut) 
496     return kFALSE;
497   if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut) 
498     return kFALSE;
499
500   return kTRUE;
501 }
502 //----------------------------------------------------------------------------
503 void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice,
504                                            Int_t evFirst,Int_t evLast) const {
505   // Create a file with simulation info for the reconstructed tracks
506   
507   TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
508   TFile *esdFile = TFile::Open("AliESDs.root");
509
510   AliMC *mc = gAlice->GetMCApp();
511   
512   Int_t      label;
513   TParticle *part;  
514   TParticle *mumpart;
515   REFTRACK   reftrk;
516   
517   AliESD* event = new AliESD;
518   TTree* tree = (TTree*) esdFile->Get("esdTree");
519   tree->SetBranchAddress("ESD",&event);
520   // loop on events in file
521   for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
522     if(iEvent>evLast) break;
523     tree->GetEvent(iEvent);
524     Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
525
526     gAlice->GetEvent(ev);
527
528     Int_t nentr=(Int_t)event->GetNumberOfTracks();
529
530     // Tree for true track parameters
531     char ttname[100];
532     sprintf(ttname,"Tree_Ref_%d",ev);
533     TTree *reftree = new TTree(ttname,"Tree with true track params");
534     reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
535
536     for(Int_t i=0; i<nentr; i++) {
537       AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
538       label = TMath::Abs(esdtrack->GetLabel());
539
540       part = (TParticle*)mc->Particle(label); 
541       reftrk.lab = label;
542       reftrk.pdg = part->GetPdgCode();
543       reftrk.mumlab = part->GetFirstMother();
544       if(part->GetFirstMother()>=0) {
545         mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
546         reftrk.mumpdg = mumpart->GetPdgCode();
547         reftrk.mumprongs = mumpart->GetNDaughters();
548       } else {
549         reftrk.mumpdg=-1;
550         reftrk.mumprongs=-1;
551       }
552       reftrk.Vx = part->Vx();
553       reftrk.Vy = part->Vy();
554       reftrk.Vz = part->Vz();
555       reftrk.Px = part->Px();
556       reftrk.Py = part->Py();
557       reftrk.Pz = part->Pz();
558       
559       reftree->Fill();
560
561     } // loop on tracks   
562
563     outFile->cd();
564     reftree->Write();
565
566     delete reftree;
567   } // loop on events
568
569   esdFile->Close();
570   outFile->Close();
571
572   return;
573 }
574 //-----------------------------------------------------------------------------
575 void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
576   // add pdg codes to candidate decay tracks (for sim)
577
578   TString refFileName("D0TracksRefFile.root");
579   if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) { 
580     printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n"); 
581     return;
582   }
583   TFile *refFile = TFile::Open(refFileName.Data());
584
585   Char_t refTreeName[100];
586   Int_t  event;
587   Int_t  pdg[2],mumpdg[2],mumlab[2];
588   REFTRACK reftrk;
589
590   // read-in reference tree for event 0 (the only event for Pb-Pb)
591   sprintf(refTreeName,"Tree_Ref_%d",0);
592   TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
593   refTree0->SetBranchAddress("rectracks",&reftrk);
594
595   AliD0toKpi *theD0 = 0; 
596   treeD0in->SetBranchAddress("D0toKpi",&theD0);
597   treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
598
599   Int_t entries = (Int_t)treeD0in->GetEntries();
600
601   for(Int_t i=0; i<entries; i++) {
602     if(i%100==0) printf("  done %d candidates of %d\n",i,entries);    
603
604     treeD0in->GetEvent(i);
605     event = theD0->EventNo();
606
607     if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
608       refTree0->GetEvent(theD0->GetTrkNum(0));
609       pdg[0]    = reftrk.pdg;
610       mumpdg[0] = reftrk.mumpdg;
611       mumlab[0] = reftrk.mumlab;
612       refTree0->GetEvent(theD0->GetTrkNum(1));
613       pdg[1]    = reftrk.pdg;
614       mumpdg[1] = reftrk.mumpdg;
615       mumlab[1] = reftrk.mumlab;
616     } else {
617       sprintf(refTreeName,"Tree_Ref_%d",event);
618       TTree *refTree = (TTree*)refFile->Get(refTreeName);
619       refTree->SetBranchAddress("rectracks",&reftrk);
620       refTree->GetEvent(theD0->GetTrkNum(0));
621       pdg[0]    = reftrk.pdg;
622       mumpdg[0] = reftrk.mumpdg;
623       mumlab[0] = reftrk.mumlab;
624       refTree->GetEvent(theD0->GetTrkNum(1));
625       pdg[1]    = reftrk.pdg;
626       mumpdg[1] = reftrk.mumpdg;
627       mumlab[1] = reftrk.mumlab;
628       delete refTree;
629     }
630     
631     theD0->SetPdgCodes(pdg);
632     theD0->SetMumPdgCodes(mumpdg);
633     
634     if(TMath::Abs(mumpdg[0])==421 && 
635        TMath::Abs(mumpdg[1])==421 && 
636        mumlab[0]==mumlab[1] &&
637        reftrk.mumprongs==2 && 
638        ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||  
639         (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
640        ) theD0->SetSignal();
641     
642     if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
643
644   }
645
646   delete refTree0;
647
648   refFile->Close();
649
650   return;
651 }
652
653
654
655
656
657