[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEBtoJPSItoEle.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
//*************************************************************************
//                   Class AliAnalysisTaskSEBtoJPSItoEle
//                AliAnalysisTaskSE for the reconstruction 
//                   of heavy-flavour decay candidates
//            Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
//*************************************************************************
class AliAnalysisTaskSE;
class AliESDEvent;
class AliVEvent;
class iostream;
class TSystem;
class TROOT;
#include <TFile.h>
#include <TClonesArray.h>
#include <TDatabasePDG.h>

#include "AliAnalysisManager.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "AliAODRecoDecayHF2Prong.h"

#include "AliAnalysisBtoJPSItoEle.h"
#include "AliAnalysisTaskSEBtoJPSItoEle.h"

ClassImp(AliAnalysisTaskSEBtoJPSItoEle)

//______________________________________________________________________________
AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() :
AliAnalysisTaskSE(),
fOutput(0),
fCdfFit(0),     
fNtupleJPSI(0),
fhDecayTimeMCjpsifromB(0),
fhDecayTime(0),                         
fhInvMass(0),                           
fhD0(0),                                
fhD0D0(0),                              
fhCosThetaStar(0),                      
fhCosThetaPointing(0),                  
fhCtsVsD0D0(0),
fOkAODMC(kFALSE),
fNameMCfile(""),
fOkMinimize(kFALSE)
{
  // default constructor
}                 
//_________________________________________________________________________________________________
AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) :
AliAnalysisTaskSE(name),
fOutput(0),
fCdfFit(0),
fNtupleJPSI(0),
fhDecayTimeMCjpsifromB(0),
fhDecayTime(0),                         
fhInvMass(0),                           
fhD0(0),                                
fhD0D0(0),                              
fhCosThetaStar(0),                      
fhCosThetaPointing(0),                  
fhCtsVsD0D0(0),
fOkAODMC(kFALSE),
fNameMCfile(""),
fOkMinimize(kFALSE)
{
  // standard constructor

  // Output slot #1 writes into a TList container
  DefineOutput(1,TList::Class());  //My private output
}
//_________________________________________________________________________________________________
AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle()
{
  // destructor

  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::Init()
{
  // Initialization

  if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n");

  return;

}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() 
{
  // Create the output container

  if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n");

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList();
  fOutput->SetOwner();

  if(fOkAODMC){
  fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
  fhDecayTimeMCjpsifromB->Sumw2();
  fhDecayTimeMCjpsifromB->SetMinimum(0);
  fOutput->Add(fhDecayTimeMCjpsifromB);
  }

  fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
  fhDecayTime->Sumw2();
  fhDecayTime->SetMinimum(0);
  fOutput->Add(fhDecayTime);

  fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2);
  fhInvMass->Sumw2();
  fhInvMass->SetMinimum(0);
  fOutput->Add(fhInvMass);

  fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.);
  fhD0->Sumw2();
  fhD0->SetMinimum(0);
  fOutput->Add(fhD0);

  fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
  fhD0D0->Sumw2();
  fhD0D0->SetMinimum(0);
  fOutput->Add(fhD0D0);

  fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
  fhCosThetaStar->Sumw2();
  fhCosThetaStar->SetMinimum(0);
  fOutput->Add(fhCosThetaStar);

  fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
  fhCosThetaPointing->Sumw2();
  fhCosThetaPointing->SetMinimum(0);
  fOutput->Add(fhCosThetaPointing);

  fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
  fhCtsVsD0D0->Sumw2();
  fhCtsVsD0D0->SetMinimum(0);
  fOutput->Add(fhCtsVsD0D0);

  fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
  fOutput->Add(fNtupleJPSI);

  return;
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event:

  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());

  TClonesArray *inputArrayJPSItoEle = 0;

  if(!aod && AODEvent() && IsStandardAOD()) {
    // In case there is an AOD handler writing a standard AOD, use the AOD 
    // event in memory rather than the input (ESD) event.    
    aod = dynamic_cast<AliAODEvent*> (AODEvent());
    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
    // have to taken from the AOD event hold by the AliAODExtension
    AliAODHandler* aodHandler = (AliAODHandler*) 
      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
    if(aodHandler->GetExtensions()) {
      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
      AliAODEvent *aodFromExt = ext->GetAOD();
      // load Jpsi candidates   
      inputArrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
    }
  } else {
    // load Jpsi candidates                                                   
    inputArrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
  }

  if(!inputArrayJPSItoEle) {
    printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n");
    return;
  } 

  // fix for temporary bug in ESDfilter 
  // the AODs with null vertex pointer didn't pass the PhysSel
  if(!aod->GetPrimaryVertex()) return;


  // AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  
  // load MC particles
  TClonesArray* mcArray = 
     dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
  if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
     AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));

  // Read AOD Monte Carlo information
  if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);

  // loop over J/Psi->ee candidates
  Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
  printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);

  for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
    AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);

    Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI

    //Secondary vertex
    Double_t vtxSec[3] = {0.,0.,0.};
    Double_t vtxPrim[3] = {0.,0.,0.};
    d->GetSecondaryVtx(vtxSec); 
    Bool_t unsetvtx=kFALSE;
    if(!d->GetOwnPrimaryVtx()) {
      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
      unsetvtx=kTRUE;
    }
    //Int_t okJPSI=0;
    // here analyze only if cuts are passed

    if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
      if (mcLabel == -1)
         {
          AliDebug(2,"No MC particle found");
         }
      else {

      fhInvMass->Fill(d->InvMassJPSIee()); 
      fhD0->Fill(10000*d->ImpParXY());
      fhD0D0->Fill(1e8*d->Prodd0d0());
      fhCosThetaStar->Fill(d->CosThetaStarJPSI());      
      fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0());
      fhCosThetaPointing->Fill(d->CosPointingAngle());

      // compute X variable
      AliAODVertex* primVertex = d->GetOwnPrimaryVtx();
      vtxPrim[0] = primVertex->GetX();
      vtxPrim[1] = primVertex->GetY();
      vtxPrim[2] = primVertex->GetZ();
      Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt();
      Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt();

      fhDecayTime->Fill(10000*pseudoProperDecayTime);
      // Post the data already here
      PostData(1,fOutput);
    
      fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values

      }//
  
     } // end of JPSItoEle candidates selection according to cuts

    if(unsetvtx) d->UnsetOwnPrimaryVtx();

   }// end loop on JPSI to ele candidates

}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/)
{
  //
  // Terminate analysis
  //

  if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n");

  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {
    printf("ERROR: fOutput not available\n");
    return;
  }

  if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
  fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
  fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
  fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
  fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
  fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
  fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
  fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
  fNtupleJPSI = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleJPSI"));

  if(fOkMinimize){
     Double_t* pseudoProperValues=0x0;
     Double_t* invMassValues=0x0;
     Int_t ncand=0;
     fCdfFit = new AliAnalysisBtoJPSItoEle();

     printf("+++\n+++  AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n");
     fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand);
     printf("+++\n+++  X and M vectors filled starting from N-tuple ---> OK\n+++\n");
     printf("+++\n+++  Number of candidates ---> %d J/psi \n+++\n",ncand);
     fCdfFit->SetPtBin(0);
     printf("+++\n+++  Pt bin setted n. ---> %d  \n+++\n",fCdfFit->GetPtBin());

     if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);} 
        else { 
          //SetPathMCfile(".");
          TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin()); 
          fCdfFit->CloneMCtemplate(hLocal); 
        }

     printf("+++\n+++  MC template histo copied ---> OK\n+++\n");
     //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
     fCdfFit->DoMinimization();
   } 

  return;
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
{
  //
  // apply Pt cuts (lower cut, upper cut)
  //
  for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9])
{
  //
  // apply D0 and JPSI cuts 
  //
  for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray) 
{
  //
  // Reads MC information if needed (for example in case of parametrized PID)
  //

  // load MC particles
  TClonesArray* mcArray =
     dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
  if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
     AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));

  // load MC header 
  AliAODMCHeader* mcHeader =
    (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
  if(!mcHeader){
    printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
  }

  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();

  //Double_t rmax = 0.000005;
  Double_t mcPrimVtx[3];

  // get MC primary Vtx
  mcHeader->GetVertex(mcPrimVtx);

  Int_t nInCandidates = inputArray->GetEntriesFast();
  printf("Number of Candidates for MC analysis: %d\n",nInCandidates);

  Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
  Int_t labJPSIdaugh0=0;
  Int_t labJPSIdaugh1=0;

  for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {

    AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
    Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI

    Double_t vtxPrim[3] = {0.,0.,0.};
    Double_t vtxSec[3] = {0.,0.,0.};
    dd->GetSecondaryVtx(vtxSec);
    Bool_t unsetvtx=kFALSE;
    if(!dd->GetOwnPrimaryVtx()) {
      dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
      unsetvtx=kTRUE;
    }

    if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only

      AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
      AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
      lab0 = trk0->GetLabel();
      lab1 = trk1->GetLabel();

      AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
      if(!part0) {
        printf("no MC particle\n");
        continue;
      }

      while(part0->GetMother()>=0) {
       labMother=part0->GetMother();
       part0 = (AliAODMCParticle*)mcArray->At(labMother);
       if(!part0) {
         printf("no MC mother particle\n");
         break;
       }
       pdgMother = TMath::Abs(part0->GetPdgCode());
       if(pdgMother==443) {//this for JPSI
         labJPSIdaugh0=labMother;
         break;
       }
      }

      AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
      if(!part1) {
        printf("no MC particle\n");
        continue;
      }

      while(part1->GetMother()>=0) {
       labMother=part1->GetMother();
       part1 = (AliAODMCParticle*)mcArray->At(labMother);
       if(!part1) {
         printf("no MC mother particle\n");
         break;
       }
       pdgMother = TMath::Abs(part1->GetPdgCode());
       if(pdgMother==443) {//this for JPSI
         labJPSIdaugh1=labMother;
         break;
       }
      }

      if (mcLabelSec == -1)
         {
          AliDebug(2,"No MC particle found");
         }
      else {

      if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
          AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
          pdgJPSI = partJPSI->GetPdgCode();
          Int_t pdaugh0 = partJPSI->GetDaughter(0);
          Int_t pdaugh1 = partJPSI->GetDaughter(1);
          AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
          AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
          pdg0 = partDaugh0->GetPdgCode();
          pdg1 = partDaugh1->GetPdgCode();
        if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
           labJPSIMother = partJPSI->GetMother();
           AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
           pdgJPSIMother = partJPSIMother->GetPdgCode();

        // chech if JPSI comes from B 
           if( pdgJPSIMother==511   || pdgJPSIMother==521   ||
               pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
               pdgJPSIMother==513   || pdgJPSIMother==523   ||
               pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
               pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
               pdgJPSIMother==515   || pdgJPSIMother==525   ||
               pdgJPSIMother==531   || pdgJPSIMother==10531 ||
               pdgJPSIMother==533   || pdgJPSIMother==10533 ||
               pdgJPSIMother==20533 || pdgJPSIMother==535   ||
               pdgJPSIMother==541   || pdgJPSIMother==10541 ||
               pdgJPSIMother==543   || pdgJPSIMother==10543 ||
               pdgJPSIMother==20543 || pdgJPSIMother==545)
               { //this is for MC JPSI -> ee from B-hadrons

                  fhInvMass->Fill(dd->InvMassJPSIee()); 
                  fhD0->Fill(10000*dd->ImpParXY());
                  fhD0D0->Fill(1e8*dd->Prodd0d0());
                  fhCosThetaStar->Fill(dd->CosThetaStarJPSI());      
                  fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
                  fhCosThetaPointing->Fill(dd->CosPointingAngle());
  
                  // compute X variable
                  AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
                  vtxPrim[0] = primVertex->GetX();
                  vtxPrim[1] = primVertex->GetY();
                  vtxPrim[2] = primVertex->GetZ();
                  Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
                  Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
  
                  fhDecayTime->Fill(10000*pseudoProperDecayTime);
                  // Post the data already here
                  PostData(1,fOutput);
    
                  fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values

                  Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
                  Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
                  fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);

                  // Post the data already here
                  PostData(1,fOutput);

               } //this is for MC JPSI -> ee from B-hadrons
            } //this is for MC JPSI -> ee
        }
       } // end of check if JPSI is signal
      }
    }

}
//_________________________________________________________________________________________________
TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum) 
{
  //
  // Open a local file with MC x distribution for JPSI from B-hadron
  //

  TH1F* hCsiMC = new TH1F();
  TString useFile = datadir.Data();
  useFile.Append("/CsiMCfromKine_10PtBins.root");
  TFile *f = new TFile(useFile);
  Double_t scale = 0.;
  char processCsiMCXHisto[1024];
  if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper");
  else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum);
  hCsiMC = (TH1F*)f->Get(processCsiMCXHisto);
  scale=1/hCsiMC->Integral();
  hCsiMC->Scale(scale);
  printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries());

  return hCsiMC;

}
Commit	Line	Data
7b54a483	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
b557eb43	31	#include "AliAnalysisManager.h"
b557eb43	32	#include "AliAODHandler.h"
7b54a483	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
7b54a483	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	{
9d6d35b0	81	// standard constructor
7b54a483	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
7b54a483	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
f9fd1412	124	fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
7b54a483	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
b557eb43	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	}
8931c313	191
7b54a483	192	if(!inputArrayJPSItoEle) {
	193	printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n");
	194	return;
	195	}
	196
