GetESD.sh becomes runtest.sh
[u/mrichter/AliRoot.git] / ANALYSIS / AliD0toKpiAnalysis.cxx
index e3a244c03c3bc8e4a14c9efb0b1086ebcd6925b5..e9294d2dd6c455b52076f888f1ae50f32a04c159 100644 (file)
@@ -18,7 +18,7 @@
 // Note: the two decay tracks are labelled: 0 (positive track)
 //                                          1 (negative track)
 // An example of usage can be found in the macro AliD0toKpiTest.C
-//            Origin: A. Dainese    andrea.dainese@pd.infn.it            
+//            Origin: A. Dainese    andrea.dainese@lnl.infn.it            
 //----------------------------------------------------------------------------
 #include <TKey.h>
 #include <TFile.h>
 #include <TSystem.h>
 #include <TParticle.h>
 #include "AliESD.h"
-#include "AliStack.h"
+#include "AliMC.h"
+#include "AliRun.h"
 #include "AliRunLoader.h"
-#include "AliITSVertexerTracks.h"
+#include "AliVertexerTracks.h"
 #include "AliESDVertex.h"
 #include "AliESDv0.h"
 #include "AliD0toKpi.h"
@@ -41,6 +42,7 @@ typedef struct {
   Int_t pdg;
   Int_t mumlab;
   Int_t mumpdg;
+  Int_t mumprongs;
   Float_t Vx,Vy,Vz;
   Float_t Px,Py,Pz;
 } REFTRACK;
@@ -50,6 +52,7 @@ ClassImp(AliD0toKpiAnalysis)
 //----------------------------------------------------------------------------
 AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
   // Default constructor
+
   SetPtCut();
   Setd0Cut();
   SetMassCut();
@@ -131,251 +134,24 @@ Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
 
   return mom*a;
 }
