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   // fix for temporary bug in ESDfilter
 198   // the AODs with null vertex pointer didn't pass the PhysSel
 199   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 200
 201
 202   // AOD primary vertex
 203   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 204
 205   // load MC particles
 206   TClonesArray* mcArray =
 207      dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 208   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 209      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
 210
 211   // Read AOD Monte Carlo information
 212   if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);
 213
 214   // loop over J/Psi->ee candidates
 215   Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
 216   printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);
 217
 218   for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
 219     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
 220
 221     Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 222
 223     //Secondary vertex
 224     Double_t vtxSec[3] = {0.,0.,0.};
 225     Double_t vtxPrim[3] = {0.,0.,0.};
 226     d->GetSecondaryVtx(vtxSec);
 227     Bool_t unsetvtx=kFALSE;
 228     if(!d->GetOwnPrimaryVtx()) {
 229       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 230       unsetvtx=kTRUE;
 231     }
 232     //Int_t okJPSI=0;
 233     // here analyze only if cuts are passed
 234
 235     if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
 236       if (mcLabel == -1)
 237          {
 238           AliDebug(2,"No MC particle found");
 239          }
 240       else {
 241
 242       fhInvMass->Fill(d->InvMassJPSIee());
 243       fhD0->Fill(10000*d->ImpParXY());
 244       fhD0D0->Fill(1e8*d->Prodd0d0());
 245       fhCosThetaStar->Fill(d->CosThetaStarJPSI());
 246       fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0());
 247       fhCosThetaPointing->Fill(d->CosPointingAngle());
 248
 249       // compute X variable
 250       AliAODVertex* primVertex = d->GetOwnPrimaryVtx();
 251       vtxPrim[0] = primVertex->GetX();
 252       vtxPrim[1] = primVertex->GetY();
 253       vtxPrim[2] = primVertex->GetZ();
 254       Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt();
 255       Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt();
 256
 257       fhDecayTime->Fill(10000*pseudoProperDecayTime);
 258       // Post the data already here
 259       PostData(1,fOutput);
 260
 261       fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values
 262
 263       }//
 264
 265      } // end of JPSItoEle candidates selection according to cuts
 266
 267     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 268
 269    }// end loop on JPSI to ele candidates
 270
 271 }
 272 //_________________________________________________________________________________________________
 273 void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/)
 274 {
 275   //
 276   // Terminate analysis
 277   //
 278
 279   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n");
 280
 281   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 282   if (!fOutput) {
 283     printf("ERROR: fOutput not available\n");
 284     return;
 285   }
 286
 287   if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
 288   fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
 289   fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
 290   fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
 291   fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
 292   fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
 293   fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
 294   fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
 295   fNtupleJPSI = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleJPSI"));
 296
 297   if(fOkMinimize){
 298      Double_t* pseudoProperValues=0x0;
 299      Double_t* invMassValues=0x0;
 300      Int_t ncand=0;
 301      fCdfFit = new AliAnalysisBtoJPSItoEle();
 302
 303      printf("+++\n+++  AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n");
 304      fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand);
 305      printf("+++\n+++  X and M vectors filled starting from N-tuple ---> OK\n+++\n");
 306      printf("+++\n+++  Number of candidates ---> %d J/psi \n+++\n",ncand);
 307      fCdfFit->SetPtBin(0);
 308      printf("+++\n+++  Pt bin setted n. ---> %d  \n+++\n",fCdfFit->GetPtBin());
 309
 310      if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);}
 311         else {
 312           //SetPathMCfile(".");
 313           TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin());
 314           fCdfFit->CloneMCtemplate(hLocal);
 315         }
 316
 317      printf("+++\n+++  MC template histo copied ---> OK\n+++\n");
 318      //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
 319      fCdfFit->DoMinimization();
 320    }
 321
 322   return;
 323 }
 324 //_________________________________________________________________________________________________
 325 void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
 326 {
 327   //
 328   // apply Pt cuts (lower cut, upper cut)
 329   //
 330   for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
 331 }
 332 //_________________________________________________________________________________________________
 333 void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9])
 334 {
 335   //
 336   // apply D0 and JPSI cuts
 337   //
 338   for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
 339 }
 340 //_________________________________________________________________________________________________
 341 void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
 342 {
 343   //
 344   // Reads MC information if needed (for example in case of parametrized PID)
 345   //
 346
 347   // load MC particles
 348   TClonesArray* mcArray =
 349      dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 350   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 351      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
 352
 353   // load MC header
 354   AliAODMCHeader* mcHeader =
 355     (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 356   if(!mcHeader){
 357     printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
 358   }
 359
 360   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 361
 362   //Double_t rmax = 0.000005;
 363   Double_t mcPrimVtx[3];
 364
 365   // get MC primary Vtx
 366   mcHeader->GetVertex(mcPrimVtx);
 367
 368   Int_t nInCandidates = inputArray->GetEntriesFast();
 369   printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
 370
 371   Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
 372   Int_t labJPSIdaugh0=0;
 373   Int_t labJPSIdaugh1=0;
 374
 375   for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
 376
 377     AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
 378     Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 379
 380     Double_t vtxPrim[3] = {0.,0.,0.};
 381     Double_t vtxSec[3] = {0.,0.,0.};
 382     dd->GetSecondaryVtx(vtxSec);
 383     Bool_t unsetvtx=kFALSE;
 384     if(!dd->GetOwnPrimaryVtx()) {
 385       dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 386       unsetvtx=kTRUE;
 387     }
 388
 389     if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
 390
 391       AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
 392       AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
 393       lab0 = trk0->GetLabel();
 394       lab1 = trk1->GetLabel();
 395
 396       AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
 397       if(!part0) {
 398         printf("no MC particle\n");
 399         continue;
 400       }
 401
 402       while(part0->GetMother()>=0) {
 403        labMother=part0->GetMother();
 404        part0 = (AliAODMCParticle*)mcArray->At(labMother);
 405        if(!part0) {
 406          printf("no MC mother particle\n");
 407          break;
 408        }
 409        pdgMother = TMath::Abs(part0->GetPdgCode());
 410        if(pdgMother==443) {//this for JPSI
 411          labJPSIdaugh0=labMother;
 412          break;
 413        }
 414       }
 415
 416       AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
 417       if(!part1) {
 418         printf("no MC particle\n");
 419         continue;
 420       }
 421
 422       while(part1->GetMother()>=0) {
 423        labMother=part1->GetMother();
 424        part1 = (AliAODMCParticle*)mcArray->At(labMother);
 425        if(!part1) {
 426          printf("no MC mother particle\n");
 427          break;
 428        }
 429        pdgMother = TMath::Abs(part1->GetPdgCode());
 430        if(pdgMother==443) {//this for JPSI
 431          labJPSIdaugh1=labMother;
 432          break;
 433        }
 434       }
 435
 436       if (mcLabelSec == -1)
 437          {
 438           AliDebug(2,"No MC particle found");
 439          }
 440       else {
 441
 442       if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
 443           AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
 444           pdgJPSI = partJPSI->GetPdgCode();
 445           Int_t pdaugh0 = partJPSI->GetDaughter(0);
 446           Int_t pdaugh1 = partJPSI->GetDaughter(1);
 447           AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
 448           AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
 449           pdg0 = partDaugh0->GetPdgCode();
 450           pdg1 = partDaugh1->GetPdgCode();
 451         if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
 452            labJPSIMother = partJPSI->GetMother();
 453            AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
 454            pdgJPSIMother = partJPSIMother->GetPdgCode();
 455
 456         // chech if JPSI comes from B
 457            if( pdgJPSIMother==511   || pdgJPSIMother==521   ||
 458                pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
 459                pdgJPSIMother==513   || pdgJPSIMother==523   ||
 460                pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
 461                pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
 462                pdgJPSIMother==515   || pdgJPSIMother==525   ||
 463                pdgJPSIMother==531   || pdgJPSIMother==10531 ||
 464                pdgJPSIMother==533   || pdgJPSIMother==10533 ||
 465                pdgJPSIMother==20533 || pdgJPSIMother==535   ||
 466                pdgJPSIMother==541   || pdgJPSIMother==10541 ||
 467                pdgJPSIMother==543   || pdgJPSIMother==10543 ||
 468                pdgJPSIMother==20543 || pdgJPSIMother==545)
 469                { //this is for MC JPSI -> ee from B-hadrons
 470
 471                   fhInvMass->Fill(dd->InvMassJPSIee());
 472                   fhD0->Fill(10000*dd->ImpParXY());
 473                   fhD0D0->Fill(1e8*dd->Prodd0d0());
 474                   fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
 475                   fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
 476                   fhCosThetaPointing->Fill(dd->CosPointingAngle());
 477
 478                   // compute X variable
 479                   AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
 480                   vtxPrim[0] = primVertex->GetX();
 481                   vtxPrim[1] = primVertex->GetY();
 482                   vtxPrim[2] = primVertex->GetZ();
 483                   Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
 484                   Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
 485
 486                   fhDecayTime->Fill(10000*pseudoProperDecayTime);
 487                   // Post the data already here
 488                   PostData(1,fOutput);
 489
 490                   fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values
 491
 492                   Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
 493                   Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
 494                   fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
 495
 496                   // Post the data already here
 497                   PostData(1,fOutput);
 498
 499                } //this is for MC JPSI -> ee from B-hadrons
 500             } //this is for MC JPSI -> ee
 501         }
 502        } // end of check if JPSI is signal
 503       }
 504     }
 505
 506 }
 507 //_________________________________________________________________________________________________
 508 TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum)
 509 {
 510   //
 511   // Open a local file with MC x distribution for JPSI from B-hadron
 512   //
 513
 514   TH1F* hCsiMC = new TH1F();
 515   TString useFile = datadir.Data();
 516   useFile.Append("/CsiMCfromKine_10PtBins.root");
 517   TFile *f = new TFile(useFile);
 518   Double_t scale = 0.;
 519   char processCsiMCXHisto[1024];
 520   if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper");
 521   else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum);
 522   hCsiMC = (TH1F*)f->Get(processCsiMCXHisto);
 523   scale=1/hCsiMC->Integral();
 524   hCsiMC->Scale(scale);
 525   printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries());
 526
 527   return hCsiMC;
 528
 529 }
 530