added possibility to use MC centrality from generated tracks (Tim)
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Sep 2013 06:53:08 +0000 (06:53 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Sep 2013 06:53:08 +0000 (06:53 +0000)
PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg
PWGCF/Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
PWGCF/Correlations/DPhi/FourierDecomposition/AliDhcTask.h
PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.cxx [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.h [new file with mode: 0644]
PWGCF/Correlations/macros/fd/AddTaskMCTruthCent.C [new file with mode: 0644]

index 40a48ce..9a4301b 100644 (file)
@@ -45,6 +45,7 @@ set ( SRCS
     Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
     Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx
     Correlations/DPhi/FourierDecomposition/AliMCTruthTrackMaker.cxx
+    Correlations/DPhi/FourierDecomposition/AliMCTruthCent.cxx
     Correlations/DPhi/AliLeadingV0Correlation.cxx 
     Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.cxx 
     Correlations/DPhi/MuonHadron/AliAnalysisTaskMuonHadronCorrelations.cxx 
index 3988cf6..9344c43 100644 (file)
@@ -38,7 +38,7 @@ AliDhcTask::AliDhcTask()
   fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
   fPtTACrit(kTRUE), fAllTAHists(kFALSE), fMixInEtaT(kFALSE),
   fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0), fOmitFirstEv(kTRUE),
-  fDoFillSame(kFALSE), fDoMassCut(kFALSE), fClassName(),
+  fDoFillSame(kFALSE), fDoMassCut(kFALSE), fCheckVertex(kTRUE), fClassName(),
   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36), fTriggerMatch(kTRUE),
   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
   fMixBCent(0x0), fMixBZvtx(0x0), fHEffT(0x0), fHEffA(0x0),
@@ -61,7 +61,7 @@ AliDhcTask::AliDhcTask(const char *name, Bool_t def)
   fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
   fPtTACrit(kTRUE), fAllTAHists(kFALSE), fMixInEtaT(kFALSE),
   fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0), fOmitFirstEv(kTRUE),
-  fDoFillSame(kFALSE), fDoMassCut(kFALSE), fClassName(),
+  fDoFillSame(kFALSE), fDoMassCut(kFALSE), fCheckVertex(kTRUE), fClassName(),
   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36), fTriggerMatch(kTRUE),
   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
   fMixBCent(0x0), fMixBZvtx(0x0), fHEffT(0x0), fHEffA(0x0),
@@ -143,6 +143,7 @@ void AliDhcTask::PrintDhcSettings()
   AliInfo(Form(" fill same event in any case? %d", fDoFillSame));
   AliInfo(Form(" do invariant mass cut? %d", fDoMassCut));
   AliInfo(Form(" omit first event? %d\n", fOmitFirstEv));
+  AliInfo(Form(" check the vertex? %d\n", fCheckVertex));
 }
 
 //________________________________________________________________________
