PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx

   1 //
   2 // Calculate the multiplicity in the forward regions event-by-event
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODForwardMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliForwardMultiplicityTask.h"
  15 #include "AliTriggerAnalysis.h"
  16 #include "AliPhysicsSelection.h"
  17 #include "AliLog.h"
  18 #include "AliESDEvent.h"
  19 #include "AliAODHandler.h"
  20 #include "AliMultiplicity.h"
  21 #include "AliInputEventHandler.h"
  22 #include "AliForwardCorrectionManager.h"
  23 #include "AliAnalysisManager.h"
  24 #include <TH1.h>
  25 #include <TDirectory.h>
  26 #include <TTree.h>
  27 #include <TROOT.h>
  28
  29
  30 //====================================================================
  31 AliForwardMultiplicityTask::AliForwardMultiplicityTask()
  32   : AliForwardMultiplicityBase(),
  33     fHData(0),
  34     fESDFMD(),
  35     fHistos(),
  36     fAODFMD(),
  37     fAODEP(),
  38     fRingSums(),
  39     fEventInspector(),
  40     fSharingFilter(),
  41     fDensityCalculator(),
  42     fCorrections(),
  43     fHistCollector(),
  44     fEventPlaneFinder(),
  45     fList(0)
  46 {
  47   //
  48   // Constructor
  49   //
  50   DGUARD(0,0,"Default construction of AliForwardMultiplicityTask");
  51 }
  52
  53 //____________________________________________________________________
  54 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
  55   : AliForwardMultiplicityBase(name),
  56     fHData(0),
  57     fESDFMD(),
  58     fHistos(),
  59     fAODFMD(false),
  60     fAODEP(false),
  61     fRingSums(),
  62     fEventInspector("event"),
  63     fSharingFilter("sharing"),
  64     fDensityCalculator("density"),
  65     fCorrections("corrections"),
  66     fHistCollector("collector"),
  67     fEventPlaneFinder("eventplane"),
  68     fList(0)
  69 {
  70   //
  71   // Constructor
  72   //
  73   // Parameters:
  74   //    name Name of task
  75   //
  76   DGUARD(0,0,"named construction of AliForwardMultiplicityTask: %s", name);
  77   DefineOutput(1, TList::Class());
  78 }
  79
  80 //____________________________________________________________________
  81 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
  82   : AliForwardMultiplicityBase(o),
  83     fHData(o.fHData),
  84     fESDFMD(o.fESDFMD),
  85     fHistos(o.fHistos),
  86     fAODFMD(o.fAODFMD),
  87     fAODEP(o.fAODEP),
  88     fRingSums(o.fRingSums),
  89     fEventInspector(o.fEventInspector),
  90     fSharingFilter(o.fSharingFilter),
  91     fDensityCalculator(o.fDensityCalculator),
  92     fCorrections(o.fCorrections),
  93     fHistCollector(o.fHistCollector),
  94     fEventPlaneFinder(o.fEventPlaneFinder),
  95     fList(o.fList)
  96 {
  97   //
  98   // Copy constructor
  99   //
 100   // Parameters:
 101   //    o Object to copy from
 102   //
 103   DGUARD(0,0,"Copy construction of AliForwardMultiplicityTask");
 104   DefineOutput(1, TList::Class());
 105 }
 106
 107 //____________________________________________________________________
 108 AliForwardMultiplicityTask&
 109 AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
 110 {
 111   //
 112   // Assignment operator
 113   //
 114   // Parameters:
 115   //    o Object to assign from
 116   //
 117   // Return:
 118   //    Reference to this object
 119   //
 120   DGUARD(fDebug,3,"Assignment to AliForwardMultiplicityTask");
 121   if (&o == this) return *this;
 122   AliForwardMultiplicityBase::operator=(o);
 123
 124   fHData             = o.fHData;
 125   fEventInspector    = o.fEventInspector;
 126   fSharingFilter     = o.fSharingFilter;
 127   fDensityCalculator = o.fDensityCalculator;
 128   fCorrections       = o.fCorrections;
 129   fHistCollector     = o.fHistCollector;
 130   fEventPlaneFinder  = o.fEventPlaneFinder;
 131   fHistos            = o.fHistos;
 132   fAODFMD            = o.fAODFMD;
 133   fAODEP             = o.fAODEP;
 134   fRingSums          = o.fRingSums;
 135   fList              = o.fList;
 136
 137   return *this;
 138 }
 139
 140 //____________________________________________________________________
 141 void
 142 AliForwardMultiplicityTask::SetDebug(Int_t dbg)
 143 {
 144   //
 145   // Set debug level
 146   //
 147   // Parameters:
 148   //    dbg Debug level
 149   //
 150   fEventInspector.SetDebug(dbg);
 151   fSharingFilter.SetDebug(dbg);
 152   fDensityCalculator.SetDebug(dbg);
 153   fCorrections.SetDebug(dbg);
 154   fHistCollector.SetDebug(dbg);
 155   fEventPlaneFinder.SetDebug(dbg);
 156 }
 157
 158 //____________________________________________________________________
 159 Bool_t
 160 AliForwardMultiplicityTask::InitializeSubs()
 161 {
 162   //
 163   // Initialise the sub objects and stuff.  Called on first event
 164   //
 165   //
 166   DGUARD(fDebug,1,"Initialize sub-algorithms");
 167   const TAxis* pe = 0;
 168   const TAxis* pv = 0;
 169
 170   if (!ReadCorrections(pe,pv)) return false;
 171
 172   fHistos.Init(*pe);
 173   fAODFMD.Init(*pe);
 174   fAODEP.Init(*pe);
 175   fRingSums.Init(*pe);
 176
 177   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 178   fHData->SetStats(0);
 179   fHData->SetDirectory(0);
 180   fList->Add(fHData);
 181
 182   TList* rings = new TList;
 183   rings->SetName("ringSums");
 184   rings->SetOwner();
 185   fList->Add(rings);
 186
 187   rings->Add(fRingSums.Get(1, 'I'));
 188   rings->Add(fRingSums.Get(2, 'I'));
 189   rings->Add(fRingSums.Get(2, 'O'));
 190   rings->Add(fRingSums.Get(3, 'I'));
 191   rings->Add(fRingSums.Get(3, 'O'));
 192   fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 193   fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 194   fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 195   fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 196   fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 197
 198   fEventInspector.Init(*pv);
 199   fSharingFilter.Init(*pe);
 200   fDensityCalculator.Init(*pe);
 201   fCorrections.Init(*pe);
 202   fHistCollector.Init(*pv,*pe);
 203   fEventPlaneFinder.Init(*pe);
 204
 205   this->Print();
 206   return true;
 207 }
 208
 209 //____________________________________________________________________
 210 void
 211 AliForwardMultiplicityTask::UserCreateOutputObjects()
 212 {
 213   //
 214   // Create output objects
 215   //
 216   //
 217   DGUARD(fDebug,1,"Create user ouput");
 218   fList = new TList;
 219   fList->SetOwner();
 220
 221   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 222   AliAODHandler*      ah =
 223     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 224   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 225
 226
 227   TObject* obj = &fAODFMD;
 228   ah->AddBranch("AliAODForwardMult", &obj);
 229   TObject* epobj = &fAODEP;
 230   ah->AddBranch("AliAODForwardEP", &epobj);
 231
 232   fEventInspector.DefineOutput(fList);
 233   fSharingFilter.DefineOutput(fList);
 234   fDensityCalculator.DefineOutput(fList);
 235   fCorrections.DefineOutput(fList);
 236   fHistCollector.DefineOutput(fList);
 237   fEventPlaneFinder.DefineOutput(fList);
 238
 239   PostData(1, fList);
 240 }
 241 //____________________________________________________________________
 242 void
 243 AliForwardMultiplicityTask::UserExec(Option_t*)
 244 {
 245   //
 246   // Process each event
 247   //
 248   // Parameters:
 249   //    option Not used
 250   //
 251
 252   DGUARD(fDebug,1,"Process the input event");
 253   // static Int_t cnt = 0;
 254   // cnt++;
 255   // Get the input data
 256   AliESDEvent* esd = GetESDEvent();
 257   if (!esd) return;
 258
 259   // Clear stuff
 260   fHistos.Clear();
 261   fESDFMD.Clear();
 262   fAODFMD.Clear();
 263   fAODEP.Clear();
 264
 265   Bool_t   lowFlux   = kFALSE;
 266   UInt_t   triggers  = 0;
 267   UShort_t ivz       = 0;
 268   Double_t vz        = 0;
 269   Double_t cent      = -1;
 270   UShort_t nClusters = 0;
 271   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 272                                                ivz, vz, cent, nClusters);
 273
 274   if (found & AliFMDEventInspector::kNoEvent)    return;
 275   if (found & AliFMDEventInspector::kNoTriggers) return;
 276
 277   // Set trigger bits, and mark this event for storage
 278   fAODFMD.SetTriggerBits(triggers);
 279   fAODFMD.SetSNN(fEventInspector.GetEnergy());
 280   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 281   fAODFMD.SetCentrality(cent);
 282   fAODFMD.SetNClusters(nClusters);
 283   MarkEventForStore();
 284
 285   if (found & AliFMDEventInspector::kNoSPD)      return;
 286   if (found & AliFMDEventInspector::kNoFMD)      return;
 287   if (found & AliFMDEventInspector::kNoVertex)   return;
 288
 289   if (triggers & AliAODForwardMult::kPileUp) return;
 290
 291   fAODFMD.SetIpZ(vz);
 292
 293   if (found & AliFMDEventInspector::kBadVertex) return;
 294
 295   // We we do not want to use low flux specific code, we disable it here.
 296   if (!fEnableLowFlux) lowFlux = false;
 297
 298   // Get FMD data
 299   AliESDFMD* esdFMD = esd->GetFMDData();
 300   //  // Apply the sharing filter (or hit merging or clustering if you like)
 301   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 302     AliWarning("Sharing filter failed!");
 303     return;
 304   }
 305
 306   // Calculate the inclusive charged particle density
 307   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent,vz)) {
 308     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
 309     AliWarning("Density calculator failed!");
 310     return;
 311   }
 312
 313   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
 314     if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos))
 315       AliWarning("Eventplane finder failed!");
 316   }
 317
 318   // Do the secondary and other corrections.
 319   if (!fCorrections.Correct(fHistos, ivz)) {
 320     AliWarning("Corrections failed");
 321     return;
 322   }
 323
 324   if (!fHistCollector.Collect(fHistos, fRingSums,
 325                               ivz, fAODFMD.GetHistogram())) {
 326     AliWarning("Histogram collector failed");
 327     return;
 328   }
 329
 330   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 331     fHData->Add(&(fAODFMD.GetHistogram()));
 332
 333   PostData(1, fList);
 334 }
 335
 336 //____________________________________________________________________
 337 void
 338 AliForwardMultiplicityTask::Terminate(Option_t*)
 339 {
 340   //
 341   // End of job
 342   //
 343   // Parameters:
 344   //    option Not used
 345   //
 346   DGUARD(fDebug,1,"Processing the merged results");
 347
 348   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 349   if (!list) {
 350     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 351     if (GetOutputData(1)) GetOutputData(1)->Print();
 352     return;
 353   }
 354
 355   // Get our histograms from the container
 356   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 357   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 358   TH1I* hTriggers    = 0;
 359   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 360                                        hEventsTrVtx, hTriggers)) {
 361     AliError(Form("Didn't get histograms from event selector "
 362                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 363                   hEventsTr, hEventsTrVtx));
 364     list->ls();
 365     return;
 366   }
 367
 368   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 369   if (!hData) {
 370     AliError(Form("Couldn't get our summed histogram from output "
 371                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 372     list->ls();
 373     return;
 374   }
 375
 376   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 377   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 378   TH1D* norm   = hData->ProjectionX("norm",   0,  0,  "");
 379   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 380   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 381   dNdeta->Divide(norm);
 382   dNdeta->SetStats(0);
 383   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 384                 "width");
 385   list->Add(dNdeta);
 386   list->Add(norm);
 387
 388   MakeRingdNdeta(list, "ringSums", list, "ringResults");
 389
 390   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
 391   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 392   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 393 }
 394
 395 //
 396 // EOF
 397 //