Add efficiency histos in the task for combinatorics (Francesco) + updated macro for...
authorfprino <prino@to.infn.it>
Fri, 28 Feb 2014 23:32:37 +0000 (00:32 +0100)
committerfprino <prino@to.infn.it>
Fri, 28 Feb 2014 23:32:37 +0000 (00:32 +0100)
PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h
PWGHF/vertexingHF/macros/AddTaskCombinHF.C

index 8ba5591..be79632 100644 (file)
@@ -24,6 +24,7 @@
 
 #include <TList.h>
 #include <TH1F.h>
+#include <TH2F.h>
 #include <TH3F.h>
 #include <THnSparse.h>
 
@@ -36,6 +37,7 @@
 #include "AliAODMCHeader.h"
 #include "AliAODVertex.h"
 #include "AliAODTrack.h"
+#include "AliVertexingHFUtils.h"
 #include "AliAnalysisTaskCombinHF.h"
 
 ClassImp(AliAnalysisTaskCombinHF)
@@ -47,6 +49,14 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
   fOutput(0x0), 
   fHistNEvents(0x0),
   fHistTrackStatus(0x0),
+  fHistCheckOrigin(0x0),
+  fHistCheckOriginSel(0x0),
+  fHistCheckDecChan(0x0),
+  fHistCheckDecChanAcc(0x0),
+  fPtVsYGen(0x0),
+  fPtVsYGenLimAcc(0x0),
+  fPtVsYGenAcc(0x0),
+  fPtVsYReco(0x0),
   fMassVsPtVsY(0x0),
   fMassVsPtVsYRot(0x0),
   fMassVsPtVsYLSpp(0x0),
@@ -66,6 +76,8 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
   fAnalysisCuts(0x0),
   fMinMass(1.750),
   fMaxMass(2.150),
+  fEtaAccCut(0.9),
+  fPtAccCut(0.1),
   fNRotations(9),
   fMinAngleForRot(5*TMath::Pi()/6),
   fMaxAngleForRot(7*TMath::Pi()/6),
@@ -74,6 +86,8 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
   fCounter(0x0),
   fMeson(kDzero),
   fReadMC(kFALSE),
+  fPromptFeeddown(kPrompt),
+  fGoUpToQuark(kTRUE),
   fFullAnalysis(0)  
 {
   // default constructor
@@ -85,6 +99,14 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy
   fOutput(0x0), 
   fHistNEvents(0x0),
   fHistTrackStatus(0x0),
+  fHistCheckOrigin(0x0),
+  fHistCheckOriginSel(0x0),
+  fHistCheckDecChan(0x0),
+  fHistCheckDecChanAcc(0x0),
+  fPtVsYGen(0x0),
+  fPtVsYGenLimAcc(0x0),
+  fPtVsYGenAcc(0x0),
+  fPtVsYReco(0x0),
   fMassVsPtVsY(0x0),
   fMassVsPtVsYRot(0x0),
   fMassVsPtVsYLSpp(0x0),
@@ -104,6 +126,8 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy
   fAnalysisCuts(analysiscuts),
   fMinMass(1.750),
   fMaxMass(2.150),
+  fEtaAccCut(0.9),
+  fPtAccCut(0.1),
   fNRotations(9),
   fMinAngleForRot(5*TMath::Pi()/6),
   fMaxAngleForRot(7*TMath::Pi()/6),
@@ -112,6 +136,8 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy
   fCounter(0x0),
   fMeson(meson),
   fReadMC(kFALSE),
+  fPromptFeeddown(1),
+  fGoUpToQuark(kTRUE),
   fFullAnalysis(0)
 
 {
@@ -130,6 +156,14 @@ AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
   delete fOutput;
   delete fHistNEvents;
   delete fHistTrackStatus;
+  delete fHistCheckOrigin;
+  delete fHistCheckOriginSel;
+  delete fHistCheckDecChan;
+  delete fHistCheckDecChanAcc;
+  delete fPtVsYGen;
+  delete fPtVsYGenLimAcc;
+  delete fPtVsYGenAcc;
+  delete fPtVsYReco;
   delete fMassVsPtVsY; 
   delete fMassVsPtVsYLSpp;
   delete fMassVsPtVsYLSmm;
@@ -188,6 +222,50 @@ void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
   fHistTrackStatus->SetMinimum(0);
   fOutput->Add(fHistTrackStatus);
 
+  if(fReadMC){
+
+    fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
+    fHistCheckOrigin->Sumw2();
+    fHistCheckOrigin->SetMinimum(0);
+    fOutput->Add(fHistCheckOrigin);
+
+    fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
+    fHistCheckOriginSel->Sumw2();
+    fHistCheckOriginSel->SetMinimum(0);
+    fOutput->Add(fHistCheckOriginSel);
+
+    fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
+    fHistCheckDecChan->Sumw2();
+    fHistCheckDecChan->SetMinimum(0);
+    fOutput->Add(fHistCheckDecChan);
+
+    fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
+    fHistCheckDecChanAcc->Sumw2();
+    fHistCheckDecChanAcc->SetMinimum(0);
+    fOutput->Add(fHistCheckDecChanAcc);
+
+    fPtVsYGen= new TH2F("hPtVsYGen","",20,0.,10.,20,-1.,1.);
+    fPtVsYGen->Sumw2();
+    fPtVsYGen->SetMinimum(0);
+    fOutput->Add(fPtVsYGen);
+
+    fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",20,0.,10.,20,-1.,1.);
+    fPtVsYGenLimAcc->Sumw2();
+    fPtVsYGenLimAcc->SetMinimum(0);
+    fOutput->Add(fPtVsYGenLimAcc);
+
+    fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",20,0.,10.,20,-1.,1.);
+    fPtVsYGenAcc->Sumw2();
+    fPtVsYGenAcc->SetMinimum(0);
+    fOutput->Add(fPtVsYGenAcc);
+
+    fPtVsYReco= new TH2F("hPtVsYReco","",20,0.,10.,20,-1.,1.);
+    fPtVsYReco->Sumw2();
+    fPtVsYReco->SetMinimum(0);
+    fOutput->Add(fPtVsYReco);
+  }
+
+
  Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
   Double_t maxm=fMinMass+nMassBins*0.001;
   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
@@ -312,6 +390,7 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
       return;
     }
+    FillGenHistos(arrayMC);
   }
 
 