-/*
 //----------------------------------------------------------------------------
 void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
                                        const Char_t *outName) {
   // Find D0 candidates and calculate parameters
 
-  TString trkName("AliITStracksV2.root");
-  if(gSystem->AccessPathName(trkName.Data(),kFileExists)) {
-    printf("AliD0toKpiAnalysis::FindCandidates(): No tracks file found!\n"); 
-    return;
-  }
-
-  TString outName1=outName;
-  TString outName2("nTotEvents.dat");
-
-  Int_t    ev;
-  Int_t    nTotEv=0,nD0rec=0,nD0rec1ev=0;
-  Double_t dca;
-  Double_t v2[3],mom[6],d0[2];
-  Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
-  Int_t    iTrkP,iTrkN,trkEntries;
-  Int_t    nTrksP=0,nTrksN=0;
-  Int_t    trkNum[2];
-  Int_t    okD0=0,okD0bar=0;
-  Char_t   trkTreeName[100],vtxName[100];
-  AliITStrackV2 *postrack = 0;
-  AliITStrackV2 *negtrack = 0;
-
-  // create the AliITSVertexerTracks object
-  // (it will be used only if fVertexOnTheFly=kTrue)
-  AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
-  vertexer1->SetMinTracks(2);
-  Int_t  skipped[2];
-  Bool_t goodVtx1;
-  
-
-  // define the cuts for vertexing
-  Double_t vtxcuts[]={50., // max. allowed chi2
-                     0.0, // min. allowed negative daughter's impact param 
-                     0.0, // min. allowed positive daughter's impact param 
-                     1.0, // max. allowed DCA between the daughter tracks
-                    -1.0, // min. allowed cosine of pointing angle
-                     0.0, // min. radius of the fiducial volume
-                     2.9};// max. radius of the fiducial volume
-  
-  // create the AliV0vertexer object
-  AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
-
-  // create tree for reconstructed D0s
-  AliD0toKpi *ioD0toKpi=0;
-  TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
-  treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
-
-  // open file with tracks
-  TFile *trkFile = TFile::Open(trkName.Data());
-
-  // loop on events in file
-  for(ev=evFirst; ev<=evLast; ev++) {
-    printf(" --- Looking for D0->Kpi in event  %d ---\n",ev);
-
-    // retrieve primary vertex from file
-    sprintf(vtxName,"VertexTracks_%d",ev);
-    AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
-    vertex1stored->GetXYZ(fV1);
-    delete vertex1stored;
-
-    // retrieve tracks from file
-    sprintf(trkTreeName,"TreeT_ITS_%d",ev);
-    TTree *trkTree=(TTree*)trkFile->Get(trkTreeName);
-    if(!trkTree) { 
-      printf("AliD0toKpiAnalysis::FindCandidates():\n  tracks tree not found for evet %d\n",ev); 
-      continue; 
-    }
-    trkEntries = (Int_t)trkTree->GetEntries();
-    printf(" Number of tracks: %d\n",trkEntries);
-
-    // count the total number of events
-    nTotEv++;
-
-    // call function which applies sigle-track selection and
-    // separetes positives and negatives
-    TObjArray trksP(trkEntries/2);
-    Int_t *trkEntryP = new Int_t[trkEntries];
-    TObjArray trksN(trkEntries/2);
-    Int_t *trkEntryN = new Int_t[trkEntries];
-    SelectTracks(*trkTree,trksP,trkEntryP,nTrksP,trksN,trkEntryN,nTrksN);
-      
-    nD0rec1ev = 0;
-
-    // loop on positive tracks
-    for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
-      if(iTrkP%10==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
-         
-      // get track from track array
-      postrack = (AliITStrackV2*)trksP.At(iTrkP);
-      trkNum[0] = trkEntryP[iTrkP];      
-
-      // loop on negative tracks 
-      for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
-       // get track from tracks array
-       negtrack = (AliITStrackV2*)trksN.At(iTrkN);
-       trkNum[1] = trkEntryN[iTrkN];      
-
-       AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
-
-       //
-       // ----------- DCA MINIMIZATION ------------------
-       //
-       // find the DCA and propagate the tracks to the DCA 
-       dca = vertexer2->PropagateToDCA(pnt,ppt);
-
-       // define the AliV0vertex object
-       AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
-         
-       // get position of the secondary vertex
-       vertex2->GetXYZ(v2[0],v2[1],v2[2]);
-       
-       delete vertex2;
-  
-       // momenta of the tracks at the vertex
-       ptP = 1./TMath::Abs(ppt->Get1Pt());
-       alphaP = ppt->GetAlpha();
-       phiP = alphaP+TMath::ASin(ppt->GetSnp());
-       mom[0] = ptP*TMath::Cos(phiP); 
-       mom[1] = ptP*TMath::Sin(phiP);
-       mom[2] = ptP*ppt->GetTgl();
-         
-       ptN = 1./TMath::Abs(pnt->Get1Pt());
-       alphaN = pnt->GetAlpha();
-       phiN = alphaN+TMath::ASin(pnt->GetSnp());
-       mom[3] = ptN*TMath::Cos(phiN); 
-       mom[4] = ptN*TMath::Sin(phiN);
-       mom[5] = ptN*pnt->GetTgl();
-         
-       goodVtx1 = kTRUE;
-       // no vertexing if DeltaMass > fMassCut 
-       if(fVertexOnTheFly) {
-         goodVtx1 = kFALSE;
-         if(SelectInvMass(mom)) {
-           // primary vertex from *other* tracks in event
-           vertexer1->SetVtxStart(fV1[0],fV1[1]);
-           skipped[0] = trkEntryP[iTrkP];
-           skipped[1] = trkEntryN[iTrkN];
-           vertexer1->SetSkipTracks(2,skipped);
-           AliESDVertex *vertex1onfly = 
-             (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree); 
-           if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
-           vertex1onfly->GetXYZ(fV1);
-           //vertex1onfly->PrintStatus();
-           delete vertex1onfly;        
-         }
-       }         
-
-       // impact parameters of the tracks w.r.t. the primary vertex
-       d0[0] =  10000.*ppt->GetD(fV1[0],fV1[1]);
-       d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
-
-       // create the object AliD0toKpi
-       AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
-
-       // select D0s
-       if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
-             
-         // fill the tree
-         ioD0toKpi=&theD0;
-         treeD0->Fill();
-
-         nD0rec++; nD0rec1ev++;
-         ioD0toKpi=0;  
-       }
-       
-       negtrack = 0;
-      } // loop on negative tracks
-      postrack = 0;
-    }   // loop on positive tracks
-
-    trksP.Delete();
-    trksN.Delete();
-    delete [] trkEntryP;
-    delete [] trkEntryN;
-    delete trkTree;
-
-    printf(" Number of D0 candidates: %d\n",nD0rec1ev);
-  }    // loop on events in file
-
-
-  printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
-  printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
-
-  delete vertexer1;
-  delete vertexer2;
-
-  trkFile->Close();
-
-  // create a copy of this class to be written to output file
-  AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
-  copy->PrintStatus();
-
-  // add PDG codes to decay tracks in found candidates (in simulation mode)
-  // and store tree in the output file
-  if(!fSim) {
-    TFile *outroot = new TFile(outName1.Data(),"recreate");
-    treeD0->Write();
-    copy->Write();
-    outroot->Close();
-    delete outroot;
-  } else {
-    printf(" Now adding information from simulation (PDG codes) ...\n");
-    TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
-    MakeTracksRefFile(evFirst,evLast);
-    SimulationInfo(treeD0,treeD0sim);
-    delete treeD0;
-    TFile *outroot = new TFile(outName1.Data(),"recreate");
-    treeD0sim->Write();
-    copy->Write();
-    outroot->Close();
-    delete outroot;
-  }
-
-  // write to a file the total number of events
-  FILE *outdat = fopen(outName2.Data(),"w");
-  fprintf(outdat,"%d\n",nTotEv);
-  fclose(outdat);
-
-  return;
-}
-*/
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
-                                          const Char_t *outName) {
-  // Find D0 candidates and calculate parameters
 
-  TString esdName("AliESDs.root");
+  TString esdName="AliESDs.root";
   if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
     printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n"); 
     return;
   }
 
   TString outName1=outName;
