PWG3/AliBtoJPSItoEleAnalysis.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 BtoJPSItoEle 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 AliBtoJPSItoEleTest.C
  21 //            Origin: G.E. Bruno    giuseppe.bruno@ba.infn.it
  22 //  based on Class for charm golden channel (D0->Kpi)
  23 //----------------------------------------------------------------------------
  24 #include <TKey.h>
  25 #include <TFile.h>
  26 #include <TTree.h>
  27 #include <TString.h>
  28 #include <TSystem.h>
  29 #include <TParticle.h>
  30 #include "AliESD.h"
  31 #include "AliMC.h"
  32 #include "AliRun.h"
  33 #include "AliRunLoader.h"
  34 #include "AliVertexerTracks.h"
  35 #include "AliESDVertex.h"
  36 #include "AliESDv0.h"
  37 #include "AliBtoJPSItoEle.h"
  38 #include "AliBtoJPSItoEleAnalysis.h"
  39 #include "AliLog.h"
  40
  41 typedef struct {
  42   Int_t lab;
  43   Int_t pdg;
  44   Int_t mumlab;
  45   Int_t mumpdg;
  46   Int_t gmumlab;
  47   Int_t gmumpdg;
  48   Int_t mumprongs;
  49   Float_t Vx,Vy,Vz;
  50   Float_t Px,Py,Pz;
  51 } REFTRACK;
  52
  53 ClassImp(AliBtoJPSItoEleAnalysis)
  54
  55 //----------------------------------------------------------------------------
  56 AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis():
  57 fVertexOnTheFly(kFALSE),
  58 fSim(kFALSE),
  59 fOnlySignal(kFALSE),
  60 fOnlyPrimaryJpsi(kFALSE),
  61 fPID("TRDTPCparam"),
  62 fPtCut(0.),
  63 fd0Cut(0.),
  64 fMassCut(1000.),
  65 fPidCut(0.)
  66 {
  67   // Default constructor
  68
  69   SetBCuts();
  70   SetVertex1();
  71 }
  72 //----------------------------------------------------------------------------
  73 AliBtoJPSItoEleAnalysis::~AliBtoJPSItoEleAnalysis() {}
  74 //----------------------------------------------------------------------------
  75 void AliBtoJPSItoEleAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
  76   // select candidates that pass fBCuts and write them to a new file
  77
  78   TFile *inFile = TFile::Open(inName);
  79
  80   TTree *treeBin=(TTree*)inFile->Get("TreeB");
  81   AliBtoJPSItoEleAnalysis *inAnalysis = (AliBtoJPSItoEleAnalysis*)inFile->Get("BtoJPSItoEleAnalysis");
  82   printf("+++\n+++  I N P U T   S T A T U S:\n+++\n");
  83   inAnalysis->PrintStatus();
  84
  85
  86   AliBtoJPSItoEle *d = 0;
  87   treeBin->SetBranchAddress("BtoJPSItoEle",&d);
  88   Int_t entries = (Int_t)treeBin->GetEntries();
  89
  90   printf("+++\n+++ Number of B candidates in input tree:  %d\n+++\n",entries);
  91
  92   TTree *treeBout = new TTree("TreeB","Tree with selected B candidates");
  93   treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&d,200000,0);
  94
  95
  96   Int_t okB=0;
  97   Int_t nSel = 0;
  98
  99   for(Int_t i=0; i<entries; i++) {
 100     // get event from tree
 101     treeBin->GetEvent(i);
 102
 103     if(fSim && fOnlySignal && !d->IsSignal()) continue;
 104
 105     // check if candidate passes selection (as B or Bbar)
 106     if(d->Select(fBCuts,okB)) {
 107       nSel++;
 108       treeBout->Fill();
 109     }
 110
 111   }
 112
 113   AliBtoJPSItoEleAnalysis *outAnalysis = (AliBtoJPSItoEleAnalysis*)inAnalysis->Clone("BtoJPSItoEleAnalysis");
 114   outAnalysis->SetBCuts(fBCuts);
 115   printf("------------------------------------------\n");
 116   printf("+++\n+++  O U T P U T   S T A T U S:\n+++\n");
 117   outAnalysis->PrintStatus();
 118
 119   printf("+++\n+++ Number of B mesons in output tree:  %d\n+++\n",nSel);
 120
 121   TFile* outFile = new TFile(outName,"recreate");
 122   treeBout->Write();
 123   outAnalysis->Write();
 124   outFile->Close();
 125
 126   return;
 127 }
 128 //----------------------------------------------------------------------------
 129 Double_t AliBtoJPSItoEleAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
 130                                               Double_t time) const {
 131   // calculated the mass from momentum, track length from vertex to TOF
 132   // and time measured by the TOF
 133   if(length==0.) return -1000.;
 134   Double_t a = time*time/length/length;
 135   if(a > 1.) {
 136     a = TMath::Sqrt(a-1.);
 137   } else {
 138     a = -TMath::Sqrt(1.-a);
 139   }
 140
 141   return mom*a;
 142 }
 143 //----------------------------------------------------------------------------
 144 void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
 145                                         const Char_t *outName) {
 146   // Find candidates and calculate parameters
 147
 148
 149   TString esdName="AliESDs.root";
 150   if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
 151     printf("AliBtoJPSItoEleAnalysis::FindCandidatesESD(): No ESDs file found!\n");
 152     return;
 153   }
 154
 155   TString outName1=outName;
 156   TString outName2="nTotEvents.dat";
 157
 158   Int_t    nTotEv=0,nBrec=0,nBrec1ev=0;
 159   Double_t dca;
 160   Double_t v2[3],mom[6],d0[2];
 161   Int_t    iTrkP,iTrkN,trkEntries;
 162   Int_t    nTrksP=0,nTrksN=0;
 163   Int_t    trkNum[2];
 164   Double_t tofmass[2];
 165   Int_t    okB=0;
 166   AliESDtrack *postrack = 0;
 167   AliESDtrack *negtrack = 0;
 168
 169   // create the AliVertexerTracks object
 170   // (it will be used only if fVertexOnTheFly=kTrue)
 171   AliVertexerTracks *vertexer1 = new AliVertexerTracks();
 172   if(fVertexOnTheFly) {
 173     // open the mean vertex
 174     TFile *invtx = new TFile("AliESDVertexMean.root");
 175     AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
 176     invtx->Close();
 177     vertexer1->SetVtxStart(initVertex);
 178     delete invtx;
 179   }
 180   Int_t  skipped[2];
 181   Bool_t goodVtx1;
 182
 183   // create tree for reconstructed decayes
 184   AliBtoJPSItoEle *ioBtoJPSItoEle=0;
 185   TTree *treeB = new TTree("TreeB","Tree with candidates");
 186   treeB->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&ioBtoJPSItoEle,200000,0);
 187
 188   // open file with tracks
 189   TFile *esdFile = TFile::Open(esdName.Data());
 190   AliESD* event = new AliESD;
 191   TTree* tree = (TTree*) esdFile->Get("esdTree");
 192   if(!tree) {
 193     Error("FindCandidatesESD", "no ESD tree found");
 194     return;
 195   }
 196   tree->SetBranchAddress("ESD",&event);
 197
 198   // loop on events in file
 199   for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
 200     if(iEvent > evLast) break;
 201     tree->GetEvent(iEvent);
 202     Int_t ev = (Int_t)event->GetEventNumberInFile();
 203     printf("--- Finding B -> JPSI -> e+ e- in event %d\n",ev);
 204     // count the total number of events
 205     nTotEv++;
 206
 207     trkEntries = (Int_t)event->GetNumberOfTracks();
 208     printf(" Number of tracks: %d\n",trkEntries);
 209     if(trkEntries<2) continue;
 210
 211     // retrieve primary vertex from file
 212     if(!event->GetPrimaryVertex()) {
 213       printf(" No vertex\n");
 214       continue;
 215     }
 216     event->GetPrimaryVertex()->GetXYZ(fV1);
 217
 218     // call function which applies sigle-track selection and
 219     // separetes positives and negatives
 220     TObjArray trksP(trkEntries/2);
 221     Int_t    *trkEntryP   = new Int_t[trkEntries];
 222     TObjArray trksN(trkEntries/2);
 223     Int_t    *trkEntryN   = new Int_t[trkEntries];
 224     TTree *trkTree = new TTree();
 225     SelectTracks(event,trksP,trkEntryP,nTrksP,
 226                        trksN,trkEntryN,nTrksN);
 227
 228     printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
 229
 230
 231     nBrec1ev = 0;
 232
 233     // LOOP ON  POSITIVE  TRACKS
 234     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 235       if(iTrkP%1==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
 236
 237       // get track from track array
 238       postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
 239       trkNum[0] = trkEntryP[iTrkP];
 240
 241       // LOOP ON  NEGATIVE  TRACKS
 242       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 243
 244         // get track from tracks array
 245         negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
 246         trkNum[1] = trkEntryN[iTrkN];
 247
 248         //
 249         // ----------- DCA MINIMIZATION ------------------
 250         //
 251         // find the DCA and propagate the tracks to the DCA
 252         Double_t b=event->GetMagneticField();
 253         AliESDtrack nt(*negtrack), pt(*postrack);
 254         dca = nt.PropagateToDCA(&pt,b);
 255
 256         // define the AliESDv0 object
 257         AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
 258
 259         // get position of the secondary vertex
 260         vertex2.GetXYZ(v2[0],v2[1],v2[2]);
 261         vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
 262         vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
 263         // impact parameters of the tracks w.r.t. the primary vertex
 264         // d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
 265         // d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
 266         goodVtx1 = kTRUE;
 267
 268         // no vertexing if DeltaMass > fMassCut
 269         if(fVertexOnTheFly) {
 270           goodVtx1 = kFALSE;
 271           if(SelectInvMass(mom)) {
 272             // primary vertex from *other* tracks in the event
 273             vertexer1->SetFieldkG(event->GetMagneticField());
 274             skipped[0] = trkEntryP[iTrkP];
 275             skipped[1] = trkEntryN[iTrkN];
 276             vertexer1->SetSkipTracks(2,skipped);
 277             AliESDVertex *vertex1onfly =
 278               (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
 279             if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
 280             vertex1onfly->GetXYZ(fV1);
 281             //vertex1onfly->PrintStatus();
 282             delete vertex1onfly;
 283           }
 284         }
 285
 286         // impact parameters of the tracks w.r.t. the primary vertex
 287         d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
 288         d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
 289
 290         // create the object AliBtoJPSItoEle
 291         AliBtoJPSItoEle theB(ev,trkNum,fV1,v2,dca,mom,d0);
 292         // select B's
 293         if(goodVtx1 && theB.Select(fBCuts,okB)) {
 294           // get PID info from ESD
 295           AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
 296           Double_t esdpid0[5];
 297           t0->GetESDpid(esdpid0);
 298           if(t0->GetStatus()&AliESDtrack::kTOFpid) {
 299             tofmass[0] = CalculateTOFmass(t0->GetP(),
 300                                           t0->GetIntegratedLength(),
 301                                           t0->GetTOFsignal());
 302           } else {
 303             tofmass[0] = -1000.;
 304           }
 305           AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
 306           Double_t esdpid1[5];
 307           t1->GetESDpid(esdpid1);
 308           if(t1->GetStatus()&AliESDtrack::kTOFpid) {
 309             tofmass[1] = CalculateTOFmass(t1->GetP(),
 310                                           t1->GetIntegratedLength(),
 311                                           t1->GetTOFsignal());
 312           } else {
 313             tofmass[1] = -1000.;
 314           }
 315
 316           theB.SetPIDresponse(esdpid0,esdpid1);
 317           theB.SetTOFmasses(tofmass);
 318
 319           // fill the tree
 320           ioBtoJPSItoEle=&theB;
 321           treeB->Fill();
 322
 323           nBrec++; nBrec1ev++;
 324           ioBtoJPSItoEle=0;
 325         }
 326
 327         negtrack = 0;
 328       } // loop on negative tracks
 329       postrack = 0;
 330     }   // loop on positive tracks
 331
 332     delete [] trkEntryP;
 333     delete [] trkEntryN;
 334     delete trkTree;
 335
 336     printf(" Number of B candidates: %d\n",nBrec1ev);
 337   }    // loop on events in file
 338
 339
 340   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
 341   printf("\n+++\n+++ Total number of B candidates: %d\n+++\n",nBrec);
 342
 343   delete vertexer1;
 344
 345   esdFile->Close();
 346
 347   // create a copy of this class to be written to output file
 348   AliBtoJPSItoEleAnalysis *copy = (AliBtoJPSItoEleAnalysis*)this->Clone("BtoJPSItoEleAnalysis");
 349
 350   // add PDG codes to decay tracks in found candidates (in simulation mode)
 351   // and store tree in the output file
 352   if(!fSim) {
 353     TFile *outroot = new TFile(outName1.Data(),"recreate");
 354     treeB->Write();
 355     copy->Write();
 356     outroot->Close();
 357     delete outroot;
 358   } else {
 359     printf(" Now adding information from simulation (PDG codes) ...\n");
 360     TTree *treeBsim = new TTree("TreeB","Tree with B candidates");
 361     SimulationInfo(treeB,treeBsim);
 362     delete treeB;
 363     TFile *outroot = new TFile(outName1.Data(),"recreate");
 364     treeBsim->Write();
 365     copy->Write();
 366     outroot->Close();
 367     delete outroot;
 368   }
 369
 370   // write to a file the total number of events
 371   FILE *outdat = fopen(outName2.Data(),"w");
 372   fprintf(outdat,"%d\n",nTotEv);
 373   fclose(outdat);
 374
 375   return;
 376 }
 377 //-----------------------------------------------------------------------------
 378 void AliBtoJPSItoEleAnalysis::PrintStatus() const {
 379   // Print parameters being used
 380
 381   printf("Preselections:\n");
 382   printf("    fPtCut   = %f GeV\n",fPtCut);
 383   printf("    fd0Cut   = %f micron\n",fd0Cut);
 384   printf("    fMassCut = %f GeV\n",fMassCut);
 385   printf("    fPidCut  > %f \n",fPidCut);
 386   if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
 387   if(fSim) {
 388     printf("Simulation mode\n");
 389     if(fOnlySignal && !(fOnlyPrimaryJpsi)) printf("  Only signal goes to file\n");
 390     if(fOnlyPrimaryJpsi && !(fOnlySignal)) printf("  Only primary Jpsi go to file\n");
 391     if(fOnlyPrimaryJpsi && fOnlySignal) printf("  Both signal and primary Jpsi go to file\n");
 392   }
 393   printf("Cuts on candidates:\n");
 394   printf("    |M-MJPsi| [GeV]  < %f\n",fBCuts[0]);
 395   printf("    dca    [micron]  < %f\n",fBCuts[1]);
 396   printf("    cosThetaStar     < %f\n",fBCuts[2]);
 397   printf("    pTP     [GeV]    > %f\n",fBCuts[3]);
 398   printf("    pTN     [GeV]    > %f\n",fBCuts[4]);
 399   printf("    |d0P|  [micron]  < %f\n",fBCuts[5]);
 400   printf("    |d0N|  [micron]  < %f\n",fBCuts[6]);
 401   printf("    d0d0  [micron^2] < %f\n",fBCuts[7]);
 402   printf("    cosThetaPoint    > %f\n",fBCuts[8]);
 403
 404   return;
 405 }
 406 //-----------------------------------------------------------------------------
 407 Bool_t AliBtoJPSItoEleAnalysis::SelectInvMass(const Double_t p[6]) const {
 408   // Apply preselection in the invariant mass of the pair
 409
 410   Double_t mJPsi = 3.096916;
 411   Double_t mel   = 0.00510998902;
 412
 413   Double_t energy[2];
 414   Double_t mom2[2],momTot2;
 415
 416   mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
 417   mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
 418
 419   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
 420             (p[1]+p[4])*(p[1]+p[4])+
 421             (p[2]+p[5])*(p[2]+p[5]);
 422
 423   // J/Psi -> e+ e-
 424   energy[1] = TMath::Sqrt(mel*mel+mom2[1]);
 425   energy[0] = TMath::Sqrt(mel*mel+mom2[0]);
 426
 427   Double_t minvJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
 428
 429   if(TMath::Abs(minvJPsi-mJPsi)  < fMassCut) return kTRUE;
 430   return kFALSE;
 431 }
 432 //-----------------------------------------------------------------------------
 433 void AliBtoJPSItoEleAnalysis::SelectTracks(AliESD *event,
 434         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
 435         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
 436   // Create two TObjArrays with positive and negative tracks and
 437   // apply single-track preselection
 438
 439   nTrksP=0,nTrksN=0;
 440
 441   Int_t entr = event->GetNumberOfTracks();
 442
 443   // transfer ITS tracks from ESD to arrays and to a tree
 444   for(Int_t i=0; i<entr; i++) {
 445
 446     AliESDtrack *esdtrack = event->GetTrack(i);
 447     UInt_t status = esdtrack->GetStatus();
 448
 449     if(!(status&AliESDtrack::kITSin)) continue;
 450
 451     // single track selection
 452     if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
 453
 454     if(esdtrack->GetSign()<0) { // negative track
 455       trksN.AddLast(esdtrack);
 456       trkEntryN[nTrksN] = i;
 457       nTrksN++;
 458     } else {                 // positive track
 459       trksP.AddLast(esdtrack);
 460       trkEntryP[nTrksP] = i;
 461       nTrksP++;
 462     }
 463
 464   } // loop on ESD tracks
 465
 466   return;
 467 }
 468 //-----------------------------------------------------------------------------
 469 void AliBtoJPSItoEleAnalysis::SetBCuts(Double_t cut0,Double_t cut1,
 470                                    Double_t cut2,Double_t cut3,Double_t cut4,
 471                                    Double_t cut5,Double_t cut6,
 472                                    Double_t cut7,Double_t cut8) {
 473   // Set the cuts for B selection
 474   fBCuts[0] = cut0;
 475   fBCuts[1] = cut1;
 476   fBCuts[2] = cut2;
 477   fBCuts[3] = cut3;
 478   fBCuts[4] = cut4;
 479   fBCuts[5] = cut5;
 480   fBCuts[6] = cut6;
 481   fBCuts[7] = cut7;
 482   fBCuts[8] = cut8;
 483
 484   return;
 485 }
 486 //-----------------------------------------------------------------------------
 487 void AliBtoJPSItoEleAnalysis::SetBCuts(const Double_t cuts[9]) {
 488   // Set the cuts for B selection
 489
 490   for(Int_t i=0; i<9; i++) fBCuts[i] = cuts[i];
 491
 492   return;
 493 }
 494 //-----------------------------------------------------------------------------
 495 Bool_t
 496 AliBtoJPSItoEleAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
 497   // Check if track passes some kinematical cuts
 498   // Magnetic field "b" (kG)
 499
 500   if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
 501     return kFALSE;
 502   if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
 503     return kFALSE;
 504   //select only tracks with the "combined PID"
 505   UInt_t status = trk.GetStatus();
 506   if ((status&AliESDtrack::kESDpid)==0) return kTRUE;
 507   Double_t r[5];
 508   trk.GetESDpid(r);
 509   if(r[0] < fPidCut) return kFALSE;
 510
 511   return kTRUE;
 512 }
 513 //----------------------------------------------------------------------------
 514 void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *gAlice,
 515                                            Int_t evFirst,Int_t evLast) const {
 516   // Create a file with simulation info for the reconstructed tracks
 517
 518   TFile *outFile = TFile::Open("BTracksRefFile.root","recreate");
 519   TFile *esdFile = TFile::Open("AliESDs.root");
 520
 521   AliMC *mc = gAlice->GetMCApp();
 522
 523   Int_t      label;
 524   TParticle *part;
 525   TParticle *mumpart;
 526   TParticle *gmumpart;
 527   REFTRACK   reftrk;
 528
 529   AliESD* event = new AliESD;
 530   TTree* tree = (TTree*) esdFile->Get("esdTree");
 531   tree->SetBranchAddress("ESD",&event);
 532   // loop on events in file
 533   for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
 534     if(iEvent>evLast) break;
 535     tree->GetEvent(iEvent);
 536     Int_t ev = (Int_t)event->GetEventNumberInFile();
 537
 538     gAlice->GetEvent(ev);
 539
 540     Int_t nentr=(Int_t)event->GetNumberOfTracks();
 541
 542     // Tree for true track parameters
 543     char ttname[100];
 544     sprintf(ttname,"Tree_Ref_%d",ev);
 545     TTree *reftree = new TTree(ttname,"Tree with true track params");
 546     reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:gmumlab:gmumpdg:mumprongs:Vx/F:Vy:Vz:Px:Py:Pz");
 547 //    reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
 548
 549     for(Int_t i=0; i<nentr; i++) {
 550       AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
 551       label = TMath::Abs(esdtrack->GetLabel());
 552
 553       part = (TParticle*)mc->Particle(label);
 554       reftrk.lab = label;
 555       reftrk.pdg = part->GetPdgCode();
 556       reftrk.mumlab = part->GetFirstMother();
 557       if(part->GetFirstMother()>=0) {
 558         mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
 559         reftrk.mumpdg = mumpart->GetPdgCode();
 560         reftrk.mumprongs = mumpart->GetNDaughters();
 561         reftrk.gmumlab = mumpart->GetFirstMother();
 562         if(mumpart->GetFirstMother()>=0) {
 563           gmumpart = (TParticle*)gAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
 564           reftrk.gmumpdg = gmumpart->GetPdgCode();
 565         }
 566       } else {
 567         reftrk.mumpdg=-1;
 568         reftrk.mumprongs=-1;
 569         reftrk.gmumpdg=-1;
 570         reftrk.gmumlab=part->GetFirstMother(); // If it hasn't any mother, then it has neither Gmother!
 571         // reftrk.gmumlab=-1; // If it hasn't any mother, then it has neither Gmother!
 572       }
 573       reftrk.Vx = part->Vx();
 574       reftrk.Vy = part->Vy();
 575       reftrk.Vz = part->Vz();
 576       reftrk.Px = part->Px();
 577       reftrk.Py = part->Py();
 578       reftrk.Pz = part->Pz();
 579
 580       reftree->Fill();
 581
 582     } // loop on tracks
 583
 584     outFile->cd();
 585     reftree->Write();
 586
 587     delete reftree;
 588   } // loop on events
 589
 590   esdFile->Close();
 591   outFile->Close();
 592
 593   return;
 594 }
 595 //-----------------------------------------------------------------------------
 596 void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) const {
 597   // add pdg codes to candidate decay tracks (for sim)
 598
 599   TString refFileName("BTracksRefFile.root");
 600   if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
 601     printf("AliBtoJPSItoEleAnalysis::SimulationInfo: no reference file found!\n");
 602     return;
 603   }
 604   TFile *refFile = TFile::Open(refFileName.Data());
 605
 606   Char_t refTreeName[100];
 607   Int_t  event;
 608   Int_t  pdg[2],mumpdg[2],mumlab[2],gmumpdg[2],gmumlab[2];
 609   REFTRACK reftrk;
 610
 611   // read-in reference tree for event 0 (the only event for Pb-Pb)
 612   sprintf(refTreeName,"Tree_Ref_%d",0);
 613   TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
 614   refTree0->SetBranchAddress("rectracks",&reftrk);
 615
 616   AliBtoJPSItoEle *theB = 0;
 617   treeBin->SetBranchAddress("BtoJPSItoEle",&theB);
 618   treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&theB,200000,0);
 619
 620   Int_t entries = (Int_t)treeBin->GetEntries();
 621
 622   for(Int_t i=0; i<entries; i++) {
 623     if(i%100==0) printf("  done %d candidates of %d\n",i,entries);
 624
 625     treeBin->GetEvent(i);
 626     event = theB->EventNo();
 627
 628     if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
 629       refTree0->GetEvent(theB->GetTrkNum(0));
 630       pdg[0]    = reftrk.pdg;
 631       mumpdg[0] = reftrk.mumpdg;
 632       mumlab[0] = reftrk.mumlab;
 633       gmumpdg[0] = reftrk.gmumpdg;
 634       gmumlab[0] = reftrk.gmumlab;
 635       refTree0->GetEvent(theB->GetTrkNum(1));
 636       pdg[1]    = reftrk.pdg;
 637       mumpdg[1] = reftrk.mumpdg;
 638       mumlab[1] = reftrk.mumlab;
 639       gmumpdg[1] = reftrk.gmumpdg;
 640       gmumlab[1] = reftrk.gmumlab;
 641     } else {
 642       sprintf(refTreeName,"Tree_Ref_%d",event);
 643       TTree *refTree = (TTree*)refFile->Get(refTreeName);
 644       refTree->SetBranchAddress("rectracks",&reftrk);
 645       refTree->GetEvent(theB->GetTrkNum(0));
 646       pdg[0]    = reftrk.pdg;
 647       mumpdg[0] = reftrk.mumpdg;
 648       mumlab[0] = reftrk.mumlab;
 649       gmumpdg[0] = reftrk.gmumpdg;
 650       gmumlab[0] = reftrk.gmumlab;
 651       refTree->GetEvent(theB->GetTrkNum(1));
 652       pdg[1]    = reftrk.pdg;
 653       mumpdg[1] = reftrk.mumpdg;
 654       mumlab[1] = reftrk.mumlab;
 655       gmumpdg[1] = reftrk.gmumpdg;
 656       gmumlab[1] = reftrk.gmumlab;
 657       delete refTree;
 658     }
 659
 660     theB->SetPdgCodes(pdg);
 661     theB->SetMumPdgCodes(mumpdg);
 662     theB->SetGMumPdgCodes(gmumpdg);
 663
 664     if(gmumpdg[0]==gmumpdg[1] &&              // Both GrandMothers are of the same sign
 665        (TMath::Abs(gmumpdg[0])==521   || TMath::Abs(gmumpdg[0])==511   ||  // GrandMother Bplus/Bminus or B0/B0bar
 666         TMath::Abs(gmumpdg[0])==523   || TMath::Abs(gmumpdg[0])==513   ||  // B0s/B0sbar
 667         TMath::Abs(gmumpdg[0])==515   || TMath::Abs(gmumpdg[0])==525   ||  //
 668         TMath::Abs(gmumpdg[0])==531   || TMath::Abs(gmumpdg[0])==533   ||  //
 669         TMath::Abs(gmumpdg[0])==535   || TMath::Abs(gmumpdg[0])==541   ||  //
 670         TMath::Abs(gmumpdg[0])==543   || TMath::Abs(gmumpdg[0])==545   ||  //
 671         TMath::Abs(gmumpdg[0])==10521 || TMath::Abs(gmumpdg[0])==10511 ||  //   all possible
 672         TMath::Abs(gmumpdg[0])==10523 || TMath::Abs(gmumpdg[0])==10513 ||  //    B mesons
 673         TMath::Abs(gmumpdg[0])==20523 || TMath::Abs(gmumpdg[0])==20513 ||  //
 674         TMath::Abs(gmumpdg[0])==10531 || TMath::Abs(gmumpdg[0])==10533 ||  //
 675         TMath::Abs(gmumpdg[0])==20533 || TMath::Abs(gmumpdg[0])==10541 ||  //
 676         TMath::Abs(gmumpdg[0])==20543 || TMath::Abs(gmumpdg[0])==10543 ||  //
 677         TMath::Abs(gmumpdg[0])==4122  || TMath::Abs(gmumpdg[0])==4222  ||  // All possible B baryons
 678         TMath::Abs(gmumpdg[0])==4212  || TMath::Abs(gmumpdg[0])==4112  ||  // All possible B baryons
 679         TMath::Abs(gmumpdg[0])==4224  || TMath::Abs(gmumpdg[0])==4214  ||  // All possible B baryons
 680         TMath::Abs(gmumpdg[0])==4114  || TMath::Abs(gmumpdg[0])==4232  ||  // All possible B baryons
 681         TMath::Abs(gmumpdg[0])==4132  || TMath::Abs(gmumpdg[0])==4322  ||  // All possible B baryons
 682         TMath::Abs(gmumpdg[0])==4312  || TMath::Abs(gmumpdg[0])==4324  ||  // All possible B baryons
 683         TMath::Abs(gmumpdg[0])==4314  || TMath::Abs(gmumpdg[0])==4332  ||  // All possible B baryons
 684         TMath::Abs(gmumpdg[0])==4334  || TMath::Abs(gmumpdg[0])==4412  ||  // All possible B baryons
 685         TMath::Abs(gmumpdg[0])==4422  || TMath::Abs(gmumpdg[0])==4414  ||  // All possible B baryons
 686         TMath::Abs(gmumpdg[0])==4424  || TMath::Abs(gmumpdg[0])==4432  ||  // All possible B baryons
 687         TMath::Abs(gmumpdg[0])==4434  || TMath::Abs(gmumpdg[0])==4444      // All possible B baryons
 688        ) &&
 689        mumpdg[0]==443 && mumpdg[1]== 443 &&
 690        mumlab[0]==mumlab[1] &&
 691        reftrk.mumprongs==2 &&
 692        pdg[0]==-11 && pdg[1]==11
 693      ) theB->SetSignal();
 694
 695     else if (  // here consider the case of primary J/psi
 696       gmumlab[0]<0 && gmumlab[1]<0 &&
 697        mumpdg[0]==443 && mumpdg[1]== 443 &&
 698        mumlab[0]==mumlab[1] &&
 699        reftrk.mumprongs==2 &&
 700        pdg[0]==-11 && pdg[1]==11
 701      ) theB->SetJpsiPrimary();
 702
 703     // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();
 704
 705  // write it out 1) always if you have not asked for onlySignal or OnlyPrimaryJpsi (or both)
 706  //          or  2) if you have asked for Signal and it is Signal
 707  //          or  3) if you have asked for Primary Jpsi and it is a Primary Jpsi
 708     if ( (!fOnlySignal && !fOnlyPrimaryJpsi) || (fOnlySignal && theB->IsSignal())
 709         || (fOnlyPrimaryJpsi && theB->IsJpsiPrimary()) ) treeBout->Fill();
 710   }
 711
 712   delete refTree0;
 713
 714   refFile->Close();
 715
 716   return;
 717 }
 718
 719
 720
 721
 722
 723