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