-  TString outName2("nTotEvents.dat");
+  TString outName2="nTotEvents.dat";
 
   Int_t    nTotEv=0,nD0rec=0,nD0rec1ev=0;
   Double_t dca;
   Double_t v2[3],mom[6],d0[2];
-  //Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
   Int_t    iTrkP,iTrkN,trkEntries;
   Int_t    nTrksP=0,nTrksN=0;
   Int_t    trkNum[2];
@@ -384,27 +160,20 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
   AliESDtrack *postrack = 0;
   AliESDtrack *negtrack = 0;
 
-  // create the AliITSVertexerTracks object
+  // create the AliVertexerTracks object
   // (it will be used only if fVertexOnTheFly=kTrue)
-  AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
-  vertexer1->SetMinTracks(2);
+  AliVertexerTracks *vertexer1 = new AliVertexerTracks;
+  if(fVertexOnTheFly) {
+    // open the mean vertex
+    TFile *invtx = new TFile("AliESDVertexMean.root");
+    AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+    invtx->Close();
+    vertexer1->SetVtxStart(initVertex);
+    delete invtx;
+  }
   Int_t  skipped[2];
   Bool_t goodVtx1;
   
-  /*
-  // define the cuts for vertexing
-  Double_t vtxcuts[]={50., // max. allowed chi2
-                     0.0, // min. allowed negative daughter's impact param 
-                     0.0, // min. allowed positive daughter's impact param 
-                     1.0, // max. allowed DCA between the daughter tracks
-                    -1.0, // min. allowed cosine of pointing angle
-                     0.0, // min. radius of the fiducial volume
-                     2.9};// max. radius of the fiducial volume
-  
-  // create the AliV0vertexer object
-  AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
-  */
-
   // create tree for reconstructed D0s
   AliD0toKpi *ioD0toKpi=0;
   TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
