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