@@ -434,7 +435,7 @@ void AliDhcTask::UserExec(Option_t *)
     if (dType == kESD) {
       const AliESDVertex* vertex = fESD->GetPrimaryVertex();
       fZVertex = vertex->GetZ();
-      if (!VertexOk()) {
+      if (fCheckVertex && !VertexOk()) {
         if (fVerbosity > 1)
           AliInfo(Form("Event REJECTED (ESD vertex not OK). z = %.1f", fZVertex));
         return;
@@ -446,7 +447,7 @@ void AliDhcTask::UserExec(Option_t *)
     } else if (dType == kAOD) {
       const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
       fZVertex = vertex->GetZ();
-      if (!VertexOk()) {
+      if (fCheckVertex && !VertexOk()) {
         if (fVerbosity > 1)
           Info("Exec()", "Event REJECTED (AOD vertex not OK). z = %.1f", fZVertex);
         return;
index cda9105..112a227 100644 (file)
@@ -40,6 +40,7 @@ class AliDhcTask : public AliAnalysisTaskSE {
   void         SetDEtaDPhiBins(Int_t nbe, Int_t nbp)  { fNBdeta=nbe; fNBdphi=nbp; }
   void         SetDoFillSame(Bool_t b)                { fDoFillSame = b;          }
   void         SetDoMassCut(Bool_t b)                 { fDoMassCut = b;           }
+  void         SetCheckVertex(Bool_t b)               { fCheckVertex = b;         }
   void         SetDoWeights(Bool_t b)                 { fDoWeights = b;           }
   void         SetEtaMax(Double_t eta)                { fEtaMax = eta;            }
   void         SetEtaTRange(Double_t eL, Double_t eH) { fEtaTLo=eL; fEtaTHi=eH;   }
@@ -107,6 +108,7 @@ class AliDhcTask : public AliAnalysisTaskSE {
   Bool_t             fOmitFirstEv;     //  if true skip first event in chunk
   Bool_t             fDoFillSame;      //  If true fill same event immediately (not waiting for pool)
   Bool_t             fDoMassCut;       //  If true cut on invariant mass
+  Bool_t             fCheckVertex;     //  switch to flase for MC generator-level analysis
   TString            fClassName;       //  If not null only process events with given class
   TString            fCentMethod;      //  centrality selection method
   Int_t              fNBdeta;          //  no. deta bins
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.cxx b/PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.cxx
new file mode 100644 (file)
index 0000000..4b2ba89
--- /dev/null
@@ -0,0 +1,221 @@
+//
+// Determine MC truth centrality for kinematics-only productions
+// pass it via centrality object
+//
+
+#include "AliMCTruthCent.h"
+
+#include <TString.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TChain.h>
+#include <TList.h>
+
+#include "AliAnalysisManager.h"
+#include "AliVEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliVParticle.h"
+#include "AliLog.h"
+#include "AliESDEvent.h"
+#include "AliCentrality.h"
+
+ClassImp(AliMCTruthCent)
+
+//________________________________________________________________________
+AliMCTruthCent::AliMCTruthCent() :
+AliAnalysisTaskSE("AliMCTruthCent"),
+fOutputList(0x0),
+fFillHistos(kFALSE),
+fHMultV0A(0x0), fHMultV0C(0x0), fHMultV0M(0x0), fHMultV0AvsV0C(0x0),
+fHCentV0A(0x0), fHCentV0C(0x0), fHCentV0M(0x0), fHCentV0AvsV0C(0x0),
+fV0ALo(0.0), fV0AHi(100.0),
+fV0CLo(0.0), fV0CHi(100.0),
+fV0MLo(0.0), fV0MHi(200.0)
+{
+  // Constructor.
+}
+
+//________________________________________________________________________
+AliMCTruthCent::AliMCTruthCent(const char *name) :
+AliAnalysisTaskSE(name),
+fOutputList(0x0),
+fFillHistos(kFALSE),
+fHMultV0A(0x0), fHMultV0C(0x0), fHMultV0M(0x0), fHMultV0AvsV0C(0x0),
+fHCentV0A(0x0), fHCentV0C(0x0), fHCentV0M(0x0), fHCentV0AvsV0C(0x0),
+fV0ALo(0.0), fV0AHi(100.0),
+fV0CLo(0.0), fV0CHi(100.0),
+fV0MLo(0.0), fV0MHi(200.0)
+{
+  // Constructor.
+  
+  // Define input slot
+  DefineInput(0, TChain::Class());
+
+}
+
+//________________________________________________________________________
+AliMCTruthCent::~AliMCTruthCent()
+{
+  // Destructor.
+  if (fOutputList)
+    delete fOutputList;
+}
+
+//________________________________________________________________________
+void AliMCTruthCent::UserCreateOutputObjects()
+{
+  // Create my user objects.
+  
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    AliFatal("No analysis manager!");
+    return;
+  }
+  
+  AliVEventHandler *handler = mgr->GetInputEventHandler();
+  if (!handler) {
+    AliFatal("No input handler!");
+    return;
+  }
+  
+  if (!handler->InheritsFrom("AliESDInputHandler")) {
+    AliError("Sorry, this only works for kinematics trees ... return");
+    return;
+  }
+  
+  if (fFillHistos) {
+    fOutputList = new TList();
+    fOutputList->SetOwner();
+    
+    // histograms
+    fHMultV0A      = new TH1D("fHMultV0A","V0A multiplicity",500,0.0,500.0);
+    fHMultV0C      = new TH1D("fHMultV0C","V0C multiplicity",500,0.0,500.0);
+    fHMultV0M      = new TH1D("fHMultV0M","V0M multiplicity",500,0.0,500.0);
+    fHMultV0AvsV0C = new TH2D("fHMultV0AvsV0C","V0A vs. C multiplicity",500,0.0,500.0,500,0.0,500.0);
+    
+    fOutputList->Add(fHMultV0A);
+    fOutputList->Add(fHMultV0C);
+    fOutputList->Add(fHMultV0M);
+    fOutputList->Add(fHMultV0AvsV0C);
+    
+    fHCentV0A      = new TH1D("fHCentV0A","V0A centrality",100,0.0,100.0);
+    fHCentV0C      = new TH1D("fHCentV0C","V0C centrality",100,0.0,100.0);
+    fHCentV0M      = new TH1D("fHCentV0M","V0M centrality",100,0.0,100.0);
+    fHCentV0AvsV0C = new TH2D("fHCentV0AvsV0C","V0A vs. C centrality",100,0.0,100.0,100,0.0,100.0);
+    
+    fOutputList->Add(fHCentV0A);
+    fOutputList->Add(fHCentV0C);
+    fOutputList->Add(fHCentV0M);
+    fOutputList->Add(fHCentV0AvsV0C);
+    
+    PostData(1, fOutputList);
+  }
+}
+
+//________________________________________________________________________
+void AliMCTruthCent::UserExec(Option_t *)
+{
+  // Main loop, called for each event.
+  if (!InputEvent()) {
+    AliError("Could not retrieve event! Returning");
+    return;
+  }
+  
+  if (!MCEvent()) {
+    AliError("Could not retrieve ESD MC event! Returning");
+    return;
+  }
+  
+  // can I get the centrality object from the InputEvent?
+  // then use it and overwrite the centrality
+  AliCentrality *esdCent = InputEvent()->GetCentrality();
+  esdCent->SetQuality(0);
+  
+  // define experiment-like estimators: multiplicity in range of
+  // V0A  2.8 < eta <  5.1
+  // V0C -3.7 < eta < -1.7
+  // V0M: sum of V0A and V0C
+  Double_t etaV0ALo = 2.8;
+  Double_t etaV0AHi = 5.1;
+  Double_t etaV0CLo = -3.7;
+  Double_t etaV0CHi = -1.7;
+  
+  Int_t nV0A = 0;
+  Int_t nV0C = 0;
+  
+  // loop over MC tracks
+  Double_t       etaTrack = 0.0;
+  AliVParticle * track    = 0x0;
+  const Int_t Ntracks = MCEvent()->GetNumberOfTracks();
+  for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
+    track = GetTrack(iTracks);
+    if (!track)
+      continue;
+    
+    if (track->Charge()==0)
+      continue;
+
+    etaTrack = track->Eta();
+    
+    if (etaTrack>etaV0ALo && etaTrack<etaV0AHi) {
+      ++nV0A;
+    } else if (etaTrack>etaV0CLo && etaTrack<etaV0CHi) {
+      ++nV0C;
+    }
+  }
+  Int_t nV0M = nV0A+nV0C;
+
+  
+  // linear centrality approximation
+  Double_t centV0A = (fV0AHi-nV0A) / (fV0AHi-fV0ALo);
+  centV0A *= 100.0;
+  if (centV0A >= 100.0)
+    centV0A = 99.9;
+  if (centV0A < 0.0)
+    centV0A = 0.0;
+  esdCent->SetCentralityV0A(centV0A);
+
+  Double_t centV0C = (fV0CHi-nV0C) / (fV0CHi-fV0CLo);
+  centV0C *= 100.0;
+  if (centV0C >= 100.0)
+    centV0C = 99.9;
+  if (centV0C < 0.0)
+    centV0C = 0.0;
+  esdCent->SetCentralityV0C(centV0C);
+  
+  Double_t centV0M = (fV0MHi-nV0M) / (fV0MHi-fV0MLo);
+  centV0M *= 100.0;
+  if (centV0M >= 100.0)
+    centV0M = 99.9;
+  if (centV0M < 0.0)
+    centV0M = 0.0;
+  esdCent->SetCentralityV0M(centV0M);
+
+  if (fFillHistos) {
+    // write out distributions to histograms
+    fHMultV0A->Fill(nV0A);
+    fHMultV0C->Fill(nV0C);
+    fHMultV0M->Fill(nV0M);
+    fHMultV0AvsV0C->Fill(nV0A,nV0C);
+    
+    fHCentV0A->Fill(centV0A);
+    fHCentV0C->Fill(centV0C);
+    fHCentV0M->Fill(centV0M);
+    fHCentV0AvsV0C->Fill(centV0A,centV0C);
+    
+    PostData(1, fOutputList);
+  }
+
+  // next step: get centrality from histogram...
+  
+}
+
+//________________________________________________________________________
+AliVParticle* AliMCTruthCent::GetTrack(Int_t i)
+{
+  if (!MCEvent()->IsPhysicalPrimary(i))
+    return 0;
+  
+  return MCEvent()->GetTrack(i);
+}
+
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.h b/PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.h
new file mode 100644 (file)
index 0000000..a7156a2
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIMCTRUTHCENT_H
+#define ALIMCTRUTHCENT_H
+
+class TList;
+class TH1;
+class TH2;
+
+class AliESDEvent;
+class AliMCEvent;
+class AliVParticle;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliMCTruthCent : public AliAnalysisTaskSE {
+ public:
+  AliMCTruthCent();
+  AliMCTruthCent(const char *name);
+  virtual ~AliMCTruthCent();
+
+  void SetV0ARange(Double_t mL, Double_t mH)     { fV0ALo = mL; fV0AHi = mH; }
+  void SetV0CRange(Double_t mL, Double_t mH)     { fV0CLo = mL; fV0CHi = mH; }
+  void SetV0MRange(Double_t mL, Double_t mH)     { fV0MLo = mL; fV0MHi = mH; }
+  void SetFillHistos()                           { fFillHistos=kTRUE; DefineOutput(1, TList::Class()); }
+  void UserCreateOutputObjects();
+  void UserExec(Option_t *option);
+
+ protected:
+  AliVParticle      *GetTrack(Int_t i);
+  
+  TList            * fOutputList;           //! Output list
+  Bool_t             fFillHistos;           //! flag to fill the QA histos
+  TH1D             * fHMultV0A;             //!
+  TH1D             * fHMultV0C;             //!
+  TH1D             * fHMultV0M;             //!
+  TH2D             * fHMultV0AvsV0C;        //!
+  TH1D             * fHCentV0A;             //!
+  TH1D             * fHCentV0C;             //!
+  TH1D             * fHCentV0M;             //!
+  TH2D             * fHCentV0AvsV0C;        //!
+  Double_t           fV0ALo;                //! for linear centrality approximation
+  Double_t           fV0AHi;                //!
+  Double_t           fV0CLo;                //!
+  Double_t           fV0CHi;                //!
+  Double_t           fV0MLo;                //!
+  Double_t           fV0MHi;                //!
+  
+ private:
+  AliMCTruthCent(const AliMCTruthCent&);            // not implemented
+  AliMCTruthCent &operator=(const AliMCTruthCent&); // not implemented
+
+  ClassDef(AliMCTruthCent, 1); // Task to select tracks in MC events
+};
+#endif
diff --git a/PWGCF/Correlations/macros/fd/AddTaskMCTruthCent.C b/PWGCF/Correlations/macros/fd/AddTaskMCTruthCent.C
new file mode 100644 (file)
index 0000000..a035223
--- /dev/null
@@ -0,0 +1,47 @@
+AliMCTruthCent* AddTaskMCTruthCent(Bool_t doHists = kFALSE)
+{
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskMCTruthCent", "No analysis manager to connect to.");
+    return NULL;
+  }
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskMCTruthCent", "This task requires an input event handler");
+    return NULL;
+  }
+  
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+  TString name("AliMCTruthCent");
+  AliMCTruthCent *eTask = new AliMCTruthCent(name);
+  if (doHists) {
+    eTask->SetFillHistos();
+  }
+  
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+  mgr->AddTask(eTask);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+  mgr->ConnectInput  (eTask, 0,  cinput1 );
+  
+  if (doHists) {
+    AliAnalysisDataContainer *coutput = mgr->CreateContainer("MCTruthCentStat",
+                                                             TList::Class(),
+                                                             AliAnalysisManager::kOutputContainer,
+                                                             "MCTruthCent.root");
+    mgr->ConnectOutput(eTask,1,coutput);
+  }
+  
+  return eTask;
+}