@@ -414,31 +183,31 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
   TFile *esdFile = TFile::Open(esdName.Data());
   AliESD* event = new AliESD;
   TTree* tree = (TTree*) esdFile->Get("esdTree");
-  if (!tree) {
+  if(!tree) {
     Error("FindCandidatesESD", "no ESD tree found");
     return;
   }
-  tree->SetBranchAddress("ESD", &event);
+  tree->SetBranchAddress("ESD",&event);
 
   // loop on events in file
-  for (Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
-    if (iEvent > evLast) break;
+  for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
+    if(iEvent > evLast) break;
     tree->GetEvent(iEvent);
-    Int_t ev = (Int_t)event->GetEventNumber();
+    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.
     printf("--- Finding D0 -> Kpi in event %d\n",ev);
+    // count the total number of events
+    nTotEv++;
 
-    // retrieve primary vertex from file
-    //sprintf(vtxName,"Vertex_%d",ev);
-    //AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
-    //vertex1stored->GetXYZ(fV1);
-    //delete vertex1stored;
-    event->GetVertex()->GetXYZ(fV1);
-
-    trkEntries = event->GetNumberOfTracks();
+    trkEntries = (Int_t)event->GetNumberOfTracks();
     printf(" Number of tracks: %d\n",trkEntries);
+    if(trkEntries<2) continue;
 
-    // count the total number of events
-    nTotEv++;
+    // retrieve primary vertex from file
+    if(!event->GetPrimaryVertex()) { 
+      printf(" No vertex\n");
+      continue;
+    }
+    event->GetPrimaryVertex()->GetXYZ(fV1);
 
     // call function which applies sigle-track selection and
     // separetes positives and negatives
@@ -447,34 +216,29 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
     TObjArray trksN(trkEntries/2);
     Int_t    *trkEntryN   = new Int_t[trkEntries];
     TTree *trkTree = new TTree();
-    if(fVertexOnTheFly) {
-      SelectTracksESDvtx(*event,trkTree,trksP,trkEntryP,nTrksP,
-                                       trksN,trkEntryN,nTrksN);
-    } else {
-      SelectTracksESD(*event,trksP,trkEntryP,nTrksP,
-                            trksN,trkEntryN,nTrksN);
-    }      
+    SelectTracks(event,trksP,trkEntryP,nTrksP,
+                      trksN,trkEntryN,nTrksN);
+
+    printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
 
-    AliDebugClass(1,Form(" pos. tracks: %d    neg .tracks: %d",nTrksP,nTrksN));
 
     nD0rec1ev = 0;
 
-    // loop on positive tracks
+    // LOOP ON  POSITIVE  TRACKS
     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
-      if(iTrkP%10==0) AliDebugClass(1,Form("  Processing positive track number %d of %d",iTrkP,nTrksP));
+      if(iTrkP%1==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
          
       // get track from track array
       postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
       trkNum[0] = trkEntryP[iTrkP];      
 
-      // loop on negative tracks 
+      // LOOP ON  NEGATIVE  TRACKS
       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 
        // get track from tracks array
        negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
        trkNum[1] = trkEntryN[iTrkN];      
 
-        {
        //
        // ----------- DCA MINIMIZATION ------------------
        //
@@ -491,53 +255,32 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
         vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
         vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
        // impact parameters of the tracks w.r.t. the primary vertex
-
        d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
        d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
-       }
-        /*
-       // momenta of the tracks at the vertex
-        //Double_t x,par[5]; postrack->GetExternalParameters(x,par); 
-       //ptP = 1./TMath::Abs(par[4]);
-       //alphaP = postrack->GetAlpha();
-       //phiP = alphaP+TMath::ASin(par[2]);
-         postrack->GetPxPyPz();
-       mom[0] = ptP*TMath::Cos(phiP); 
-       mom[1] = ptP*TMath::Sin(phiP);
-       mom[2] = ptP*par[3];
-         
-       ptN = 1./TMath::Abs(pnt->Get1Pt());
-       alphaN = pnt->GetAlpha();
-       phiN = alphaN+TMath::ASin(pnt->GetSnp());
-       mom[3] = ptN*TMath::Cos(phiN); 
-       mom[4] = ptN*TMath::Sin(phiN);
-       mom[5] = ptN*pnt->GetTgl();
-       */
        goodVtx1 = kTRUE;
+       
        // no vertexing if DeltaMass > fMassCut 
        if(fVertexOnTheFly) {
          goodVtx1 = kFALSE;
          if(SelectInvMass(mom)) {
-           // primary vertex from *other* tracks in event
-           vertexer1->SetVtxStart(fV1[0],fV1[1]);
+           // primary vertex from *other* tracks in the event
            skipped[0] = trkEntryP[iTrkP];
            skipped[1] = trkEntryN[iTrkN];
            vertexer1->SetSkipTracks(2,skipped);
            AliESDVertex *vertex1onfly = 
-             (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree); 
+             (AliESDVertex*)vertexer1->FindPrimaryVertex(event); 
            if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
            vertex1onfly->GetXYZ(fV1);
            //vertex1onfly->PrintStatus();
            delete vertex1onfly;        
          }
-       }         
+       }
+       
 
        // create the object AliD0toKpi
        AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
-
        // select D0s
        if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
-             
          // get PID info from ESD
          AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
          Double_t esdpid0[5];
@@ -575,14 +318,11 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
       } // loop on negative tracks
       postrack = 0;
     }   // loop on positive tracks
