PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 //                   Class AliAnalysisTaskSEBtoJPSItoEle
  17 //                AliAnalysisTaskSE for the reconstruction
  18 //                   of heavy-flavour decay candidates
  19 //            Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
  20 //*************************************************************************
  21 class AliAnalysisTaskSE;
  22 class AliESDEvent;
  23 class AliVEvent;
  24 class iostream;
  25 class TSystem;
  26 class TROOT;
  27 #include <TFile.h>
  28 #include <TClonesArray.h>
  29 #include <TDatabasePDG.h>
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliAODHandler.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODMCParticle.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODRecoDecayHF2Prong.h"
  37
  38 #include "AliAnalysisBtoJPSItoEle.h"
  39 #include "AliAnalysisTaskSEBtoJPSItoEle.h"
  40
  41 ClassImp(AliAnalysisTaskSEBtoJPSItoEle)
  42
  43 //______________________________________________________________________________
  44 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() :
  45 AliAnalysisTaskSE(),
  46 fOutput(0),
  47 fCdfFit(0),
  48 fNtupleJPSI(0),
  49 fhDecayTimeMCjpsifromB(0),
  50 fhDecayTime(0),
  51 fhInvMass(0),
  52 fhD0(0),
  53 fhD0D0(0),
  54 fhCosThetaStar(0),
  55 fhCosThetaPointing(0),
  56 fhCtsVsD0D0(0),
  57 fOkAODMC(kFALSE),
  58 fNameMCfile(""),
  59 fOkMinimize(kFALSE)
  60 {
  61   // default constructor
  62 }
  63 //_________________________________________________________________________________________________
  64 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) :
  65 AliAnalysisTaskSE(name),
  66 fOutput(0),
  67 fCdfFit(0),
  68 fNtupleJPSI(0),
  69 fhDecayTimeMCjpsifromB(0),
  70 fhDecayTime(0),
  71 fhInvMass(0),
  72 fhD0(0),
  73 fhD0D0(0),
  74 fhCosThetaStar(0),
  75 fhCosThetaPointing(0),
  76 fhCtsVsD0D0(0),
  77 fOkAODMC(kFALSE),
  78 fNameMCfile(""),
  79 fOkMinimize(kFALSE)
  80 {
  81   // standard constructor
  82
  83   // Output slot #1 writes into a TList container
  84   DefineOutput(1,TList::Class());  //My private output
  85 }
  86 //_________________________________________________________________________________________________
  87 AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle()
  88 {
  89   // destructor
  90
  91   if (fOutput) {
  92     delete fOutput;
  93     fOutput = 0;
  94   }
  95 }
  96 //_________________________________________________________________________________________________
  97 void AliAnalysisTaskSEBtoJPSItoEle::Init()
  98 {
  99   // Initialization
 100
 101   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n");
 102
 103   return;
 104
 105 }
 106 //_________________________________________________________________________________________________
 107 void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects()
 108 {
 109   // Create the output container
 110
 111   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n");
 112
 113   // Several histograms are more conveniently managed in a TList
 114   fOutput = new TList();
 115   fOutput->SetOwner();
 116
 117   if(fOkAODMC){
 118   fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
 119   fhDecayTimeMCjpsifromB->Sumw2();
 120   fhDecayTimeMCjpsifromB->SetMinimum(0);
 121   fOutput->Add(fhDecayTimeMCjpsifromB);
 122   }
 123
 124   fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
 125   fhDecayTime->Sumw2();
 126   fhDecayTime->SetMinimum(0);
 127   fOutput->Add(fhDecayTime);
 128
 129   fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2);
 130   fhInvMass->Sumw2();
 131   fhInvMass->SetMinimum(0);
 132   fOutput->Add(fhInvMass);
 133
 134   fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.);
 135   fhD0->Sumw2();
 136   fhD0->SetMinimum(0);
 137   fOutput->Add(fhD0);
 138
 139   fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
 140   fhD0D0->Sumw2();
 141   fhD0D0->SetMinimum(0);
 142   fOutput->Add(fhD0D0);
 143
 144   fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
 145   fhCosThetaStar->Sumw2();
 146   fhCosThetaStar->SetMinimum(0);
 147   fOutput->Add(fhCosThetaStar);
 148
 149   fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
 150   fhCosThetaPointing->Sumw2();
 151   fhCosThetaPointing->SetMinimum(0);
 152   fOutput->Add(fhCosThetaPointing);
 153
 154   fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
 155   fhCtsVsD0D0->Sumw2();
 156   fhCtsVsD0D0->SetMinimum(0);
 157   fOutput->Add(fhCtsVsD0D0);
 158
 159   fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
 160   fOutput->Add(fNtupleJPSI);
 161
 162   return;
 163 }
 164 //_________________________________________________________________________________________________
 165 void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/)
 166 {
 167   // Execute analysis for current event:
 168
 169   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 170
 171   TClonesArray *inputArrayJPSItoEle = 0;
 172
 173   if(!aod && AODEvent() && IsStandardAOD()) {
 174     // In case there is an AOD handler writing a standard AOD, use the AOD
 175     // event in memory rather than the input (ESD) event.
 176     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 177     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 178     // have to taken from the AOD event hold by the AliAODExtension
 179     AliAODHandler* aodHandler = (AliAODHandler*)
 180       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 181     if(aodHandler->GetExtensions()) {
 182       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 183       AliAODEvent *aodFromExt = ext->GetAOD();
 184       // load Jpsi candidates
 185       inputArrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
 186     }
 187   } else {
 188     // load Jpsi candidates
 189     inputArrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
 190   }
 191
 192   if(!inputArrayJPSItoEle) {
 193     printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n");
 194     return;
 195   }
 196
 197   // AOD primary vertex
 198   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 199
 200   // load MC particles
 201   TClonesArray* mcArray =
 202      dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 203   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 204      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
 205
 206   // Read AOD Monte Carlo information
 207   if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);
 208
 209   // loop over J/Psi->ee candidates
 210   Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
 211   printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);
 212
 213   for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
 214     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
 215
 216     Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 217
 218     //Secondary vertex
 219     Double_t vtxSec[3] = {0.,0.,0.};
 220     Double_t vtxPrim[3] = {0.,0.,0.};
 221     d->GetSecondaryVtx(vtxSec);
 222     Bool_t unsetvtx=kFALSE;
 223     if(!d->GetOwnPrimaryVtx()) {
 224       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 225       unsetvtx=kTRUE;
 226     }
 227     //Int_t okJPSI=0;
 228     // here analyze only if cuts are passed
 229
 230     if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
 231       if (mcLabel == -1)
 232          {
 233           AliDebug(2,"No MC particle found");
 234          }
 235       else {
 236
 237       fhInvMass->Fill(d->InvMassJPSIee());
 238       fhD0->Fill(10000*d->ImpParXY());
 239       fhD0D0->Fill(1e8*d->Prodd0d0());
 240       fhCosThetaStar->Fill(d->CosThetaStarJPSI());
 241       fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0());
 242       fhCosThetaPointing->Fill(d->CosPointingAngle());
 243
 244       // compute X variable
 245       AliAODVertex* primVertex = d->GetOwnPrimaryVtx();
 246       vtxPrim[0] = primVertex->GetX();
 247       vtxPrim[1] = primVertex->GetY();
 248       vtxPrim[2] = primVertex->GetZ();
 249       Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt();
 250       Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt();
 251
 252       fhDecayTime->Fill(10000*pseudoProperDecayTime);
 253       // Post the data already here
 254       PostData(1,fOutput);
 255
 256       fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values
 257
 258       }//
 259
 260      } // end of JPSItoEle candidates selection according to cuts
 261
 262     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 263
 264    }// end loop on JPSI to ele candidates
 265
 266 }
 267 //_________________________________________________________________________________________________
 268 void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/)
 269 {
 270   //
 271   // Terminate analysis
 272   //
 273
 274   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n");
 275
 276   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 277   if (!fOutput) {
 278     printf("ERROR: fOutput not available\n");
 279     return;
 280   }
 281
 282   if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
 283   fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
 284   fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
 285   fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
 286   fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
 287   fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
 288   fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
 289   fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
 290   fNtupleJPSI = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleJPSI"));
 291
 292   if(fOkMinimize){
 293      Double_t* pseudoProperValues=0x0;
 294      Double_t* invMassValues=0x0;
 295      Int_t ncand=0;
 296      fCdfFit = new AliAnalysisBtoJPSItoEle();
 297
 298      printf("+++\n+++  AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n");
 299      fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand);
 300      printf("+++\n+++  X and M vectors filled starting from N-tuple ---> OK\n+++\n");
 301      printf("+++\n+++  Number of candidates ---> %d J/psi \n+++\n",ncand);
 302      fCdfFit->SetPtBin(0);
 303      printf("+++\n+++  Pt bin setted n. ---> %d  \n+++\n",fCdfFit->GetPtBin());
 304
 305      if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);}
 306         else {
 307           //SetPathMCfile(".");
 308           TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin());
 309           fCdfFit->CloneMCtemplate(hLocal);
 310         }
 311
 312      printf("+++\n+++  MC template histo copied ---> OK\n+++\n");
 313      //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
 314      fCdfFit->DoMinimization();
 315    }
 316
 317   return;
 318 }
 319 //_________________________________________________________________________________________________
 320 void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
 321 {
 322   //
 323   // apply Pt cuts (lower cut, upper cut)
 324   //
 325   for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
 326 }
 327 //_________________________________________________________________________________________________
 328 void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9])
 329 {
 330   //
 331   // apply D0 and JPSI cuts
 332   //
 333   for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
 334 }
 335 //_________________________________________________________________________________________________
 336 void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
 337 {
 338   //
 339   // Reads MC information if needed (for example in case of parametrized PID)
 340   //
 341
 342   // load MC particles
 343   TClonesArray* mcArray =
 344      dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 345   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 346      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
 347
 348   // load MC header
 349   AliAODMCHeader* mcHeader =
 350     (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 351   if(!mcHeader){
 352     printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
 353   }
 354
 355   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 356
 357   //Double_t rmax = 0.000005;
 358   Double_t mcPrimVtx[3];
 359
 360   // get MC primary Vtx
 361   mcHeader->GetVertex(mcPrimVtx);
 362
 363   Int_t nInCandidates = inputArray->GetEntriesFast();
 364   printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
 365
 366   Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
 367   Int_t labJPSIdaugh0=0;
 368   Int_t labJPSIdaugh1=0;
 369
 370   for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
 371
 372     AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
 373     Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 374
 375     Double_t vtxPrim[3] = {0.,0.,0.};
 376     Double_t vtxSec[3] = {0.,0.,0.};
 377     dd->GetSecondaryVtx(vtxSec);
 378     Bool_t unsetvtx=kFALSE;
 379     if(!dd->GetOwnPrimaryVtx()) {
 380       dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 381       unsetvtx=kTRUE;
 382     }
 383
 384     if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
 385
 386       AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
 387       AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
 388       lab0 = trk0->GetLabel();
 389       lab1 = trk1->GetLabel();
 390
 391       AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
 392       if(!part0) {
 393         printf("no MC particle\n");
 394         continue;
 395       }
 396
 397       while(part0->GetMother()>=0) {
 398        labMother=part0->GetMother();
 399        part0 = (AliAODMCParticle*)mcArray->At(labMother);
 400        if(!part0) {
 401          printf("no MC mother particle\n");
 402          break;
 403        }
 404        pdgMother = TMath::Abs(part0->GetPdgCode());
 405        if(pdgMother==443) {//this for JPSI
 406          labJPSIdaugh0=labMother;
 407          break;
 408        }
 409       }
 410
 411       AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
 412       if(!part1) {
 413         printf("no MC particle\n");
 414         continue;
 415       }
 416
 417       while(part1->GetMother()>=0) {
 418        labMother=part1->GetMother();
 419        part1 = (AliAODMCParticle*)mcArray->At(labMother);
 420        if(!part1) {
 421          printf("no MC mother particle\n");
 422          break;
 423        }
 424        pdgMother = TMath::Abs(part1->GetPdgCode());
 425        if(pdgMother==443) {//this for JPSI
 426          labJPSIdaugh1=labMother;
 427          break;
 428        }
 429       }
 430
 431       if (mcLabelSec == -1)
 432          {
 433           AliDebug(2,"No MC particle found");
 434          }
 435       else {
 436
 437       if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
 438           AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
 439           pdgJPSI = partJPSI->GetPdgCode();
 440           Int_t pdaugh0 = partJPSI->GetDaughter(0);
 441           Int_t pdaugh1 = partJPSI->GetDaughter(1);
 442           AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
 443           AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
 444           pdg0 = partDaugh0->GetPdgCode();
 445           pdg1 = partDaugh1->GetPdgCode();
 446         if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
 447            labJPSIMother = partJPSI->GetMother();
 448            AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
 449            pdgJPSIMother = partJPSIMother->GetPdgCode();
 450
 451         // chech if JPSI comes from B
 452            if( pdgJPSIMother==511   || pdgJPSIMother==521   ||
 453                pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
 454                pdgJPSIMother==513   || pdgJPSIMother==523   ||
 455                pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
 456                pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
 457                pdgJPSIMother==515   || pdgJPSIMother==525   ||
 458                pdgJPSIMother==531   || pdgJPSIMother==10531 ||
 459                pdgJPSIMother==533   || pdgJPSIMother==10533 ||
 460                pdgJPSIMother==20533 || pdgJPSIMother==535   ||
 461                pdgJPSIMother==541   || pdgJPSIMother==10541 ||
 462                pdgJPSIMother==543   || pdgJPSIMother==10543 ||
 463                pdgJPSIMother==20543 || pdgJPSIMother==545)
 464                { //this is for MC JPSI -> ee from B-hadrons
 465
 466                   fhInvMass->Fill(dd->InvMassJPSIee());
 467                   fhD0->Fill(10000*dd->ImpParXY());
 468                   fhD0D0->Fill(1e8*dd->Prodd0d0());
 469                   fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
 470                   fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
 471                   fhCosThetaPointing->Fill(dd->CosPointingAngle());
 472
 473                   // compute X variable
 474                   AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
 475                   vtxPrim[0] = primVertex->GetX();
 476                   vtxPrim[1] = primVertex->GetY();
 477                   vtxPrim[2] = primVertex->GetZ();
 478                   Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
 479                   Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
 480
 481                   fhDecayTime->Fill(10000*pseudoProperDecayTime);
 482                   // Post the data already here
 483                   PostData(1,fOutput);
 484
 485                   fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values
 486
 487                   Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
 488                   Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
 489                   fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
 490
 491                   // Post the data already here
 492                   PostData(1,fOutput);
 493
 494                } //this is for MC JPSI -> ee from B-hadrons
 495             } //this is for MC JPSI -> ee
 496         }
 497        } // end of check if JPSI is signal
 498       }
 499     }
 500
 501 }
 502 //_________________________________________________________________________________________________
 503 TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum)
 504 {
 505   //
 506   // Open a local file with MC x distribution for JPSI from B-hadron
 507   //
 508
 509   TH1F* hCsiMC = new TH1F();
 510   TString useFile = datadir.Data();
 511   useFile.Append("/CsiMCfromKine_10PtBins.root");
 512   TFile *f = new TFile(useFile);
 513   Double_t scale = 0.;
 514   char processCsiMCXHisto[1024];
 515   if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper");
 516   else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum);
 517   hCsiMC = (TH1F*)f->Get(processCsiMCXHisto);
 518   scale=1/hCsiMC->Integral();
 519   hCsiMC->Scale(scale);
 520   printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries());
 521
 522   return hCsiMC;
 523
 524 }
 525