HLT/trigger/AliD0vtxFinderSgn_pp_VTX.C

   1 /****************************************************************************
   2  *                                                                          *
   3  * This macro computes the position of the D0 decay vertex using            *
   4  * helix DCA minimization by J.Belikov                                      *
   5  *                                                                          *
   6  * Reconstructed D0 are stored in a tree as AliD0toKpi objects              *
   7  *                                                                          *
   8  ****************************************************************************/
   9 //#define __COMPILE__
  10 //#ifdef __COMPILE__
  11 #if !defined(__CINT__) || defined(__MAKECINT__)
  12 //-- --- standard headers-------------
  13 #include <iostream.h>
  14 //--------Root headers ---------------
  15 #include <TSystem.h>
  16 #include <TFile.h>
  17 #include <TString.h>
  18 #include <TStopwatch.h>
  19 #include <TObject.h>
  20 #include <TVector3.h>
  21 #include <TTree.h>
  22 #include <TParticle.h>
  23 #include <TArray.h>
  24 //----- AliRoot headers ---------------
  25 #include "alles.h"
  26 #include "AliRun.h"
  27 #include "AliKalmanTrack.h"
  28 #include "AliITStrackV2.h"
  29 #include "AliHeader.h"
  30 #include "AliGenEventHeader.h"
  31 #include "AliV0vertex.h"
  32 #include "AliV0vertexer.h"
  33 #include "AliITSVertex.h"
  34 #include "AliITSVertexer.h"
  35 #include "AliITSVertexerTracks.h"
  36 #include "AliD0toKpi.h"
  37 #endif
  38 //-------------------------------------
  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 } RECTRACK;
  48
  49 // field (T)
  50 const Double_t kBz = 0.4;
  51
  52 // primary vertex
  53 Double_t gv1[3] = {0.,0.,0.};
  54
  55 // sigle track cuts
  56 //const Double_t kPtCut = 0.5;  // GeV/c
  57 //const Double_t kd0Cut = 50.; // micron
  58 const Double_t kPtCut = 0.;  // GeV/c
  59 const Double_t kd0Cut = 0.; // micron
  60
  61 // cuts on D0 candidate (to be passed to function AliD0toKpi::Select())
  62 //
  63 // cuts[0] = inv. mass half width [GeV]
  64 // cuts[1] = dca [micron]
  65 // cuts[2] = cosThetaStar
  66 // cuts[3] = pTK [GeV/c]
  67 // cuts[4] = pTPi [GeV/c]
  68 // cuts[5] = d0K [micron]   upper limit!
  69 // cuts[6] = d0Pi [micron]  upper limit!
  70 // cuts[7] = d0d0 [micron^2]
  71 // cuts[8] = cosThetaPoint
  72 const Double_t kD0cuts[9] = {0.012,
  73                              500.,
  74                              0.8,
  75                              0.5,
  76                              0.5,
  77                              500.,
  78                              500.,
  79                             -5000.,
  80                              0.8};
  81 const Double_t kD0cutsAll[9] = {.1,
  82                                 1.e5,
  83                                 1.1,
  84                                 0.,
  85                                 0.,
  86                                 1.e10,
  87                                 1.e10,
  88                                 1.e10,
  89                                 -1.1};
  90 // mass cut for 1st level rejection
  91 const Double_t kMassCut = 0.1; // GeV
  92
  93 // this function ckecks files existence
  94 Bool_t GetInputFiles();
  95
  96 // 1st level rejection based on inv. mass
  97 Bool_t SelectInvMass(const Double_t p[6]);
  98
  99 // this function creates TObjArrays with positive and negative tracks
 100 void   SelectTracks(TTree& itsTree,
 101                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 102                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
 103
 104 // this function applies single track cuts
 105 Bool_t TrkCuts(const AliITStrackV2& trk);
 106
 107 void   AliD0vtxFinderSgn_pp_VTX(Int_t evFirst=0,Int_t evLast=999) {
 108
 109   const Char_t *name="AliD0vtxFinderSgn_pp_VTX";
 110   cerr<<'\n'<<name<<"...\n";
 111   gBenchmark->Start(name);
 112
 113   // load library for D0 candidates class
 114   gSystem->Load("libD0toKpi.so");
 115
 116   // check existence of input files
 117   if(!GetInputFiles()) { cerr<<"No tracks file found"<<endl; return; }
 118
 119   TString outNameSgn("AliD0toKpiSgn.root");
 120   TString outNameBkg("AliD0toKpiBkgS.root");
 121
 122   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
 123
 124   Bool_t isSignal;
 125   Int_t ev;
 126   Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
 127   Double_t dca;
 128   Double_t v2[3],mom[6],d0[2];
 129   Int_t pdg[2],mum[2],mumlab[2];
 130   Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
 131   Int_t l,k,iTrkP,iTrkN,itsEntries;
 132   Int_t mesonD0,mesonD0label[10];
 133   Int_t nTrksP=0,nTrksN=0;
 134   Int_t okD0=0,okD0bar=0;
 135   Char_t vtxName[100],trksName[100],refsName[100];
 136   RECTRACK rectrk;
 137   AliITStrackV2 *postrack = 0;
 138   AliITStrackV2 *negtrack = 0;
 139
 140   // create the AliITSVertexerTracks object
 141   AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
 142   vertexer1->SetMinTracks(2);
 143   vertexer1->SetDebug(0);
 144   Int_t skipped[2];
 145
 146   // define the cuts for vertexing
 147   Double_t vtxcuts[]={33., // max. allowed chi2
 148                       0.0, // min. allowed negative daughter's impact param
 149                       0.0, // min. allowed positive daughter's impact param
 150                       1.0, // max. allowed DCA between the daughter tracks
 151                      -1.0, // min. allowed cosine of V0's pointing angle
 152                       0.0, // min. radius of the fiducial volume
 153                       2.9};// max. radius of the fiducial volume
 154
 155   // create the AliV0vertexer object
 156   AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
 157
 158
 159   // Create trees for reconstructed D0s
 160   AliD0toKpi *ioD0toKpi=0;
 161   // Signal
 162   TTree D0TreeSgn("TreeD0","Tree with D0 candidates");
 163   D0TreeSgn.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
 164   // Background
 165   TTree D0TreeBkg("TreeD0","Tree with D0 candidates");
 166   D0TreeBkg.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
 167
 168   // Open file with ITS tracks
 169   TFile* itstrks = TFile::Open("AliITStracksV2.root");
 170
 171   // Open file with ITS track references
 172   TFile* itsrefs = TFile::Open("ITStracksRefFile.root");
 173
 174
 175
 176   // loop on events in file
 177   for(ev=evFirst; ev<=evLast; ev++) {
 178     printf(" --- Processing event  %d ---\n",ev);
 179     sprintf(trksName,"TreeT_ITS_%d",ev);
 180     sprintf(refsName,"Tree_Ref_%d",ev);
 181
 182
 183
 184     // tracks from ITS
 185     TTree *itsTree=(TTree*)itstrks->Get(trksName);
 186     if(!itsTree) continue;
 187     itsEntries = (Int_t)itsTree->GetEntries();
 188     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
 189
 190     // tree from reference file
 191     TTree *refTree=(TTree*)itsrefs->Get(refsName);
 192     refTree->SetBranchAddress("rectracks",&rectrk);
 193
 194     Int_t refEntries = (Int_t)refTree->GetEntries();
 195     k=0;
 196     for(l=0; l<10; l++) mesonD0label[l]=-1;
 197     for(l=0; l<refEntries; l++) {
 198       refTree->GetEvent(l);
 199       if(TMath::Abs(rectrk.mumpdg)!=421) continue;
 200       //cerr<<" pdg: "<<rectrk.pdg<<"  mumpdg: "<<rectrk.mumpdg<<" mumlabel: "<<rectrk.mumlab<<endl;
 201       mesonD0label[k] = rectrk.mumlab;
 202       k++;
 203     }
 204
 205     mesonD0 = mesonD0label[0];
 206     for(l=0; l<10; l++) {
 207       if(mesonD0label[l]==-1) continue;
 208       for(k=0; k<10; k++) {
 209         if(k==l) continue;
 210         if(mesonD0label[k] == mesonD0label[l]) mesonD0 = mesonD0label[k];
 211       }
 212     }
 213     //cerr<<" mesonD0: "<<mesonD0<<endl;
 214
 215     Double_t *brWgt = new Double_t[refEntries];
 216     for(l=0; l<refEntries; l++) {
 217       refTree->GetEvent(l);
 218       // normally tracks have weight = 1
 219       brWgt[l] = 1.;
 220       // decay products of non-good D0 weighted with D0->Kpi B.R.
 221       if(TMath::Abs(rectrk.mumpdg)==421 && rectrk.mumlab!=mesonD0)
 222         brWgt[l]=0.038;
 223       // decay products of D+ weighted with D+->Kpipi B.R.
 224       if(TMath::Abs(rectrk.mumpdg)==411)
 225         brWgt[l]=0.09;
 226     }
 227
 228     // count the total number of events
 229     nTotEv++;
 230
 231     // call function which applies sigle track selection and
 232     // separetes positives and negatives
 233     TObjArray trksP(itsEntries/2);
 234     Int_t *itsEntryP = new Int_t[itsEntries];
 235     TObjArray trksN(itsEntries/2);
 236     Int_t *itsEntryN = new Int_t[itsEntries];
 237     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
 238
 239     nD0rec1ev = 0;
 240     // loop on positive tracks
 241     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 242       //if(iTrkP % 10==0) cerr<<"--- Processing positive track number: "<<iTrkP<<" of "<<nTrksP<<"\r";
 243
 244       // get track from track array
 245       postrack = (AliITStrackV2*)trksP.At(iTrkP);
 246
 247       // get info on tracks PDG and mothers PDG from reference file
 248       refTree->GetEvent(itsEntryP[iTrkP]);
 249       pdg[0] = rectrk.pdg;
 250       mum[0] = rectrk.mumpdg;
 251       mumlab[0] = rectrk.mumlab;
 252
 253       // loop on negative tracks
 254       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 255         // get track from track array
 256         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
 257
 258         // get info on tracks PDG and mothers PDG from reference file
 259         refTree->GetEvent(itsEntryN[iTrkN]);
 260         pdg[1] = rectrk.pdg;
 261         mum[1] = rectrk.mumpdg;
 262         mumlab[1] = rectrk.mumlab;
 263
 264         //
 265         // ----------- DCA MINIMIZATION ------------------
 266         //
 267         // find the DCA and propagate the tracks to the DCA
 268         dca = vertexer2->PropagateToDCA(negtrack,postrack);
 269
 270         // define the AliV0vertex object
 271         AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
 272
 273         // get position of the vertex
 274         vertex2->GetXYZ(v2[0],v2[1],v2[2]);
 275
 276         delete vertex2;
 277
 278         // momenta of the tracks at the vertex
 279         ptP = 1./TMath::Abs(postrack->Get1Pt());
 280         alphaP = postrack->GetAlpha();
 281         phiP = alphaP+TMath::ASin(postrack->GetSnp());
 282         mom[0] = ptP*TMath::Cos(phiP);
 283         mom[1] = ptP*TMath::Sin(phiP);
 284         mom[2] = ptP*postrack->GetTgl();
 285
 286         ptN = 1./TMath::Abs(negtrack->Get1Pt());
 287         alphaN = negtrack->GetAlpha();
 288         phiN = alphaN+TMath::ASin(negtrack->GetSnp());
 289         mom[3] = ptN*TMath::Cos(phiN);
 290         mom[4] = ptN*TMath::Sin(phiN);
 291         mom[5] = ptN*negtrack->GetTgl();
 292
 293         Bool_t goodVtx = kFALSE;
 294         // no vertexing if DeltaMass > kMassCut
 295         if(SelectInvMass(mom)) {
 296           // primary vertex from other tracks in event
 297           vertexer1->SetVtxStart(0.,0.);
 298           skipped[0] = itsEntryP[iTrkP];
 299           skipped[1] = itsEntryN[iTrkN];
 300           vertexer1->SetSkipTracks(2,skipped);
 301           AliITSVertex *vertex1 =
 302             (AliITSVertex*)vertexer1->VertexOnTheFly(*itsTree);
 303           if(vertex1->GetNContributors()>0) goodVtx = kTRUE;
 304           vertex1->GetXYZ(gv1);
 305           // impact parameters of the tracks w.r.t. the primary vertex
 306           d0[0] =  10000.*postrack->GetD(gv1[0],gv1[1]);
 307           d0[1] = -10000.*negtrack->GetD(gv1[0],gv1[1]);
 308           //vertex1->PrintStatus();
 309           delete vertex1;
 310         }
 311
 312         isSignal = kFALSE;
 313         if(mumlab[0]==mesonD0 && mumlab[1]==mesonD0) isSignal = kTRUE;
 314
 315
 316         // create the object TD0rec and store it in the tree
 317         AliD0toKpi theD0(isSignal,ev,gv1,v2,dca,mom,d0,pdg,mum);
 318
 319         // select D0s
 320         if(goodVtx && theD0.Select(kD0cutsAll,okD0,okD0bar)) {
 321
 322           // compute the weights
 323           theD0.ComputeWgts();
 324           theD0.CorrectWgt4BR(brWgt[itsEntryN[iTrkN]]*brWgt[itsEntryP[iTrkP]]);
 325
 326           // fill the tree
 327           ioD0toKpi=&theD0;
 328           if(isSignal) {
 329             D0TreeSgn.Fill();
 330             nD0recSgn++;
 331           } else {
 332             D0TreeBkg.Fill();
 333             if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
 334                TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
 335               nD0recBkgS++;
 336             } else {
 337               nD0recBkg++;
 338             }
 339           }
 340
 341           nD0rec1ev++;
 342
 343           ioD0toKpi=0;
 344         } // end if select
 345
 346         negtrack = 0;
 347       } // loop on negative tracks
 348       postrack = 0;
 349     }   // loop on positive tracks
 350
 351     delete [] itsEntryP;
 352     delete [] itsEntryN;
 353     delete [] brWgt;
 354     delete itsTree;
 355     delete refTree;
 356
 357     printf("\n+++\n+++ Number of D0 candidates: %d\n+++\n",nD0rec1ev);
 358     printf("----------------------------------------------------------\n");
 359
 360   }    // loop on events in file
 361
 362
 363   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
 364   printf("\n+++\n+++ Total number of D0 candidates:\n
 365                  +++    Sgn:   %d\n
 366                  +++    BkgS:  %d\n
 367                  +++    Bkg:   %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
 368
 369   delete vertexer2;
 370
 371   itstrks->Close();
 372   itsrefs->Close();
 373
 374   // store trees in files
 375   TFile* outrootSgn = new TFile(outNameSgn.Data(),"recreate");
 376   D0TreeSgn.Write();
 377   outrootSgn->Close();
 378   delete outrootSgn;
 379   TFile* outrootBkg = new TFile(outNameBkg.Data(),"recreate");
 380   D0TreeBkg.Write();
 381   outrootBkg->Close();
 382   delete outrootBkg;
 383
 384   gBenchmark->Stop(name);
 385   gBenchmark->Show(name);
 386
 387   return;
 388 }
 389 //___________________________________________________________________________
 390 Bool_t   GetInputFiles() {
 391   TString itsName("AliITStracksV2.root");
 392   TString refName("ITStracksRefFile.root");
 393
 394   if(gSystem->AccessPathName(itsName.Data(),kFileExists) ||
 395      gSystem->AccessPathName(refName.Data(),kFileExists)) return kFALSE;
 396
 397   return kTRUE;
 398 }
 399 //___________________________________________________________________________
 400 Bool_t SelectInvMass(const Double_t p[6]) {
 401
 402   Double_t mD0 = 1.8645;
 403   Double_t mPi = 0.13957;
 404   Double_t mKa = 0.49368;
 405
 406   Double_t energy[2];
 407   Double_t mom2[2],momTot2;
 408
 409   mom2[0] = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 410   mom2[1] = p[3]*p[3]+p[4]*p[4]+p[5]*p[5];
 411
 412   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
 413             (p[1]+p[4])*(p[1]+p[4])+
 414             (p[2]+p[5])*(p[2]+p[5]);
 415
 416   // D0 -> K- Pi+
 417   energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
 418   energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
 419
 420   Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
 421
 422
 423   // D0bar -> K+ Pi-
 424   energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
 425   energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
 426
 427   Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
 428
 429   if(TMath::Abs(minvD0-mD0) < kMassCut)    return kTRUE;
 430   if(TMath::Abs(minvD0bar-mD0) < kMassCut) return kTRUE;
 431   return kFALSE;
 432 }
 433 //___________________________________________________________________________
 434 void   SelectTracks(TTree& itsTree,
 435                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 436                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
 437 //
 438 // this function creates two TObjArrays with positive and negative tracks
 439 //
 440   nTrksP=0,nTrksN=0;
 441
 442
 443   Int_t entr = (Int_t)itsTree.GetEntries();
 444
 445   // trasfer tracks from tree to arrays
 446   for(Int_t i=0; i<entr; i++) {
 447
 448     AliITStrackV2 *itstrack = new AliITStrackV2;
 449     itsTree.SetBranchAddress("tracks",&itstrack);
 450
 451     itsTree.GetEvent(i);
 452
 453     // single track selection
 454     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
 455
 456     if(itstrack->Get1Pt()>0.) { // negative track
 457       trksN.AddLast(itstrack);
 458       itsEntryN[nTrksN] = i;
 459       nTrksN++;
 460     } else {                    // positive track
 461       trksP.AddLast(itstrack);
 462       itsEntryP[nTrksP] = i;
 463       nTrksP++;
 464     }
 465
 466   }
 467
 468   return;
 469 }
 470 //____________________________________________________________________________
 471 Bool_t TrkCuts(const AliITStrackV2& trk) {
 472 //
 473 // this function tells if track passes some kinematical cuts
 474 //
 475   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
 476   if(TMath::Abs(10000.*trk.GetD(gv1[0],gv1[1])) < kd0Cut) return kFALSE;
 477
 478   return kTRUE;
 479 }
 480 //____________________________________________________________________________
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491