PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx

   1 //====================================================================
   2 #include "AliForwarddNdetaTask.h"
   3 #include <TMath.h>
   4 #include <TH2D.h>
   5 #include <TH1D.h>
   6 #include <THStack.h>
   7 #include <TList.h>
   8 #include <TFile.h>
   9 #include <AliAnalysisManager.h>
  10 #include <AliAODEvent.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODInputHandler.h>
  13 #include "AliForwardUtil.h"
  14 #include "AliAODForwardMult.h"
  15
  16 //____________________________________________________________________
  17 AliForwarddNdetaTask::AliForwarddNdetaTask()
  18   : AliBasedNdetaTask()
  19 {
  20   //
  21   // Constructor
  22   //
  23 }
  24
  25 //____________________________________________________________________
  26 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
  27   : AliBasedNdetaTask("Forward")
  28 {
  29   //
  30   // Constructor
  31   //
  32   // Paramters
  33   //   name    Name of task
  34 }
  35
  36 //____________________________________________________________________
  37 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
  38   : AliBasedNdetaTask(o)
  39 {
  40   //
  41   // Copy constructor
  42   //
  43 }
  44
  45 //____________________________________________________________________
  46 AliBasedNdetaTask::CentralityBin*
  47 AliForwarddNdetaTask::MakeCentralityBin(const char* name, Short_t l, Short_t h)
  48   const
  49 {
  50   //
  51   // Make a new centrality bin
  52   //
  53   // Parameters:
  54   //    name   Histogram names
  55   //    l      Lower cut
  56   //    h      Upper cut
  57   //
  58   // Return:
  59   //    Newly allocated object (of our type)
  60   //
  61   return new AliForwarddNdetaTask::CentralityBin(name, l, h);
  62 }
  63
  64 //____________________________________________________________________
  65 void
  66 AliForwarddNdetaTask::UserExec(Option_t* option)
  67 {
  68   //
  69   // Called at each event
  70   //
  71   // This is called once in the master
  72   //
  73   // Parameters:
  74   //    option Not used
  75   //
  76   AliBasedNdetaTask::UserExec(option);
  77
  78   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
  79   if (!aod) {
  80     AliError("Cannot get the AOD event");
  81     return;
  82   }
  83
  84   TObject* obj = aod->FindListObject("Forward");
  85   if (!obj) {
  86     AliWarning("No forward object found");
  87     return;
  88   }
  89   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
  90
  91   TObject* oPrimary   = aod->FindListObject("primary");
  92   if (!oPrimary) return;
  93
  94   TH2D* primary   = static_cast<TH2D*>(oPrimary);
  95
  96   // Loop over centrality bins
  97   TIter next(fListOfCentralities);
  98   CentralityBin* bin = 0;
  99   while ((bin = static_cast<CentralityBin*>(next())))
 100     bin->ProcessPrimary(forward, fTriggerMask, fVtxMin, fVtxMax, primary);
 101 }
 102
 103 //________________________________________________________________________
 104 void
 105 AliForwarddNdetaTask::Terminate(Option_t *option)
 106 {
 107   //
 108   // Called at end of event processing..
 109   //
 110   // This is called once in the master
 111   //
 112   // Parameters:
 113   //    option Not used
 114   AliBasedNdetaTask::Terminate(option);
 115
 116   THStack* truth      = new THStack("dndetaTruth",      "dN/d#eta MC Truth");
 117   THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin),
 118                                     "dN/d#eta MC Truth");
 119
 120   TIter next(fListOfCentralities);
 121   CentralityBin* bin = 0;
 122   while ((bin = static_cast<CentralityBin*>(next()))) {
 123     if (fCentAxis && bin->IsAllBin()) continue;
 124
 125     TList* results = bin->GetResults();
 126     if (!results) continue;
 127
 128     TH1* dndeta      = static_cast<TH1*>(results->FindObject("dndetaTruth"));
 129     TH1* dndetaRebin =
 130       static_cast<TH1*>(results->FindObject(Form("dndetaTruth_rebin%02d",
 131                                                 fRebin)));
 132     if (dndeta)      truth->Add(dndeta);
 133     if (dndetaRebin) truthRebin->Add(dndetaRebin);
 134   }
 135   // If available output rebinned stack
 136   if (!truth->GetHists() ||
 137       truth->GetHists()->GetEntries() <= 0) {
 138     AliWarning("No MC truth histograms found");
 139     delete truth;
 140     truth = 0;
 141   }
 142   if (truth) fOutput->Add(truth);
 143
 144   // If available output rebinned stack
 145   if (!truthRebin->GetHists() ||
 146       truthRebin->GetHists()->GetEntries() <= 0) {
 147     AliWarning("No rebinned MC truth histograms found");
 148     delete truthRebin;
 149     truthRebin = 0;
 150   }
 151   if (truthRebin) fOutput->Add(truthRebin);
 152
 153 }
 154
 155 //____________________________________________________________________
 156 TH2D*
 157 AliForwarddNdetaTask::GetHistogram(const AliAODEvent* aod, Bool_t mc)
 158 {
 159   //
 160   // Retrieve the histogram
 161   //
 162   // Parameters:
 163   //    aod AOD event
 164   //    mc  Whether to get the MC histogram or not
 165   //
 166   // Return:
 167   //    Retrieved histogram or null
 168   //
 169   TObject* obj = 0;
 170   if (mc) obj = aod->FindListObject("ForwardMC");
 171   else    obj = aod->FindListObject("Forward");
 172
 173   // We should have a forward object at least
 174   if (!obj) {
 175     if (!mc) AliWarning("No Forward object found AOD");
 176     return 0;
 177   }
 178   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 179   return &(forward->GetHistogram());
 180 }
 181
 182 //========================================================================
 183 void
 184 AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult*
 185                                                     forward,
 186                                                     Int_t triggerMask,
 187                                                     Double_t vzMin,
 188                                                     Double_t vzMax,
 189                                                     const TH2D* primary)
 190 {
 191   // Check the centrality class unless this is the 'all' bin
 192   if (!IsAllBin()) {
 193     Double_t centrality = forward->GetCentrality();
 194     if (centrality < fLow || centrality >= fHigh) return;
 195   }
 196
 197   // Create sum histogram
 198   if (!fSumPrimary) {
 199     fSumPrimary = static_cast<TH2D*>(primary->Clone("truth"));
 200     fSumPrimary->SetDirectory(0);
 201     fSumPrimary->Reset();
 202     fSums->Add(fSumPrimary);
 203   }
 204
 205   // translate real trigger mask to MC trigger mask
 206   Int_t mask = AliAODForwardMult::kB;
 207   if (triggerMask == AliAODForwardMult::kNSD) {
 208     mask ^= AliAODForwardMult::kNSD;
 209     mask =  AliAODForwardMult::kMCNSD;
 210   }
 211
 212   // Now use our normal check, but with the new mask, except
 213   vzMin = vzMax = -10000; // ignore vertex
 214   if (!forward->CheckEvent(mask, vzMin, vzMax, 0, 0, 0)) return;
 215
 216   fSumPrimary->Add(primary);
 217   Int_t n = Int_t(fSumPrimary->GetBinContent(0,0));
 218   fSumPrimary->SetBinContent(0,0, ++n);
 219 }
 220
 221 //________________________________________________________________________
 222 void
 223 AliForwarddNdetaTask::CentralityBin::End(TList*      sums,
 224                                          TList*      results,
 225                                          UShort_t    scheme,
 226                                          const TH1*  shapeCorr,
 227                                          Double_t    trigEff,
 228                                          Bool_t      symmetrice,
 229                                          Int_t       rebin,
 230                                          Bool_t      rootProj,
 231                                          Bool_t      corrEmpty,
 232                                          Bool_t      cutEdges,
 233                                          Int_t       triggerMask,
 234                                          Int_t       marker)
 235 {
 236   AliInfo(Form("In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d",
 237                GetName(), corrEmpty, cutEdges, rootProj));
 238   AliBasedNdetaTask::CentralityBin::End(sums, results, scheme,
 239                                         shapeCorr, trigEff,
 240                                         symmetrice, rebin,
 241                                         rootProj, corrEmpty, cutEdges,
 242                                         triggerMask, marker);
 243
 244   fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("truth"));
 245
 246   if (fSumPrimary) {
 247 #if 0
 248     Int_t n = fSumPrimary->GetBinContent(0,0);
 249 #else
 250     Int_t n = (triggerMask == AliAODForwardMult::kNSD ?
 251                Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) :
 252                Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
 253 #endif
 254     AliInfo(Form("Normalising MC truth to %d", n));
 255
 256     TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
 257                                                fSumPrimary->GetNbinsY(),"e");
 258     dndetaTruth->SetDirectory(0);
 259     dndetaTruth->Scale(1./n, "width");
 260
 261
 262     SetHistogramAttributes(dndetaTruth, GetColor(), 30, "Monte-Carlo truth");
 263
 264     fOutput->Add(dndetaTruth);
 265     fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
 266
 267     // Get analysis result, and form ratio
 268     TH1D* dndeta = static_cast<TH1D*>(fOutput->FindObject(Form("dndeta%s",
 269                                                                GetName())));
 270     if (!dndeta) {
 271       AliWarning(Form("No dndeta%s in the list %s",
 272                       GetName(), fOutput->GetName()));
 273     }
 274     else {
 275       TH1D* ratio = static_cast<TH1D*>(dndeta->Clone("ratio"));
 276       ratio->SetDirectory(0);
 277       ratio->Divide(dndetaTruth);
 278
 279       fOutput->Add(ratio);
 280       fOutput->Add(Rebin(ratio, rebin, cutEdges));
 281     }
 282   }
 283
 284   if (!IsAllBin()) return;
 285   TFile* file = TFile::Open("forward.root", "READ");
 286   if (!file) return;
 287
 288   TList* forward = static_cast<TList*>(file->Get("Forward"));
 289   if (!forward) {
 290     AliError("List Forward not found in forward.root");
 291     return;
 292   }
 293   TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
 294   if (!rings) {
 295     AliError("List ringResults not found in forward.root");
 296     return;
 297   }
 298   THStack* res = static_cast<THStack*>(rings->FindObject("all"));
 299   if (!res) {
 300     AliError(Form("Stack all not found in %s", rings->GetName()));
 301     return;
 302   }
 303   if (!fTriggers) {
 304     AliError("Triggers histogram not set");
 305     return;
 306   }
 307   Double_t ntotal   = 0;
 308   Double_t epsilonT = trigEff;
 309   // TEMPORARY FIX
 310   if (triggerMask == AliAODForwardMult::kNSD) {
 311     // This is a local change
 312     epsilonT = 0.92;
 313     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
 314   }
 315   AliInfo("Adding per-ring histograms to output");
 316   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
 317   TIter next(res->GetHists());
 318   TH1*  hist = 0;
 319   while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
 320   res->SetName("dndetaRings");
 321   fOutput->Add(res);
 322 }
 323
 324 //________________________________________________________________________
 325 //
 326 // EOF
 327 //