PWG3/vertexingOld/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 "AliESDEvent.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 #include "AliKFParticleBase.h"
  41 #include "AliKFParticle.h"
  42 #include "AliKFVertex.h"
  43
  44 typedef struct {
  45   Int_t lab;
  46   Int_t pdg;
  47   Int_t mumlab;
  48   Int_t mumpdg;
  49   Int_t gmumlab;
  50   Int_t gmumpdg;
  51   Int_t mumprongs;
  52   Float_t Vx,Vy,Vz;
  53   Float_t Px,Py,Pz;
  54 } REFTRACK;
  55
  56 ClassImp(AliBtoJPSItoEleAnalysis)
  57
  58 //----------------------------------------------------------------------------
  59 AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis():
  60 fVertexOnTheFly(kFALSE),
  61 fSim(kFALSE),
  62 fOnlySignal(kFALSE),
  63 fOnlyPrimaryJpsi(kFALSE),
  64 fPID("TRDTPCparam"),
  65 fPtCut(0.),
  66 fd0Cut(0.),
  67 fMassCut(1000.),
  68 fPidCut(0.),
  69 //fKFPrimVertex(kFALSE),
  70 //fKFTopConstr(kFALSE),
  71 fKFSecondVertex(kFALSE)
  72 {
  73   // Default constructor
  74
  75   SetBCuts();
  76   SetVertex1();
  77 }
  78 //----------------------------------------------------------------------------
  79 AliBtoJPSItoEleAnalysis::~AliBtoJPSItoEleAnalysis() {}
  80 //----------------------------------------------------------------------------
  81 void AliBtoJPSItoEleAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
  82   // select candidates that pass fBCuts and write them to a new file
  83
  84   TFile *inFile = TFile::Open(inName);
  85
  86   TTree *treeBin=(TTree*)inFile->Get("TreeB");
  87   AliBtoJPSItoEleAnalysis *inAnalysis = (AliBtoJPSItoEleAnalysis*)inFile->Get("BtoJPSItoEleAnalysis");
  88   printf("+++\n+++  I N P U T   S T A T U S:\n+++\n");
  89   inAnalysis->PrintStatus();
  90
  91
  92   AliBtoJPSItoEle *d = 0;
  93   treeBin->SetBranchAddress("BtoJPSItoEle",&d);
  94   Int_t entries = (Int_t)treeBin->GetEntries();
  95
  96   printf("+++\n+++ Number of B candidates in input tree:  %d\n+++\n",entries);
  97
  98   TTree *treeBout = new TTree("TreeB","Tree with selected B candidates");
  99   treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&d,200000,0);
 100
 101
 102   Int_t okB=0;
 103   Int_t nSel = 0;
 104
 105   for(Int_t i=0; i<entries; i++) {
 106     // get event from tree
 107     treeBin->GetEvent(i);
 108
 109     if(fSim && fOnlySignal && !d->IsSignal()) continue;
 110
 111     // check if candidate passes selection (as B or Bbar)
 112     if(d->Select(fBCuts,okB)) {
 113       nSel++;
 114       treeBout->Fill();
 115     }
 116
 117   }
 118
 119   AliBtoJPSItoEleAnalysis *outAnalysis = (AliBtoJPSItoEleAnalysis*)inAnalysis->Clone("BtoJPSItoEleAnalysis");
 120   outAnalysis->SetBCuts(fBCuts);
 121   printf("------------------------------------------\n");
 122   printf("+++\n+++  O U T P U T   S T A T U S:\n+++\n");
 123   outAnalysis->PrintStatus();
 124
 125   printf("+++\n+++ Number of B mesons in output tree:  %d\n+++\n",nSel);
 126
 127   TFile* outFile = new TFile(outName,"recreate");
 128   treeBout->Write();
 129   outAnalysis->Write();
 130   outFile->Close();
 131
 132   return;
 133 }
 134 //----------------------------------------------------------------------------
 135 Double_t AliBtoJPSItoEleAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
 136                                               Double_t time) const {
 137   // calculated the mass from momentum, track length from vertex to TOF
 138   // and time measured by the TOF
 139   if(length==0.) return -1000.;
 140   Double_t a = time*time/length/length;
 141   if(a > 1.) {
 142     a = TMath::Sqrt(a-1.);
 143   } else {
 144     a = -TMath::Sqrt(1.-a);
 145   }
 146
 147   return mom*a;
 148 }
 149 //----------------------------------------------------------------------------
 150 void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
 151                                         const Char_t *outName) {
 152   // Find candidates and calculate parameters
 153
 154
 155   TString esdName="AliESDs.root";
 156   if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
 157     printf("AliBtoJPSItoEleAnalysis::FindCandidatesESD(): No ESDs file found!\n");
 158     return;
 159   }
 160
 161   TString outName1=outName;
 162   TString outName2="nTotEvents.dat";
 163
 164   Int_t    nTotEv=0,nBrec=0,nBrec1ev=0;
 165   Double_t dca;
 166   Double_t v2[3],mom[6],d0[2];
 167   Int_t    iTrkP,iTrkN,trkEntries;
 168   Int_t    nTrksP=0,nTrksN=0;
 169   Int_t    trkNum[2];
 170   Double_t tofmass[2];
 171   Int_t    okB=0;
 172   AliESDtrack *postrack = 0;
 173   AliESDtrack *negtrack = 0;
 174
 175   // create the AliVertexerTracks object
 176   // (it will be used only if fVertexOnTheFly=kTrue)
 177   AliVertexerTracks *vertexer1 = new AliVertexerTracks();
 178   if(fVertexOnTheFly) {
 179     // open the mean vertex
 180     TFile *invtx = new TFile("AliESDVertexMean.root");
 181     AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
 182     invtx->Close();
 183     vertexer1->SetVtxStart(initVertex);
 184     delete invtx;
 185   }
 186   Int_t  skipped[2];
 187   Bool_t goodVtx1;
 188
 189   // create tree for reconstructed decayes
 190   AliBtoJPSItoEle *ioBtoJPSItoEle=0;
 191   TTree *treeB = new TTree("TreeB","Tree with candidates");
 192   treeB->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&ioBtoJPSItoEle,200000,0);
 193
 194   // open file with tracks
 195   TFile *esdFile = TFile::Open(esdName.Data());
 196   AliESDEvent* event = new AliESDEvent();
 197   TTree* tree = (TTree*) esdFile->Get("esdTree");
 198   if(!tree) {
 199     Error("FindCandidatesESD", "no ESD tree found");
 200     return;
 201   }
 202   event->ReadFromTree(tree);
 203
 204 /*  if (fKFPrimVertex)
 205   AliRunLoader* runLoader = 0;
 206   {
 207     if (gAlice) {
 208       delete gAlice->GetRunLoader();
 209       delete gAlice;
 210       gAlice=0;
 211     }
 212     runLoader = AliRunLoader::Open(galName.Data());
 213     if (runLoader == 0x0) {
 214       cerr<<"Can not open session"<<endl;
 215       return;
 216     }
 217     cout << "Ok open galice.root" << endl;
 218     runLoader->LoadgAlice();
 219
 220     gAlice = runLoader->GetAliRun();
 221     runLoader->LoadKinematics();
 222     runLoader->LoadHeader();
 223   } */
 224
 225   // loop on events in file
 226   for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
 227     if(iEvent > evLast) break;
 228     tree->GetEvent(iEvent);
 229     Int_t ev = (Int_t)event->GetEventNumberInFile();
 230     printf("--- Finding B -> JPSI -> e+ e- in event %d\n",ev);
 231     // count the total number of events
 232     nTotEv++;
 233
 234     trkEntries = (Int_t)event->GetNumberOfTracks();
 235     printf(" Number of tracks: %d\n",trkEntries);
 236     if(trkEntries<2) continue;
 237
 238     // retrieve primary vertex from file
 239     if(!event->GetPrimaryVertex()) {
 240       printf(" No vertex\n");
 241       continue;
 242     }
 243     event->GetPrimaryVertex()->GetXYZ(fV1);
 244
 245     // call function which applies sigle-track selection and
 246     // separetes positives and negatives
 247     TObjArray trksP(trkEntries/2);
 248     Int_t    *trkEntryP   = new Int_t[trkEntries];
 249     TObjArray trksN(trkEntries/2);
 250     Int_t    *trkEntryN   = new Int_t[trkEntries];
 251     TTree *trkTree = new TTree();
 252     SelectTracks(event,trksP,trkEntryP,nTrksP,
 253                        trksN,trkEntryN,nTrksN);
 254
 255     printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
 256
 257
 258     nBrec1ev = 0;
 259
 260 /*
 261 //===================== PRIMARY VERTEX USING KF METHODS ==========================//
 262
 263 if (fKFPrimVertex)
 264 {
 265   AliStack* stack = runLoader->Stack();
 266
 267   class TESDTrackInfo
 268     {
 269     public:
 270     TESDTrackInfo(){}
 271     AliKFParticle fParticle; // KFParticle constructed from ESD track
 272     //Bool_t fPrimUsedFlag;    // flag says that the particle was used for primary vertex fit
 273     Bool_t fOK;              // is the track good enough
 274     Int_t mcPDG;             // Monte Carlo PDG code of the particle
 275     Int_t mcMotherID;        // Monte Carlo ID of its mother
 276     };
 277
 278  // TESDTrackInfo ESDTrackInfo[trkEntries];
 279   TESDTrackInfo ESDTrackInfo[1000];
 280   for (Int_t iTr=0; iTr<trkEntries; iTr++)
 281     {
 282       TESDTrackInfo &info = ESDTrackInfo[iTr];
 283       info.fOK = 0;
 284       //info.fPrimUsedFlag = 0;
 285       info.mcPDG = -1;
 286       info.mcMotherID = -1;
 287
 288       // track quality check
 289
 290       AliESDtrack *pTrack = event->GetTrack(iTr);
 291       if( !pTrack  ) continue;
 292       if (pTrack->GetKinkIndex(0)>0) continue;
 293       if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
 294       //Int_t indi[12];
 295       //if( pTrack->GetITSclusters(indi) <5 ) continue;
 296       //Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
 297
 298       // take MC PDG
 299
 300       Int_t mcID = TMath::Abs(pTrack->GetLabel());
 301       TParticle * part = stack->Particle(TMath::Abs(mcID));
 302       info.mcPDG = part->GetPdgCode();
 303       Int_t PDG = info.mcPDG;
 304       if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
 305
 306
 307       // Construct KFParticle for the track
 308
 309       info.fParticle = AliKFParticle( *pTrack, PDG );
 310       info.fOK = 1;
 311     }
 312
 313     // Find event primary vertex with KF methods
 314
 315    AliKFVertex primVtx;
 316     {
 317    // const AliKFParticle * vSelected[trkEntries]; // Selected particles for vertex fit
 318    // Int_t vIndex[trkEntries];                    // Indices of selected particles
 319    // Bool_t vFlag[trkEntries];                    // Flags returned by the vertex finder
 320    const AliKFParticle * vSelected[1000]; // Selected particles for vertex fit
 321    Int_t vIndex[1000];                    // Indices of selected particles
 322    Bool_t vFlag[1000];                    // Flags returned by the vertex finder
 323
 324     Int_t nSelected = 0;
 325     for( Int_t i = 0; i<trkEntries; i++){
 326       if(ESDTrackInfo[i].fOK ){
 327         vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
 328         vIndex[nSelected] = i;
 329         nSelected++;
 330       }
 331     }
 332     primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
 333     //for( Int_t i = 0; i<nSelected; i++){
 334     //  if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
 335     //}
 336     if( primVtx.GetNDF() <1 ) return; // Less then two tracks in primary vertex
 337     }
 338
 339 }*/
 340
 341     // LOOP ON  POSITIVE  TRACKS
 342     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 343       if(iTrkP%1==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
 344
 345       // get track from track array
 346       postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
 347       trkNum[0] = trkEntryP[iTrkP];
 348
 349       // LOOP ON  NEGATIVE  TRACKS
 350       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 351
 352         // get track from tracks array
 353         negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
 354         trkNum[1] = trkEntryN[iTrkN];
 355
 356         //
 357         // ----------- DCA MINIMIZATION ------------------
 358         //
 359         // find the DCA and propagate the tracks to the DCA
 360         Double_t b=event->GetMagneticField();
 361         AliESDtrack nt(*negtrack), pt(*postrack);
 362         dca = nt.PropagateToDCA(&pt,b);
 363
 364         // define the AliESDv0 object
 365         AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
 366
 367         // get position of the secondary vertex
 368
 369         if (fKFSecondVertex){
 370           //Define the AliKFParticle Objects
 371           AliKFParticle trackP = AliKFParticle(pt,-11);
 372           AliKFParticle trackN = AliKFParticle(nt,11);
 373           //Construct the V0like mother
 374           AliKFParticle V0(trackP,trackN);
 375           //Get global position of the secondary vertex using KF methods
 376           v2[0] = V0.GetX();
 377           v2[1] = V0.GetY();
 378           v2[2] = V0.GetZ();
 379           mom[0] = trackP.GetPx(); mom[1] = trackP.GetPy(); mom[2] = trackP.GetPz();
 380           mom[3] = trackN.GetPx(); mom[4] = trackN.GetPy(); mom[5] = trackN.GetPz();
 381         }else{
 382           //Get position of the secondary vertex
 383           vertex2.GetXYZ(v2[0],v2[1],v2[2]);
 384           vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
 385           vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
 386         }
 387         goodVtx1 = kTRUE;
 388
 389         // no vertexing if DeltaMass > fMassCut
 390         if(fVertexOnTheFly) {
 391           goodVtx1 = kFALSE;
 392           if(SelectInvMass(mom)) {
 393             // primary vertex from *other* tracks in the event
 394             vertexer1->SetFieldkG(event->GetMagneticField());
 395             skipped[0] = trkEntryP[iTrkP];
 396             skipped[1] = trkEntryN[iTrkN];
 397             vertexer1->SetSkipTracks(2,skipped);
 398             AliESDVertex *vertex1onfly =
 399               (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
 400             if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
 401             vertex1onfly->GetXYZ(fV1);
 402             //vertex1onfly->PrintStatus();
 403             delete vertex1onfly;
 404           }
 405         }
 406
 407 /*      if (fKFSecondVertex&&fKFTopConstr&&fKFPrimVertex){
 408 //====================== TOPOLOGICAL CONSTRAINT !!=====================================
 409         //Primary vertex constructed from ESD using KF methods!!!
 410         AliKFVertex primVtxCopy(*(event->GetPrimaryVertex()));
 411
 412         //Subtract Daughters from primary vertex
 413         primVtxCopy -= trackP;
 414         primVtxCopy -= trackN;
 415
 416         //Add V0 to the vertex in order to improve primary vertex resolution
 417         primVtxCopy += V0;
 418
 419         //Set production vertex for V0
 420         V0.SetProductionVertex(primVtxCopy);
 421
 422         //Recalculate primary vertex
 423         fV1[0] = primVtxCopy.GetX();
 424         fV1[1] = primVtxCopy.GetY();
 425         fV1[2] = primVtxCopy.GetZ();
 426 //=====================================================================================
 427         }*/
 428
 429         // impact parameters of the tracks w.r.t. the primary vertex
 430         d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
 431         d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
 432
 433         // create the object AliBtoJPSItoEle
 434         AliBtoJPSItoEle theB(ev,trkNum,fV1,v2,dca,mom,d0);
 435         // select B's
 436         if(goodVtx1 && theB.Select(fBCuts,okB)) {
 437           // get PID info from ESD
 438           AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
 439           Double_t esdpid0[5];
 440           t0->GetESDpid(esdpid0);
 441           if(t0->GetStatus()&AliESDtrack::kTOFpid) {
 442             tofmass[0] = CalculateTOFmass(t0->GetP(),
 443                                           t0->GetIntegratedLength(),
 444                                           t0->GetTOFsignal());
 445           } else {
 446             tofmass[0] = -1000.;
 447           }
 448           AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
 449           Double_t esdpid1[5];
 450           t1->GetESDpid(esdpid1);
 451           if(t1->GetStatus()&AliESDtrack::kTOFpid) {
 452             tofmass[1] = CalculateTOFmass(t1->GetP(),
 453                                           t1->GetIntegratedLength(),
 454                                           t1->GetTOFsignal());
 455           } else {
 456             tofmass[1] = -1000.;
 457           }
 458
 459           theB.SetPIDresponse(esdpid0,esdpid1);
 460           theB.SetTOFmasses(tofmass);
 461
 462           // fill the tree
 463           ioBtoJPSItoEle=&theB;
 464           treeB->Fill();
 465
 466           nBrec++; nBrec1ev++;
 467           ioBtoJPSItoEle=0;
 468         }
 469
 470         negtrack = 0;
 471       } // loop on negative tracks
 472       postrack = 0;
 473     }   // loop on positive tracks
 474
 475     delete [] trkEntryP;
 476     delete [] trkEntryN;
 477     delete trkTree;
 478
 479     printf(" Number of B candidates: %d\n",nBrec1ev);
 480   }    // loop on events in file
 481
 482
 483   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
 484   printf("\n+++\n+++ Total number of B candidates: %d\n+++\n",nBrec);
 485
 486   delete vertexer1;
 487
 488   esdFile->Close();
 489
 490   // create a copy of this class to be written to output file
 491   AliBtoJPSItoEleAnalysis *copy = (AliBtoJPSItoEleAnalysis*)this->Clone("BtoJPSItoEleAnalysis");
 492
 493   // add PDG codes to decay tracks in found candidates (in simulation mode)
 494   // and store tree in the output file
 495   if(!fSim) {
 496     TFile *outroot = new TFile(outName1.Data(),"recreate");
 497     treeB->Write();
 498     copy->Write();
 499     outroot->Close();
 500     delete outroot;
 501   } else {
 502     printf(" Now adding information from simulation (PDG codes) ...\n");
 503     TTree *treeBsim = new TTree("TreeB","Tree with B candidates");
 504     SimulationInfo(treeB,treeBsim);
 505     delete treeB;
 506     TFile *outroot = new TFile(outName1.Data(),"recreate");
 507     treeBsim->Write();
 508     copy->Write();
 509     outroot->Close();
 510     delete outroot;
 511   }
 512
 513   // write to a file the total number of events
 514   FILE *outdat = fopen(outName2.Data(),"w");
 515   fprintf(outdat,"%d\n",nTotEv);
 516   fclose(outdat);
 517
 518   return;
 519 }
 520 //-----------------------------------------------------------------------------
 521 void AliBtoJPSItoEleAnalysis::PrintStatus() const {
 522   // Print parameters being used
 523
 524   printf("Preselections:\n");
 525   printf("    fPtCut   = %f GeV\n",fPtCut);
 526   printf("    fd0Cut   = %f micron\n",fd0Cut);
 527   printf("    fMassCut = %f GeV\n",fMassCut);
 528   printf("    fPidCut  > %f \n",fPidCut);
 529   if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
 530   if(fSim) {
 531     printf("Simulation mode\n");
 532     if(fOnlySignal && !(fOnlyPrimaryJpsi)) printf("  Only signal goes to file\n");
 533     if(fOnlyPrimaryJpsi && !(fOnlySignal)) printf("  Only primary Jpsi go to file\n");
 534     if(fOnlyPrimaryJpsi && fOnlySignal) printf("  Both signal and primary Jpsi go to file\n");
 535   }
 536   printf("Cuts on candidates:\n");
 537   printf("    |M-MJPsi| [GeV]  < %f\n",fBCuts[0]);
 538   printf("    dca    [micron]  < %f\n",fBCuts[1]);
 539   printf("    cosThetaStar     < %f\n",fBCuts[2]);
 540   printf("    pTP     [GeV]    > %f\n",fBCuts[3]);
 541   printf("    pTN     [GeV]    > %f\n",fBCuts[4]);
 542   printf("    |d0P|  [micron]  < %f\n",fBCuts[5]);
 543   printf("    |d0N|  [micron]  < %f\n",fBCuts[6]);
 544   printf("    d0d0  [micron^2] < %f\n",fBCuts[7]);
 545   printf("    cosThetaPoint    > %f\n",fBCuts[8]);
 546
 547   return;
 548 }
 549 //-----------------------------------------------------------------------------
 550 Bool_t AliBtoJPSItoEleAnalysis::SelectInvMass(const Double_t p[6]) const {
 551   // Apply preselection in the invariant mass of the pair
 552
 553   Double_t mJPsi = 3.096916;
 554   Double_t mel   = 0.00510998902;
 555
 556   Double_t energy[2];
 557   Double_t mom2[2],momTot2;
 558
 559   mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
 560   mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
 561
 562   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
 563             (p[1]+p[4])*(p[1]+p[4])+
 564             (p[2]+p[5])*(p[2]+p[5]);
 565
 566   // J/Psi -> e+ e-
 567   energy[1] = TMath::Sqrt(mel*mel+mom2[1]);
 568   energy[0] = TMath::Sqrt(mel*mel+mom2[0]);
 569
 570   Double_t minvJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
 571
 572   if(TMath::Abs(minvJPsi-mJPsi)  < fMassCut) return kTRUE;
 573   return kFALSE;
 574 }
 575 //-----------------------------------------------------------------------------
 576 void AliBtoJPSItoEleAnalysis::SelectTracks(AliESDEvent *event,
 577         TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
 578         TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
 579   // Create two TObjArrays with positive and negative tracks and
 580   // apply single-track preselection
 581
 582   nTrksP=0,nTrksN=0;
 583
 584   Int_t entr = event->GetNumberOfTracks();
 585
 586   // transfer ITS tracks from ESD to arrays and to a tree
 587   for(Int_t i=0; i<entr; i++) {
 588
 589     AliESDtrack *esdtrack = event->GetTrack(i);
 590     UInt_t status = esdtrack->GetStatus();
 591
 592     if(!(status&AliESDtrack::kITSin)) continue;
 593
 594     // single track selection
 595     if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
 596
 597     if(esdtrack->GetSign()<0) { // negative track
 598       trksN.AddLast(esdtrack);
 599       trkEntryN[nTrksN] = i;
 600       nTrksN++;
 601     } else {                 // positive track
 602       trksP.AddLast(esdtrack);
 603       trkEntryP[nTrksP] = i;
 604       nTrksP++;
 605     }
 606
 607   } // loop on ESD tracks
 608
 609   return;
 610 }
 611 //-----------------------------------------------------------------------------
 612 void AliBtoJPSItoEleAnalysis::SetBCuts(Double_t cut0,Double_t cut1,
 613                                    Double_t cut2,Double_t cut3,Double_t cut4,
 614                                    Double_t cut5,Double_t cut6,
 615                                    Double_t cut7,Double_t cut8) {
 616   // Set the cuts for B selection
 617   fBCuts[0] = cut0;
 618   fBCuts[1] = cut1;
 619   fBCuts[2] = cut2;
 620   fBCuts[3] = cut3;
 621   fBCuts[4] = cut4;
 622   fBCuts[5] = cut5;
 623   fBCuts[6] = cut6;
 624   fBCuts[7] = cut7;
 625   fBCuts[8] = cut8;
 626
 627   return;
 628 }
 629 //-----------------------------------------------------------------------------
 630 void AliBtoJPSItoEleAnalysis::SetBCuts(const Double_t cuts[9]) {
 631   // Set the cuts for B selection
 632
 633   for(Int_t i=0; i<9; i++) fBCuts[i] = cuts[i];
 634
 635   return;
 636 }
 637 //-----------------------------------------------------------------------------
 638 Bool_t
 639 AliBtoJPSItoEleAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
 640   // Check if track passes some kinematical cuts
 641   // Magnetic field "b" (kG)
 642
 643   if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
 644     return kFALSE;
 645   if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
 646     return kFALSE;
 647   //select only tracks with the "combined PID"
 648   UInt_t status = trk.GetStatus();
 649   if ((status&AliESDtrack::kESDpid)==0) return kTRUE;
 650   Double_t r[5];
 651   trk.GetESDpid(r);
 652   if(r[0] < fPidCut) return kFALSE;
 653
 654   return kTRUE;
 655 }
 656 //----------------------------------------------------------------------------
 657 void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *gAlice,
 658                                            Int_t evFirst,Int_t evLast) const {
 659   // Create a file with simulation info for the reconstructed tracks
 660
 661   TFile *outFile = TFile::Open("BTracksRefFile.root","recreate");
 662   TFile *esdFile = TFile::Open("AliESDs.root");
 663
 664   AliMC *mc = gAlice->GetMCApp();
 665
 666   Int_t      label;
 667   TParticle *part;
 668   TParticle *mumpart;
 669   TParticle *gmumpart;
 670   REFTRACK   reftrk;
 671
 672   AliESDEvent* event = new AliESDEvent();
 673   TTree* tree = (TTree*) esdFile->Get("esdTree");
 674   event->ReadFromTree(tree);
 675   // loop on events in file
 676   for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
 677     if(iEvent>evLast) break;
 678     tree->GetEvent(iEvent);
 679     Int_t ev = (Int_t)event->GetEventNumberInFile();
 680
 681     gAlice->GetEvent(ev);
 682
 683     Int_t nentr=(Int_t)event->GetNumberOfTracks();
 684
 685     // Tree for true track parameters
 686     char ttname[100];
 687     sprintf(ttname,"Tree_Ref_%d",ev);
 688     TTree *reftree = new TTree(ttname,"Tree with true track params");
 689     reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:gmumlab:gmumpdg:mumprongs:Vx/F:Vy:Vz:Px:Py:Pz");
 690 //    reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
 691
 692     for(Int_t i=0; i<nentr; i++) {
 693       AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
 694       label = TMath::Abs(esdtrack->GetLabel());
 695
 696       part = (TParticle*)mc->Particle(label);
 697       reftrk.lab = label;
 698       reftrk.pdg = part->GetPdgCode();
 699       reftrk.mumlab = part->GetFirstMother();
 700       if(part->GetFirstMother()>=0) {
 701         mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
 702         reftrk.mumpdg = mumpart->GetPdgCode();
 703         reftrk.mumprongs = mumpart->GetNDaughters();
 704         reftrk.gmumlab = mumpart->GetFirstMother();
 705         if(mumpart->GetFirstMother()>=0) {
 706           gmumpart = (TParticle*)gAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
 707           reftrk.gmumpdg = gmumpart->GetPdgCode();
 708         }
 709       } else {
 710         reftrk.mumpdg=-1;
 711         reftrk.mumprongs=-1;
 712         reftrk.gmumpdg=-1;
 713         reftrk.gmumlab=part->GetFirstMother(); // If it hasn't any mother, then it has neither Gmother!
 714         // reftrk.gmumlab=-1; // If it hasn't any mother, then it has neither Gmother!
 715       }
 716       reftrk.Vx = part->Vx();
 717       reftrk.Vy = part->Vy();
 718       reftrk.Vz = part->Vz();
 719       reftrk.Px = part->Px();
 720       reftrk.Py = part->Py();
 721       reftrk.Pz = part->Pz();
 722
 723       reftree->Fill();
 724
 725     } // loop on tracks
 726
 727     outFile->cd();
 728     reftree->Write();
 729
 730     delete reftree;
 731   } // loop on events
 732
 733   esdFile->Close();
 734   outFile->Close();
 735
 736   return;
 737 }
 738 //-----------------------------------------------------------------------------
 739 void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) const {
 740   // add pdg codes to candidate decay tracks (for sim)
 741
 742   TString refFileName("BTracksRefFile.root");
 743   if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
 744     printf("AliBtoJPSItoEleAnalysis::SimulationInfo: no reference file found!\n");
 745     return;
 746   }
 747   TFile *refFile = TFile::Open(refFileName.Data());
 748
 749   Char_t refTreeName[100];
 750   Int_t  event;
 751   Int_t  pdg[2],mumpdg[2],mumlab[2],gmumpdg[2],gmumlab[2];
 752   REFTRACK reftrk;
 753
 754   // read-in reference tree for event 0 (the only event for Pb-Pb)
 755   sprintf(refTreeName,"Tree_Ref_%d",0);
 756   TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
 757   refTree0->SetBranchAddress("rectracks",&reftrk);
 758
 759   AliBtoJPSItoEle *theB = 0;
 760   treeBin->SetBranchAddress("BtoJPSItoEle",&theB);
 761   treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&theB,200000,0);
 762
 763   Int_t entries = (Int_t)treeBin->GetEntries();
 764
 765   for(Int_t i=0; i<entries; i++) {
 766     if(i%100==0) printf("  done %d candidates of %d\n",i,entries);
 767
 768     treeBin->GetEvent(i);
 769     event = theB->EventNo();
 770
 771     if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
 772       refTree0->GetEvent(theB->GetTrkNum(0));
 773       pdg[0]    = reftrk.pdg;
 774       mumpdg[0] = reftrk.mumpdg;
 775       mumlab[0] = reftrk.mumlab;
 776       gmumpdg[0] = reftrk.gmumpdg;
 777       gmumlab[0] = reftrk.gmumlab;
 778       refTree0->GetEvent(theB->GetTrkNum(1));
 779       pdg[1]    = reftrk.pdg;
 780       mumpdg[1] = reftrk.mumpdg;
 781       mumlab[1] = reftrk.mumlab;
 782       gmumpdg[1] = reftrk.gmumpdg;
 783       gmumlab[1] = reftrk.gmumlab;
 784     } else {
 785       sprintf(refTreeName,"Tree_Ref_%d",event);
 786       TTree *refTree = (TTree*)refFile->Get(refTreeName);
 787       refTree->SetBranchAddress("rectracks",&reftrk);
 788       refTree->GetEvent(theB->GetTrkNum(0));
 789       pdg[0]    = reftrk.pdg;
 790       mumpdg[0] = reftrk.mumpdg;
 791       mumlab[0] = reftrk.mumlab;
 792       gmumpdg[0] = reftrk.gmumpdg;
 793       gmumlab[0] = reftrk.gmumlab;
 794       refTree->GetEvent(theB->GetTrkNum(1));
 795       pdg[1]    = reftrk.pdg;
 796       mumpdg[1] = reftrk.mumpdg;
 797       mumlab[1] = reftrk.mumlab;
 798       gmumpdg[1] = reftrk.gmumpdg;
 799       gmumlab[1] = reftrk.gmumlab;
 800       delete refTree;
 801     }
 802
 803     theB->SetPdgCodes(pdg);
 804     theB->SetMumPdgCodes(mumpdg);
 805     theB->SetGMumPdgCodes(gmumpdg);
 806
 807     if(gmumpdg[0]==gmumpdg[1] &&              // Both GrandMothers are of the same sign
 808        (TMath::Abs(gmumpdg[0])==521   || TMath::Abs(gmumpdg[0])==511   ||  // GrandMother Bplus/Bminus or B0/B0bar
 809         TMath::Abs(gmumpdg[0])==523   || TMath::Abs(gmumpdg[0])==513   ||  // B0s/B0sbar
 810         TMath::Abs(gmumpdg[0])==515   || TMath::Abs(gmumpdg[0])==525   ||  //
 811         TMath::Abs(gmumpdg[0])==531   || TMath::Abs(gmumpdg[0])==533   ||  //
 812         TMath::Abs(gmumpdg[0])==535   || TMath::Abs(gmumpdg[0])==541   ||  //
 813         TMath::Abs(gmumpdg[0])==543   || TMath::Abs(gmumpdg[0])==545   ||  //
 814         TMath::Abs(gmumpdg[0])==10521 || TMath::Abs(gmumpdg[0])==10511 ||  //   all possible
 815         TMath::Abs(gmumpdg[0])==10523 || TMath::Abs(gmumpdg[0])==10513 ||  //    B mesons
 816         TMath::Abs(gmumpdg[0])==20523 || TMath::Abs(gmumpdg[0])==20513 ||  //
 817         TMath::Abs(gmumpdg[0])==10531 || TMath::Abs(gmumpdg[0])==10533 ||  //
 818         TMath::Abs(gmumpdg[0])==20533 || TMath::Abs(gmumpdg[0])==10541 ||  //
 819         TMath::Abs(gmumpdg[0])==20543 || TMath::Abs(gmumpdg[0])==10543 ||  //
 820         TMath::Abs(gmumpdg[0])==4122  || TMath::Abs(gmumpdg[0])==4222  ||  // All possible B baryons
 821         TMath::Abs(gmumpdg[0])==4212  || TMath::Abs(gmumpdg[0])==4112  ||  // All possible B baryons
 822         TMath::Abs(gmumpdg[0])==4224  || TMath::Abs(gmumpdg[0])==4214  ||  // All possible B baryons
 823         TMath::Abs(gmumpdg[0])==4114  || TMath::Abs(gmumpdg[0])==4232  ||  // All possible B baryons
 824         TMath::Abs(gmumpdg[0])==4132  || TMath::Abs(gmumpdg[0])==4322  ||  // All possible B baryons
 825         TMath::Abs(gmumpdg[0])==4312  || TMath::Abs(gmumpdg[0])==4324  ||  // All possible B baryons
 826         TMath::Abs(gmumpdg[0])==4314  || TMath::Abs(gmumpdg[0])==4332  ||  // All possible B baryons
 827         TMath::Abs(gmumpdg[0])==4334  || TMath::Abs(gmumpdg[0])==4412  ||  // All possible B baryons
 828         TMath::Abs(gmumpdg[0])==4422  || TMath::Abs(gmumpdg[0])==4414  ||  // All possible B baryons
 829         TMath::Abs(gmumpdg[0])==4424  || TMath::Abs(gmumpdg[0])==4432  ||  // All possible B baryons
 830         TMath::Abs(gmumpdg[0])==4434  || TMath::Abs(gmumpdg[0])==4444      // All possible B baryons
 831        ) &&
 832        mumpdg[0]==443 && mumpdg[1]== 443 &&
 833        mumlab[0]==mumlab[1] &&
 834        reftrk.mumprongs==2 &&
 835        pdg[0]==-11 && pdg[1]==11
 836      ) theB->SetSignal();
 837
 838     else if (  // here consider the case of primary J/psi
 839        mumpdg[0]==443 && mumpdg[1]== 443 &&
 840        pdg[0]==-11 && pdg[1]==11 &&
 841        mumlab[0]==mumlab[1] &&
 842        reftrk.mumprongs==2 &&
 843        ( gmumlab[0]<0                     ||   // really primary J/psi (without family. e.g. from Cocktail)
 844          TMath::Abs(gmumpdg[0])==100443   ||   // from Psi(2S)
 845          TMath::Abs(gmumpdg[0])==10441    ||   // from Csi_c0(1P)
 846          TMath::Abs(gmumpdg[0])==20443    ||   // from Csi_c1(1P)
 847          TMath::Abs(gmumpdg[0])==10443    ||   // from h_c(1P)
 848          TMath::Abs(gmumpdg[0])==445      ||   // from Csi_c2(1P)
 849          (gmumpdg[0]>=81 && gmumpdg[0]<=100)   // this is for MC internal use (e.g. J/psi from string)
 850        )
 851      ) theB->SetJpsiPrimary();
 852
 853     // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();
 854
 855  // write it out 1) always if you have not asked for onlySignal or OnlyPrimaryJpsi (or both)
 856  //          or  2) if you have asked for Signal and it is Signal
 857  //          or  3) if you have asked for Primary Jpsi and it is a Primary Jpsi
 858     if ( (!fOnlySignal && !fOnlyPrimaryJpsi) || (fOnlySignal && theB->IsSignal())
 859         || (fOnlyPrimaryJpsi && theB->IsJpsiPrimary()) ) treeBout->Fill();
 860   }
 861
 862   delete refTree0;
 863
 864   refFile->Close();
 865
 866   return;
 867 }
 868
 869
 870
 871
 872