-
-    trksP.Delete();
-    trksN.Delete();
+    
     delete [] trkEntryP;
     delete [] trkEntryN;
     delete trkTree;
 
-
     printf(" Number of D0 candidates: %d\n",nD0rec1ev);
   }    // loop on events in file
 
@@ -591,7 +331,6 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
   printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
 
   delete vertexer1;
-  //delete vertexer2;
 
   esdFile->Close();
 
@@ -609,7 +348,6 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
   } else {
     printf(" Now adding information from simulation (PDG codes) ...\n");
     TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
-    MakeTracksRefFileESD();
     SimulationInfo(treeD0,treeD0sim);
     delete treeD0;
     TFile *outroot = new TFile(outName1.Data(),"recreate");
@@ -686,46 +424,8 @@ Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
   if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
   return kFALSE;
 }
-/*
 //-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracks(TTree &trkTree,
-                 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
-                 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
-  // Create two TObjArrays with positive and negative tracks and 
-  // apply single-track preselection
-
-  nTrksP=0,nTrksN=0;
-  Int_t entr = (Int_t)trkTree.GetEntries();
-
-  // trasfer tracks from tree to arrays
-  for(Int_t i=0; i<entr; i++) {
-
-    AliITStrackV2 *track = new AliITStrackV2; 
-    trkTree.SetBranchAddress("tracks",&track);
-
-    trkTree.GetEvent(i);
-
-    // single track selection
-    if(!SingleTrkCuts(*track)) { delete track; continue; }
-
-    if(track->Get1Pt()>0.) { // negative track
-      trksN.AddLast(track);
-      trkEntryN[nTrksN] = i;
-      nTrksN++;
-    } else {                 // positive track
-      trksP.AddLast(track);
-      trkEntryP[nTrksP] = i;
-      nTrksP++;
-    }
-
-  }
-
-  return;
-}
-*/
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
+void AliD0toKpiAnalysis::SelectTracks(AliESD *event,
         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
   // Create two TObjArrays with positive and negative tracks and 
@@ -733,18 +433,18 @@ void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
 
   nTrksP=0,nTrksN=0;
 
-  Int_t entr = event.GetNumberOfTracks();
+  Int_t entr = event->GetNumberOfTracks();
  
   // transfer ITS tracks from ESD to arrays and to a tree
   for(Int_t i=0; i<entr; i++) {
 
-    AliESDtrack *esdtrack = event.GetTrack(i);
+    AliESDtrack *esdtrack = event->GetTrack(i);
     UInt_t status = esdtrack->GetStatus();
 
-    if(!(status&AliESDtrack::kITSrefit)) continue;
+    if(!(status&AliESDtrack::kITSin)) continue;
 
     // single track selection
-    if(!SingleTrkCuts(*esdtrack,event.GetMagneticField())) continue;
+    if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
 
     if(esdtrack->GetSign()<0) { // negative track
       trksN.AddLast(esdtrack);
@@ -761,54 +461,6 @@ void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
   return;
 }
 //-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracksESDvtx(AliESD &event,TTree *trkTree,
-        TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
-        TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
-  // Create two TObjArrays with positive and negative tracks and 
-  // apply single-track preselection
-
-  nTrksP=0,nTrksN=0;
-
-  Int_t entr = event.GetNumberOfTracks();
-  AliESDtrack *esdtrackfortree = 0;
-  trkTree->Branch("tracks","AliESDtrack",&esdtrackfortree,entr,0);
-
-
-  // transfer the tracks from ESD to arrays and to a tree
-  for(Int_t i=0; i<entr; i++) {
-
-    AliESDtrack *esdtrack = event.GetTrack(i);
-    UInt_t status = esdtrack->GetStatus();
-
-    if(!(status&AliESDtrack::kITSrefit)) continue;
-
-    // store track in the tree to be used for primary vertex finding
-    esdtrackfortree = new AliESDtrack(*esdtrack);
-    trkTree->Fill();
-    //itstrackfortree = 0;
-    delete esdtrackfortree;
-
-    // single track selection
-    if(!SingleTrkCuts(*esdtrack,event.GetMagneticField())) continue;
-
-    if(esdtrack->GetSign()<0) { // negative track
-      trksN.AddLast(esdtrack);
-      trkEntryN[nTrksN] = i;
-      nTrksN++;
-    } else {                 // positive track
-      trksP.AddLast(esdtrack);
-      trkEntryP[nTrksP] = i;
-      nTrksP++;
-    }
-
-  } // loop on esd tracks
-
-  //delete itstrackfortree;
-
-  return;
-}
-//-----------------------------------------------------------------------------
 void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
                                   Double_t cut2,Double_t cut3,Double_t cut4,
                                   Double_t cut5,Double_t cut6,
@@ -847,110 +499,31 @@ AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
 
   return kTRUE;
 }
