ANALYSIS/AliD0toKpiAnalysis.cxx

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