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 "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 *gAlice,
 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 = gAlice->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     gAlice->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*)gAlice->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