@@ -440,6 +519,52 @@ void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD
 }
 
 //________________________________________________________________________
+void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
+  // Fill histos with generated quantities
+  Int_t totPart=arrayMC->GetEntriesFast();
+  Int_t thePDG=411;
+  Int_t nProng=3;
+  if(fMeson==kDzero){
+    thePDG=421;
+    nProng=2;
+  }
+  Int_t labDau[4];
+  for(Int_t ip=0; ip<totPart; ip++){
+    AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
+    if(TMath::Abs(part->GetPdgCode())==thePDG){
+      Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
+      fHistCheckOrigin->Fill(orig);
+      if(fPromptFeeddown==kFeeddown && orig!=5) continue;
+      else if(fPromptFeeddown==kPrompt && orig!=4) continue;
+      else if(fPromptFeeddown==kBoth && orig<4) continue;
+      fHistCheckOriginSel->Fill(orig);
+      Int_t deca=0;
+      Bool_t isGoodDecay=kFALSE;
+      if(fMeson==kDzero){
+       deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
+       if(part->GetNDaughters()!=2) continue;
+       if(deca==1) isGoodDecay=kTRUE;
+      }else if(fMeson==kDplus){ 
+       deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
+       if(deca>0) isGoodDecay=kTRUE;
+      }
+      Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
+      fHistCheckDecChan->Fill(deca);
+      if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
+      if(isGoodDecay){
+       Double_t ptgen=part->Pt();
+       Double_t ygen=part->Y();
+       if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
+         fPtVsYGen->Fill(ptgen,ygen);
+         if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
+         if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
+       }
+      }
+    }
+  }
+}
+
+//________________________________________________________________________
 Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels){
   // Fill histos for candidates with proper charge sign
 
@@ -460,11 +585,16 @@ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD
        for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
        Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
        if(labD>=0){
-         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(dgLabels[0]));
+         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
          if(part){
            Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
-           if(pdgCode==321) fMassVsPtVsYSig->Fill(mass,pt,rapid);
-           else fMassVsPtVsYRefl->Fill(mass,pt,rapid);
+           if(pdgCode==321){
+             fMassVsPtVsYSig->Fill(mass,pt,rapid);
+             AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
+             fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
+           }else{
+             fMassVsPtVsYRefl->Fill(mass,pt,rapid);
+           }
          }
        }else{
          fMassVsPtVsYBkg->Fill(mass,pt,rapid);
@@ -571,6 +701,19 @@ Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCu
 }
 
 //_________________________________________________________________
+Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
+  // check if the decay products are in the good eta and pt range
+  for (Int_t iProng = 0; iProng<nProng; iProng++){
+    AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
+    if(!mcPartDaughter) return kFALSE;
+    Double_t eta = mcPartDaughter->Eta();
+    Double_t pt = mcPartDaughter->Pt();
+    if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
+  }
+  return kTRUE;
+}
+
+//_________________________________________________________________
 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
 {
   // Terminate analysis
index f1150e7..7069d43 100644 (file)
@@ -36,6 +36,11 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE
   virtual void Terminate(Option_t *option);
 
   void SetReadMC(Bool_t read){fReadMC=read;}
+  void SelectPromptD(){fPromptFeeddown=kPrompt;}
+  void SelectFeeddownD(){fPromptFeeddown=kFeeddown;}
+  void SelectPromptAndFeeddownD(){fPromptFeeddown=kBoth;}
+  void SetGoUpToQuark(Bool_t opt){fGoUpToQuark=opt;}
+
   void SetTrackCuts(AliESDtrackCuts* cuts){
     if(fTrackCutsAll) delete fTrackCutsAll;
     fTrackCutsAll=new AliESDtrackCuts(*cuts);
@@ -68,8 +73,10 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE
   Bool_t SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts);
   Bool_t FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels);
   void FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge);
-
+  void FillGenHistos(TClonesArray* arrayMC);
+  Bool_t CheckAcceptance(TClonesArray* arrayMC, Int_t nProng, Int_t *labDau);
   enum EMesonSpecies {kDzero, kDplus, kDstar, kDs};
+  enum EPrompFd {kNone,kPrompt,kFeeddown,kBoth};
 
  private:
 
@@ -77,10 +84,17 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE
   AliAnalysisTaskCombinHF& operator=(const AliAnalysisTaskCombinHF& source); 
 
   TList   *fOutput; //! list send on output slot 0
-  TH1F *fHistNEvents; //!hist. for No. of events
-
-  TH1F *fHistTrackStatus; //!hist. of status of tracks
-  TH3F *fMassVsPtVsY;   //! hist. of Y vs. Pt vs. Mass (all cand)
+  TH1F *fHistNEvents;         //!hist. for No. of events
+  TH1F *fHistTrackStatus;     //!hist. of status of tracks
+  TH1F *fHistCheckOrigin;     //!hist. of origin (c/b) of D meson
+  TH1F *fHistCheckOriginSel;  //!hist. of origin (c/b) of D meson
+  TH1F *fHistCheckDecChan;    //!hist. of decay channel of D meson
+  TH1F *fHistCheckDecChanAcc; //!hist. of decay channel of D meson in acc.
+  TH2F *fPtVsYGen;        //! hist. of Y vs. Pt generated (all D)
+  TH2F *fPtVsYGenLimAcc;  //! hist. of Y vs. Pt generated (|y|<0.5)
+  TH2F *fPtVsYGenAcc;     //! hist. of Y vs. Pt generated (D in acc)
+  TH2F *fPtVsYReco;       //! hist. of Y vs. Pt generated (Reco D)
+  TH3F *fMassVsPtVsY;     //! hist. of Y vs. Pt vs. Mass (all cand)
   TH3F *fMassVsPtVsYRot;   //! hist. of Y vs. Pt vs. Mass (rotations)
   TH3F *fMassVsPtVsYLSpp;  //! hist. of Y vs. Pt vs. Mass (like sign ++)
   TH3F *fMassVsPtVsYLSmm;  //! hist. of Y vs. Pt vs. Mass (like sign --)
@@ -101,6 +115,10 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE
   Double_t fMinMass; // minimum value of invariant mass
   Double_t fMaxMass; // maximum value of invariant mass
 
+  Double_t fEtaAccCut; // eta limits for acceptance step
+  Double_t fPtAccCut; // pt limits for acceptance step
+
+
   Int_t fNRotations; // number of rotations
   Double_t fMinAngleForRot; // minimum angle for track rotation
   Double_t fMaxAngleForRot; // maximum angle for track rotation
@@ -108,11 +126,14 @@ class AliAnalysisTaskCombinHF : public AliAnalysisTaskSE
   Double_t fMaxAngleForRot3; // maximum angle for track rotation (3rd prong)
 
   AliNormalizationCounter *fCounter;//!Counter for normalization
+
   Int_t fMeson;          // mesonSpecies (see enum)