-/*
 //----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::MakeTracksRefFile(Int_t evFirst,Int_t evLast) 
-  const {
+void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice,
+                                          Int_t evFirst,Int_t evLast) const {
   // Create a file with simulation info for the reconstructed tracks
   
-  TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
-  TFile *trk = TFile::Open("AliITStracksV2.root");
-  AliRunLoader *rl = AliRunLoader::Open("galice.root");
-
-  // load kinematics
-  rl->LoadKinematics();
-  
-  Int_t      label;
-  TParticle *part;  
-  TParticle *mumpart;
-  REFTRACK   reftrk;
-  
-  for(Int_t ev=evFirst; ev<=evLast; ev++){
-    rl->GetEvent(ev);  
-    AliStack *stack = rl->Stack();
-
-    trk->cd();
-
-    // Tree with tracks
-    char tname[100];
-    sprintf(tname,"TreeT_ITS_%d",ev);
-
-    TTree *tracktree=(TTree*)trk->Get(tname);
-    if(!tracktree) continue;
-    AliITStrackV2 *track = new AliITStrackV2; 
-    tracktree->SetBranchAddress("tracks",&track);
-    Int_t nentr=(Int_t)tracktree->GetEntries();
-
-    // Tree for true track parameters
-    char ttname[100];
-    sprintf(ttname,"Tree_Ref_%d",ev);
-    TTree *reftree = new TTree(ttname,"Tree with true track params");
-    reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
-
-    for(Int_t i=0; i<nentr; i++) {
-      tracktree->GetEvent(i);
-      label = TMath::Abs(track->GetLabel());
-
-      part = (TParticle*)stack->Particle(label);
-      reftrk.lab = label;
-      reftrk.pdg = part->GetPdgCode();
-      reftrk.mumlab = part->GetFirstMother();
-      if(part->GetFirstMother()>=0) {
-       mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
-       reftrk.mumpdg = mumpart->GetPdgCode();
-      } else {
-       reftrk.mumpdg=-1;
-      }
-      reftrk.Vx = part->Vx();
-      reftrk.Vy = part->Vy();
-      reftrk.Vz = part->Vz();
-      reftrk.Px = part->Px();
-      reftrk.Py = part->Py();
-      reftrk.Pz = part->Pz();
-      
-      reftree->Fill();
-    } // loop on tracks   
-
-    out->cd();
-    reftree->Write();
-
-    delete track;
-    delete reftree;
-    delete tracktree;
-    delete stack;
-  } // loop on events
-
-  trk->Close();
-  out->Close();
-
-  return;
-}
-*/
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
-  // Create a file with simulation info for the reconstructed tracks
-  
-  TFile *outFile = TFile::Open("ITStracksRefFile.root","recreate");
+  TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
   TFile *esdFile = TFile::Open("AliESDs.root");
