Faster calculation of invariant mass (Yifei)

[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEJPSItoEle.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 AliAnalysisTaskSEJPSItoEle
// AliAnalysisTaskSE class to read both J/psi -> e+e- candidates and 
// like-sign pair candidates and to store them into a specific
// stand-alone AOD file for J/psi into dieletrons analysis. 
// The current task:
//
//  1) produces several histograms (including invariant mass distributions)
//     for both unlike sign (US) and like sign (LS) pairs.
//
//  2) selects only J/Psi to dielectrons candidates and copies it to a 
//     specific AOD file, namely AliAODjpsi.root. The references
//     to AliAODTrack objects are copied as well. 
//     The specific AliAODjpsi.root is thought  as the input file 
//     to the AliAnalysisBtoJPSItoEle.h class, which performs (locally)
//     the analysis of J/Psi from beauty hadrons decays. 
//     
//     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 "TROOT.h"
#include "TCanvas.h"
#include "TBits.h" //**********************

#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODHandler.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskSEJPSItoEle.h"
#include "AliAODTrack.h"
#include "AliExternalTrackParam.h"

ClassImp(AliAnalysisTaskSEJPSItoEle)

//______________________________________________________________________________
AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle() :
AliAnalysisTaskSE(),
fOutput(0),
fhDecayTimeMCjpsifromB(0),
fhDecayTime(0),                         
fhDecayTimeOut(0),
fhInvMass(0),                           
fhdEdxTPC(0),
fhD0(0),                                
fhD0D0(0),                              
fhCosThetaStar(0),                      
fhCosThetaPointing(0),                  
fhDCA(0),
fhCtsVsD0D0(0),
fHistMassLS(0),
fHistCtsLS(0),
fHistCtsLSpos(0),
fHistCtsLSneg(0),
fHistCPtaLS(0),
fHistd0d0LS(0),
fHistDCALS(0),
fVHF(0),
fTotPosPairs(0),
fTotNegPairs(0),
fLsNormalization(1.),
fOkAODMC(kFALSE),
fOkLikeSign(kFALSE),
fVerticesHFTClArr(0),
fJpsiToEleTClArr(0),
fLikeSignTClArr(0),
fTracksTClArr(0),
fOrigAOD(0),
fNewAOD(0)
{
  // default constructor
}                 
//_________________________________________________________________________________________________
AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle(const char* name) :
AliAnalysisTaskSE(name),
fOutput(0),
fhDecayTimeMCjpsifromB(0),
fhDecayTime(0),
fhDecayTimeOut(0),
fhInvMass(0),                           
fhdEdxTPC(0),
fhD0(0),                                
fhD0D0(0),                              
fhCosThetaStar(0),                      
fhCosThetaPointing(0),                  
fhDCA(0),
fhCtsVsD0D0(0),
fHistMassLS(0),
fHistCtsLS(0),
fHistCtsLSpos(0),
fHistCtsLSneg(0),
fHistCPtaLS(0),
fHistd0d0LS(0),
fHistDCALS(0),
fVHF(0),
fTotPosPairs(0),
fTotNegPairs(0),
fLsNormalization(1.),
fOkAODMC(kFALSE),
fOkLikeSign(kFALSE),
fVerticesHFTClArr(0),
fJpsiToEleTClArr(0),
fLikeSignTClArr(0),
fTracksTClArr(0),
fOrigAOD(0),
fNewAOD(0)
{
  // standard constructor

  // Input slot #0 works with an Ntuple
  DefineInput(0, TChain::Class());

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

  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  } 

  if (fVHF) {
    delete fVHF;
    fVHF = 0;
  }

}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEJPSItoEle::Init()
{
  // Initialization

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

  gROOT->LoadMacro("ConfigVertexingHF.C");

  fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
  fVHF->PrintStatus();

  return;
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() 
{
  // Create the output container

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

  if (!AODEvent()) {
     Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
     return;
  }   

  if(!fNewAOD) fNewAOD = new AliAODEvent();
  fNewAOD->CreateStdContent();

  TString filename = "AliAOD.Jpsitoele.root";
  if (!IsStandardAOD()) filename = "";

  // Array of HF vertices  
  fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
  fVerticesHFTClArr->SetName("VerticesHF");
  AddAODBranch("TClonesArray", &fVerticesHFTClArr);

  // Array of J/psi -> e+e- candidates 
  fJpsiToEleTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
  fJpsiToEleTClArr->SetName("JpsiToEleEvents");
  AddAODBranch("TClonesArray", &fJpsiToEleTClArr, filename);

  // Array of like sign candidates 
  fLikeSignTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
  fLikeSignTClArr->SetName("LikeSignEvents");
  AddAODBranch("TClonesArray", &fLikeSignTClArr, filename);

  // Array of tracks from J/psi -> e+e- candidates
  fTracksTClArr = new TClonesArray("AliAODTrack", 0);
  fTracksTClArr->SetName("Tracks");
  AddAODBranch("TClonesArray", &fTracksTClArr, filename);

  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);
  }

  // Invariant mass
  fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25);
  fhInvMass->Sumw2();
  fhInvMass->SetMinimum(0);
  fOutput->Add(fhInvMass);

  fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25);
  fHistMassLS->Sumw2();
  fHistMassLS->SetMinimum(0);
  fOutput->Add(fHistMassLS);

  // TPC signal (only for candidate tracks)
  fhdEdxTPC = new TH2F("fhdEdxTPC","TPC signal (dE/dx)",100,0.,10.,1000,0.,100.);
  fhdEdxTPC->SetXTitle("p_{T} [GeV]");
  fhdEdxTPC->SetYTitle("TPC signal (arb. units)");
  fOutput->Add(fhdEdxTPC);

  // Pseudo proper decay time
  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);
 
  // Pseudo proper decay time
  fhDecayTimeOut = new TH1F("fhDecayTimeOut", "J/#Psi pseudo proper decay time (output standalone AOD); X [#mu m]; Entries",200,-2000.,2000.);
  fhDecayTimeOut->Sumw2();
  fhDecayTimeOut->SetMinimum(0);
  fOutput->Add(fhDecayTimeOut);

  // Dictance of closest approach
  fhDCA = new TH1F("fhDCA", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
  fhDCA->Sumw2();
  fhDCA->SetMinimum(0);
  fOutput->Add(fhDCA);

  fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
  fHistDCALS->Sumw2();
  fHistDCALS->SetMinimum(0);
  fOutput->Add(fHistDCALS);

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

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

  fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
  fHistd0d0LS->Sumw2();
  fHistd0d0LS->SetMinimum(0);
  fOutput->Add(fHistd0d0LS);

  // Cosine of the decay angle
  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);

  fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
  fHistCtsLS->Sumw2();
  fHistCtsLS->SetMinimum(0);
  fOutput->Add(fHistCtsLS);

  fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
  fHistCtsLSpos->Sumw2();
  fHistCtsLSpos->SetMinimum(0);
  fOutput->Add(fHistCtsLSpos);

  fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
  fHistCtsLSneg->Sumw2();
  fHistCtsLSneg->SetMinimum(0);
  fOutput->Add(fHistCtsLSneg);

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

  fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.);
  fHistCPtaLS->Sumw2();
  fHistCPtaLS->SetMinimum(0);
  fOutput->Add(fHistCPtaLS);

  // CosThetaStar vs. d0xd0
  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);

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

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

  TClonesArray *arrayJPSItoEle = 0;
  TClonesArray *arrayLikeSign = 0;
  TClonesArray *arrayTracks = 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 tracks from AOD event   
      arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks");
      Int_t totTracks = arrayTracks->GetEntriesFast();
      if(fDebug>1) printf("Number of tracks === %d \n",totTracks);

      // load Jpsi candidates from AOD friend  
      arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
      Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
      if(fDebug>1) printf("Number of J/psi->ee candidates  === %d \n",totJPSItoEleCand);
  
      // load like sign candidates from AOD friend
      arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong");
      Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
      if(fDebug>1) printf("Number of like sign candidates  === %d \n",totLikeSignCand);
    }

  } else {
    // load tracks from AOD event
    arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks");
    Int_t totTracks = arrayTracks->GetEntriesFast();
    if(fDebug>1) printf("Number of tracks === %d \n",totTracks);

    // load Jpsi candidates                                                   
    arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
    Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
    if(fDebug>1) printf("Number of J/psi->ee candidates  === %d \n",totJPSItoEleCand);

    // load like sign candidates
    arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
    Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
    if(fDebug>1) printf("Number of like sign candidates  === %d \n",totLikeSignCand);
  }

  fOrigAOD = aod; // copy pointer to the current AliAODEvent in the data member fOrigAOD
  if (!aod) return;

  if(!arrayTracks) {
    printf("AliAnalysisTaskSEJPSItoEle::UserExec: Tracks branch not found!\n");
    return;
  }
  if(!arrayJPSItoEle) {
    printf("AliAnalysisTaskSEJPSItoEle::UserExec: JPSItoEle branch not found!\n");
    return;
  }
  if(!arrayLikeSign) {
    printf("AliAnalysisTaskSEJPSItoEle::UserExec: LikeSign2Prong branch not found!\n");
    return;
  }

  // load MC particles and read MC info (for sim only)
  TClonesArray* mcArray=0;
  if(fOkAODMC){
     mcArray =  dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
     if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
     ReadAODMCInfo(aod,arrayJPSItoEle);
     }

  // retrieve AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();

  Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0; 

  fVerticesHFTClArr->Delete();
  iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
  TClonesArray &verticesHFRef = *fVerticesHFTClArr;

  fJpsiToEleTClArr->Delete();
  iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast();
  TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr;

  fLikeSignTClArr->Delete();
  iOutLikeSign = fLikeSignTClArr->GetEntriesFast();
  TClonesArray &arrayLikeSignRef = *fLikeSignTClArr;

  fTracksTClArr->Delete();
  iOutTracks = fTracksTClArr->GetEntriesFast();
  //TClonesArray &arrayTrackRef = *fTracksTClArr;

  //
  // LOOP OVER LIKE SIGN CANDIDATES
  //

  // Access to the new AOD container of tracks

  Int_t trkIDtoEntry[100000];
  AliAODTrack* trackin=0;
  for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
    trackin = aod->GetTrack(it);
    trkIDtoEntry[trackin->GetID()]=it;
    }