-  Bool_t  fReadMC;                    //  flag for access to MC
-  Int_t fFullAnalysis; // flag to set analysis level (0 is the fastest)
+  Bool_t  fReadMC;       //  flag for access to MC
+  Int_t fPromptFeeddown; // flag to select prompt (1), feeddown (2) or all (3) 
+  Bool_t fGoUpToQuark;   // flag for definition of c,b origin
+  Int_t fFullAnalysis;   // flag to set analysis level (0 is the fastest)
 
-  ClassDef(AliAnalysisTaskCombinHF,1); // D+ task from AOD tracks
+  ClassDef(AliAnalysisTaskCombinHF,2); // D+ task from AOD tracks
 };
 
 #endif
index 547492a..7aebe5b 100644 (file)
@@ -1,4 +1,4 @@
-AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, Bool_t readMC=kTRUE)
+AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, TString containerStr="",Bool_t readMC=kTRUE, TString cutObjFile="",TString cutObjNam="")
 {
 
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
@@ -8,23 +8,41 @@ AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, Bool_t readMC=kTRUE)
 
   //Analysis Task
 
-  AliRDHFCuts* analysiscuts;
-  if(meson==0) analysiscuts=new AliRDHFCutsD0toKpi();
-  else analysiscuts=new AliRDHFCutsDplustoKpipi();
-  analysiscuts->SetStandardCutsPP2010();
-  
+
+  AliRDHFCuts* analysiscuts=0x0;
+  AliAODPidHF* pid=0x0;
+  if(!cutObjFile.IsNull()){
+    TFile *f=TFile::Open(cutObjFile.Data(),"READ");
+    if(f){
+      analysiscuts=(AliRDHFCuts*)f->Get(cutObjNam.Data());
+      pid=analysiscuts->GetPidHF();
+    }
+  }
+  else {
+    if(meson==0) analysiscuts=new AliRDHFCutsD0toKpi();
+    else analysiscuts=new AliRDHFCutsDplustoKpipi();
+    analysiscuts->SetStandardCutsPP2010();
+    pid=new AliAODPidHF();
+    pid->SetMatch(5);
+    pid->SetTPCnSigmaRangeForPions(-3.,3.);
+    pid->SetTPCnSigmaRangeForKaons(-2.,3.);
+    pid->SetTPCnSigmaRangeForProtons(-3.,3.);
+    pid->SetTOFnSigmaRangeForPions(-3.,3.);
+    pid->SetTOFnSigmaRangeForKaons(-2.,2.);
+    pid->SetTOFnSigmaRangeForProtons(-3.,3.);
+
+  }
+  if(!analysiscuts){
+    Printf("Wrong file or cut object name set");
+    return 0x0;
+  }
+
   AliAnalysisTaskCombinHF *dTask = new AliAnalysisTaskCombinHF(meson,analysiscuts);
   dTask->SetReadMC(readMC);
   dTask->SetDebugLevel(0);
-  AliAODPidHF* pid=new AliAODPidHF();
-  pid->SetMatch(5);
-  pid->SetTPCnSigmaRangeForPions(-3.,3.);
-  pid->SetTPCnSigmaRangeForKaons(-2.,3.);
-  pid->SetTPCnSigmaRangeForProtons(-3.,3.);
-  pid->SetTOFnSigmaRangeForPions(-3.,3.);
-  pid->SetTOFnSigmaRangeForKaons(-2.,2.);
-  pid->SetTOFnSigmaRangeForProtons(-3.,3.);
   dTask->SetPIDHF(pid);
+
+
   mgr->AddTask(dTask);
   
   // Create containers for input/output 
@@ -32,13 +50,13 @@ AliAnalysisTaskCombinHF *AddTaskCombinHF(Int_t meson=0, Bool_t readMC=kTRUE)
   TString mesname="Dzero";
   if(meson==1) mesname="Dplus";
   TString inname = Form("cinput%s",mesname.Data());
-  TString outname = Form("coutput%s",mesname.Data());
-  TString normname = Form("coutput%sNorm",mesname.Data());
+  TString outname = Form("coutput%s%s",mesname.Data(),containerStr.Data());
+  TString normname = Form("coutput%sNorm%s",mesname.Data(),containerStr.Data());
 
   AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),
                                                          AliAnalysisManager::kInputContainer);
   TString outputfile = AliAnalysisManager::GetCommonFileName();
-  outputfile += Form(":PWG3_D2H_InvMass%sLowPt",mesname.Data());
+  outputfile += Form(":PWG3_D2H_InvMass%sLowPt%s",mesname.Data(),containerStr.Data());
   
   
   AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),