PWGCF/Correlations/DPhi/FourierDecomposition/AliMCTruthCent.cxx

   1 //
   2 // Determine MC truth centrality for kinematics-only productions
   3 // pass it via centrality object
   4 //
   5
   6 #include "AliMCTruthCent.h"
   7
   8 #include <TString.h>
   9 #include <TH1.h>
  10 #include <TH2.h>
  11 #include <TChain.h>
  12 #include <TList.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliVEventHandler.h"
  16 #include "AliMCEvent.h"
  17 #include "AliVParticle.h"
  18 #include "AliLog.h"
  19 #include "AliESDEvent.h"
  20 #include "AliCentrality.h"
  21
  22 ClassImp(AliMCTruthCent)
  23
  24 //________________________________________________________________________
  25 AliMCTruthCent::AliMCTruthCent() :
  26 AliAnalysisTaskSE("AliMCTruthCent"),
  27 fOutputList(0x0),
  28 fFillHistos(kFALSE),
  29 fHMultV0A(0x0), fHMultV0C(0x0), fHMultV0M(0x0), fHMultV0AvsV0C(0x0),
  30 fHCentV0A(0x0), fHCentV0C(0x0), fHCentV0M(0x0), fHCentV0AvsV0C(0x0),
  31 fV0ALo(0.0), fV0AHi(100.0),
  32 fV0CLo(0.0), fV0CHi(100.0),
  33 fV0MLo(0.0), fV0MHi(200.0)
  34 {
  35   // Constructor.
  36 }
  37
  38 //________________________________________________________________________
  39 AliMCTruthCent::AliMCTruthCent(const char *name) :
  40 AliAnalysisTaskSE(name),
  41 fOutputList(0x0),
  42 fFillHistos(kFALSE),
  43 fHMultV0A(0x0), fHMultV0C(0x0), fHMultV0M(0x0), fHMultV0AvsV0C(0x0),
  44 fHCentV0A(0x0), fHCentV0C(0x0), fHCentV0M(0x0), fHCentV0AvsV0C(0x0),
  45 fV0ALo(0.0), fV0AHi(100.0),
  46 fV0CLo(0.0), fV0CHi(100.0),
  47 fV0MLo(0.0), fV0MHi(200.0)
  48 {
  49   // Constructor.
  50
  51   // Define input slot
  52   DefineInput(0, TChain::Class());
  53
  54 }
  55
  56 //________________________________________________________________________
  57 AliMCTruthCent::~AliMCTruthCent()
  58 {
  59   // Destructor.
  60   if (fOutputList)
  61     delete fOutputList;
  62 }
  63
  64 //________________________________________________________________________
  65 void AliMCTruthCent::UserCreateOutputObjects()
  66 {
  67   // Create my user objects.
  68
  69   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  70   if (!mgr) {
  71     AliFatal("No analysis manager!");
  72     return;
  73   }
  74
  75   AliVEventHandler *handler = mgr->GetInputEventHandler();
  76   if (!handler) {
  77     AliFatal("No input handler!");
  78     return;
  79   }
  80
  81   if (!handler->InheritsFrom("AliESDInputHandler")) {
  82     AliError("Sorry, this only works for kinematics trees ... return");
  83     return;
  84   }
  85
  86   if (fFillHistos) {
  87     fOutputList = new TList();
  88     fOutputList->SetOwner();
  89
  90     // histograms
  91     fHMultV0A      = new TH1D("fHMultV0A","V0A multiplicity",500,0.0,500.0);
  92     fHMultV0C      = new TH1D("fHMultV0C","V0C multiplicity",500,0.0,500.0);
  93     fHMultV0M      = new TH1D("fHMultV0M","V0M multiplicity",500,0.0,500.0);
  94     fHMultV0AvsV0C = new TH2D("fHMultV0AvsV0C","V0A vs. C multiplicity",500,0.0,500.0,500,0.0,500.0);
  95
  96     fOutputList->Add(fHMultV0A);
  97     fOutputList->Add(fHMultV0C);
  98     fOutputList->Add(fHMultV0M);
  99     fOutputList->Add(fHMultV0AvsV0C);
 100
 101     fHCentV0A      = new TH1D("fHCentV0A","V0A centrality",100,0.0,100.0);
 102     fHCentV0C      = new TH1D("fHCentV0C","V0C centrality",100,0.0,100.0);
 103     fHCentV0M      = new TH1D("fHCentV0M","V0M centrality",100,0.0,100.0);
 104     fHCentV0AvsV0C = new TH2D("fHCentV0AvsV0C","V0A vs. C centrality",100,0.0,100.0,100,0.0,100.0);
 105
 106     fOutputList->Add(fHCentV0A);
 107     fOutputList->Add(fHCentV0C);
 108     fOutputList->Add(fHCentV0M);
 109     fOutputList->Add(fHCentV0AvsV0C);
 110
 111     PostData(1, fOutputList);
 112   }
 113 }
 114
 115 //________________________________________________________________________
 116 void AliMCTruthCent::UserExec(Option_t *)
 117 {
 118   // Main loop, called for each event.
 119   if (!InputEvent()) {
 120     AliError("Could not retrieve event! Returning");
 121     return;
 122   }
 123
 124   if (!MCEvent()) {
 125     AliError("Could not retrieve ESD MC event! Returning");
 126     return;
 127   }
 128
 129   // can I get the centrality object from the InputEvent?
 130   // then use it and overwrite the centrality
 131   AliCentrality *esdCent = InputEvent()->GetCentrality();
 132   esdCent->SetQuality(0);
 133
 134   // define experiment-like estimators: multiplicity in range of
 135   // V0A  2.8 < eta <  5.1
 136   // V0C -3.7 < eta < -1.7
 137   // V0M: sum of V0A and V0C
 138   Double_t etaV0ALo = 2.8;
 139   Double_t etaV0AHi = 5.1;
 140   Double_t etaV0CLo = -3.7;
 141   Double_t etaV0CHi = -1.7;
 142
 143   Int_t nV0A = 0;
 144   Int_t nV0C = 0;
 145
 146   // loop over MC tracks
 147   Double_t       etaTrack = 0.0;
 148   AliVParticle * track    = 0x0;
 149   const Int_t Ntracks = MCEvent()->GetNumberOfTracks();
 150   for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 151     track = GetTrack(iTracks);
 152     if (!track)
 153       continue;
 154
 155     if (track->Charge()==0)
 156       continue;
 157
 158     etaTrack = track->Eta();
 159
 160     if (etaTrack>etaV0ALo && etaTrack<etaV0AHi) {
 161       ++nV0A;
 162     } else if (etaTrack>etaV0CLo && etaTrack<etaV0CHi) {
 163       ++nV0C;
 164     }
 165   }
 166   Int_t nV0M = nV0A+nV0C;
 167
 168
 169   // linear centrality approximation
 170   Double_t centV0A = (fV0AHi-nV0A) / (fV0AHi-fV0ALo);
 171   centV0A *= 100.0;
 172   if (centV0A >= 100.0)
 173     centV0A = 99.9;
 174   if (centV0A < 0.0)
 175     centV0A = 0.0;
 176   esdCent->SetCentralityV0A(centV0A);
 177
 178   Double_t centV0C = (fV0CHi-nV0C) / (fV0CHi-fV0CLo);
 179   centV0C *= 100.0;
 180   if (centV0C >= 100.0)
 181     centV0C = 99.9;
 182   if (centV0C < 0.0)
 183     centV0C = 0.0;
 184   esdCent->SetCentralityV0C(centV0C);
 185
 186   Double_t centV0M = (fV0MHi-nV0M) / (fV0MHi-fV0MLo);
 187   centV0M *= 100.0;
 188   if (centV0M >= 100.0)
 189     centV0M = 99.9;
 190   if (centV0M < 0.0)
 191     centV0M = 0.0;
 192   esdCent->SetCentralityV0M(centV0M);
 193
 194   if (fFillHistos) {
 195     // write out distributions to histograms
 196     fHMultV0A->Fill(nV0A);
 197     fHMultV0C->Fill(nV0C);
 198     fHMultV0M->Fill(nV0M);
 199     fHMultV0AvsV0C->Fill(nV0A,nV0C);
 200
 201     fHCentV0A->Fill(centV0A);
 202     fHCentV0C->Fill(centV0C);
 203     fHCentV0M->Fill(centV0M);
 204     fHCentV0AvsV0C->Fill(centV0A,centV0C);
 205
 206     PostData(1, fOutputList);
 207   }
 208
 209   // next step: get centrality from histogram...
 210
 211 }
 212
 213 //________________________________________________________________________
 214 AliVParticle* AliMCTruthCent::GetTrack(Int_t i)
 215 {
 216   if (!MCEvent()->IsPhysicalPrimary(i))
 217     return 0;
 218
 219   return MCEvent()->GetTrack(i);
 220 }
 221