/*
  for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
    vtx1->AddDaughter(trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack((*(aod->GetTrack(trkIDtoEntry[it])))));
    trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
    printf("trackin ==== %d \n", trackin);
    AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
    }
  const AliAODTrack* trackin;  
  Int_t numTracks = arrayTracks->GetEntriesFast();
  for(Int_t ntracks=0; ntracks<numTracks; ntracks++){
  trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
  printf("trackin ==== %d \n", trackin);
    AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
  }
*/

  Int_t nPosPairs=0,nNegPairs=0;
  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
  AliAODRecoDecayHF2Prong* dlsout = 0;
  //if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign);

  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
    AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
    Bool_t unsetvtx=kFALSE;
    if(!dlsin->GetOwnPrimaryVtx()) {
        dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
        unsetvtx=kTRUE;
    }
    //Int_t okBtoJPSIls=0;
    //if(dlsin->Pt() < fPtCuts[1] && dlsin->Pt() > fPtCuts[0]){ // apply pt cut only
    //if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) {
       AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0);
       fHistMassLS->Fill(dlsin->InvMassJPSIee());
       fHistCPtaLS->Fill(dlsin->CosPointingAngle());
       fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0());
       fHistDCALS->Fill(100*dlsin->GetDCA());
       if(!trk0) {
          trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]);
          printf("references to standard AOD not available \n");
       }
       if((trk0->Charge())==1) {
          nPosPairs++;
          fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11));
          fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11));
        } else {
          nNegPairs++;
          fHistCtsLS->Fill(dlsin->CosThetaStarJPSI());
          fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI());
        }
       PostData(1,fOutput);
     
    // Clone like sign candidate for output AOD
    dlsout = new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin);

    //}
    if(unsetvtx) dlsin->UnsetOwnPrimaryVtx();
  }

  //if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs);
  //if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs);

  fTotPosPairs += nPosPairs;
  fTotNegPairs += nNegPairs;

  //
  // LOOP OVER J/psi->e+e- CANDIDATES
  //

  Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast();
  //if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle);

  //totJPSIin +=  nInJPSItoEle;
  Bool_t isOkEle[2] = {kFALSE,kFALSE};

  for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {

    AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle);
    Int_t mcLabel = 0;
    if(fOkAODMC) mcLabel = dIn->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.};
    Double_t vtxSecOut[3] = {0.,0.,0.};
    Double_t vtxPrimOut[3] = {0.,0.,0.};
    dIn->GetSecondaryVtx(vtxSec); 
    Bool_t unsetvtx=kFALSE;
    if(!dIn->GetOwnPrimaryVtx()) {
      dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
      unsetvtx=kTRUE;
    }

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

         //fhInvMass->Fill(dIn->InvMassJPSIee()); 
         fhD0->Fill(10000*dIn->ImpParXY());
         fhD0D0->Fill(1e8*dIn->Prodd0d0());
         fhCosThetaStar->Fill(dIn->CosThetaStarJPSI());      
         fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0());
         fhCosThetaPointing->Fill(dIn->CosPointingAngle());
         fhDCA->Fill(100*dIn->GetDCA());

         // compute X variable
         AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx();
         vtxPrim[0] = primVertex->GetX();
         vtxPrim[1] = primVertex->GetY();
         vtxPrim[2] = primVertex->GetZ();
         Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt();
         Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt();
         fhDecayTime->Fill(10000*pseudoProperDecayTime);
  
         // retrieve daughter tracks 
         AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
         TBits clusterMap0 = trk0->GetTPCClusterMap();
         Int_t npoints0 = clusterMap0.CountBits(0); 
         AliAODPid *pid0 = trk0->GetDetPid();
         if(GetNumberOfSigmas(trk0->P(),pid0->GetTPCsignal(),npoints0,AliPID::kElectron) < 3.) isOkEle[0]=kTRUE;

         AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
         TBits clusterMap1 = trk1->GetTPCClusterMap();
         Int_t npoints1 = clusterMap1.CountBits(0);
         AliAODPid *pid1 = trk1->GetDetPid();
         if(GetNumberOfSigmas(trk1->P(),pid1->GetTPCsignal(),npoints1,AliPID::kElectron) < 3.) isOkEle[1]=kTRUE;

         if(isOkEle[0] && isOkEle[1]) {
           fhInvMass->Fill(dIn->InvMassJPSIee()); 
           fhdEdxTPC->Fill(pid0->GetTPCmomentum(),pid0->GetTPCsignal());
           fhdEdxTPC->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
         } 

         //printf("TPC momentum %f\n", pid0->GetTPCmomentum());
         //printf("TPC signal %f\n", pid0->GetTPCsignal());

         // clone candidate for output AOD
         AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++]) 
           AliAODVertex(*(dIn->GetSecondaryVtx()));
         AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++]) 
           AliAODRecoDecayHF2Prong(*dIn); // copy constructor used
         dOut->SetSecondaryVtx(v);
         dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
         v->SetParent(dOut);

         // compute X variable from the cloned candidate in the stand-alone AOD file
         dOut->GetSecondaryVtx(vtxSecOut);
         AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx();
         vtxPrimOut[0] = primVertexOut->GetX();
         vtxPrimOut[1] = primVertexOut->GetY();
         vtxPrimOut[2] = primVertexOut->GetZ();
         Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt();
         Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt();
         fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut);

         //AliAODTrack*  trk0 = (AliAODTrack*)dOut->GetDaughter(0);
         //printf("+++ Pt first daughter = %d \n +++ \n", trk0->Pt());

        }

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

    if(unsetvtx) dIn->UnsetOwnPrimaryVtx();

    PostData(0,fOutput);

 }// end loop on JPSI to ele candidates

 //printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle);
 
 //totJPSIout += iOutJPSItoEle;

 return;

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

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

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

  fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);

  if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
  fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
  fhDecayTimeOut = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeOut"));
  fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
  fhdEdxTPC = dynamic_cast<TH2F*>(fOutput->FindObject("fhdEdxTPC"));
  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"));
  fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
  fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
  fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
  fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
  fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
  fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
  fhDCA = dynamic_cast<TH1F*>(fOutput->FindObject("fhDCA"));
  fHistDCALS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));

  if(fLsNormalization>0.) {
    fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
    fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
    fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
    fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
    fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
    fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
    fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
  }

  //printf("Tot JPSI in %d\n", totJPSIin);
  //printf("Tot JPSI out %d\n", totJPSIout);

  return;
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEJPSItoEle::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 AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9])
{
  //
  // apply JPSI cuts 
  //
  for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
}
//_________________________________________________________________________________________________
void AliAnalysisTaskSEJPSItoEle::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("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n");
  }

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

  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);

                  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
    }
  }
}