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