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