7c23877d	197	// fix for temporary bug in ESDfilter
	198	// the AODs with null vertex pointer didn't pass the PhysSel
	199	if(!aod->GetPrimaryVertex()) return;
	200
	201
b557eb43	202	// AOD primary vertex
	203	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
	204
f9fd1412	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
7b54a483	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();
f9fd1412	216	printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);
7b54a483	217
	218	for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
	219	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
f9fd1412	220
	221	Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
	222
7b54a483	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
f9fd1412	234
7b54a483	235	if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
f9fd1412	236	if (mcLabel == -1)
	237	{
	238	AliDebug(2,"No MC particle found");
	239	}
	240	else {
	241
7b54a483	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
f9fd1412	263	}//
f9fd1412	264
7b54a483	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");
f9fd1412	318	//fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
f9fd1412	319	fCdfFit->DoMinimization();
7b54a483	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	//_________________________________________________________________________________________________
f9fd1412	341	void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
7b54a483	342	{
	343	//
	344	// Reads MC information if needed (for example in case of parametrized PID)
	345	//
	346
	347	// load MC particles
f9fd1412	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()));
7b54a483	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
f9fd1412	360	AliAODVertex vtx1 = (AliAODVertex)aodEvent->GetPrimaryVertex();
	361
	362	//Double_t rmax = 0.000005;
7b54a483	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++) {
f9fd1412	376
7b54a483	377	AliAODRecoDecayHF2Prong dd = (AliAODRecoDecayHF2Prong)inputArray->UncheckedAt(iCandidate);
f9fd1412	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
7b54a483	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
f9fd1412	436	if (mcLabelSec == -1)
	437	{
	438	AliDebug(2,"No MC particle found");
	439	}
	440	else {
	441
7b54a483	442	if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
f9fd1412	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
7b54a483	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
f9fd1412	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
7b54a483	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