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