Added method to recalculate the primary vertex without the daughter tracks

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

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

}