-  AliRunLoader *rl = AliRunLoader::Open("galice.root");
 
-  // load kinematics
-  rl->LoadKinematics();
+  AliMC *mc = gAlice->GetMCApp();
   
   Int_t      label;
   TParticle *part;  
   TParticle *mumpart;
   REFTRACK   reftrk;
   
-  TKey *key=0;
-  TIter next(esdFile->GetListOfKeys());
+  AliESD* event = new AliESD;
+  TTree* tree = (TTree*) esdFile->Get("esdTree");
+  tree->SetBranchAddress("ESD",&event);
   // loop on events in file
-  while ((key=(TKey*)next())!=0) {
-    AliESD *event=(AliESD*)key->ReadObj();
-    Int_t ev = (Int_t)event->GetEventNumber();
+  for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
+    if(iEvent>evLast) break;
+    tree->GetEvent(iEvent);
+    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.
 
-    rl->GetEvent(ev);  
-    AliStack *stack = rl->Stack();
+    gAlice->GetEvent(ev);
 
     Int_t nentr=(Int_t)event->GetNumberOfTracks();
 
@@ -964,15 +537,17 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
       AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
       label = TMath::Abs(esdtrack->GetLabel());
 
-      part = (TParticle*)stack->Particle(label);
+      part = (TParticle*)mc->Particle(label); 
       reftrk.lab = label;
       reftrk.pdg = part->GetPdgCode();
       reftrk.mumlab = part->GetFirstMother();
       if(part->GetFirstMother()>=0) {
-       mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
+       mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
        reftrk.mumpdg = mumpart->GetPdgCode();
+       reftrk.mumprongs = mumpart->GetNDaughters();
       } else {
        reftrk.mumpdg=-1;
+       reftrk.mumprongs=-1;
       }
       reftrk.Vx = part->Vx();
       reftrk.Vy = part->Vy();
@@ -989,8 +564,6 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
     reftree->Write();
 
     delete reftree;
-    delete event;
-    delete stack;
   } // loop on events
 
   esdFile->Close();
@@ -1002,7 +575,7 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
 void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
   // add pdg codes to candidate decay tracks (for sim)
 
-  TString refFileName("ITStracksRefFile.root");
+  TString refFileName("D0TracksRefFile.root");
   if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) { 
     printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n"); 
     return;
@@ -1058,8 +631,13 @@ void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const
     theD0->SetPdgCodes(pdg);
     theD0->SetMumPdgCodes(mumpdg);
     
-    if(TMath::Abs(mumpdg[0])==421 && TMath::Abs(mumpdg[1])==421 
-       && mumlab[0]==mumlab[1]) theD0->SetSignal();
+    if(TMath::Abs(mumpdg[0])==421 && 
+       TMath::Abs(mumpdg[1])==421 && 
+       mumlab[0]==mumlab[1] &&
+       reftrk.mumprongs==2 && 
+       ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||  
+       (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
+       ) theD0->SetSignal();
     
     if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();