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 void
 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;
 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();
 200   fDensityCalculator.Init(*pe);
 201   fCorrections.Init(*pe);
 202   fHistCollector.Init(*pv,*pe);
 203   fEventPlaneFinder.Init(*pe);
 204
 205   this->Print();
 206 }
 207
 208 //____________________________________________________________________
 209 void
 210 AliForwardMultiplicityTask::UserCreateOutputObjects()
 211 {
 212   //
 213   // Create output objects
 214   //
 215   //
 216   DGUARD(fDebug,1,"Create user ouput");
 217   fList = new TList;
 218   fList->SetOwner();
 219
 220   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 221   AliAODHandler*      ah =
 222     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 223   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 224
 225
 226   TObject* obj = &fAODFMD;
 227   ah->AddBranch("AliAODForwardMult", &obj);
 228   TObject* epobj = &fAODEP;
 229   ah->AddBranch("AliAODForwardEP", &epobj);
 230
 231   fEventInspector.DefineOutput(fList);
 232   fSharingFilter.DefineOutput(fList);
 233   fDensityCalculator.DefineOutput(fList);
 234   fCorrections.DefineOutput(fList);
 235   fHistCollector.DefineOutput(fList);
 236   fEventPlaneFinder.DefineOutput(fList);
 237
 238   PostData(1, fList);
 239 }
 240 //____________________________________________________________________
 241 void
 242 AliForwardMultiplicityTask::UserExec(Option_t*)
 243 {
 244   //
 245   // Process each event
 246   //
 247   // Parameters:
 248   //    option Not used
 249   //
 250
 251   DGUARD(fDebug,1,"Process the input event");
 252   // static Int_t cnt = 0;
 253   // cnt++;
 254   // Get the input data
 255   AliESDEvent* esd = GetESDEvent();
 256
 257   // Clear stuff
 258   fHistos.Clear();
 259   fESDFMD.Clear();
 260   fAODFMD.Clear();
 261   fAODEP.Clear();
 262
 263   Bool_t   lowFlux   = kFALSE;
 264   UInt_t   triggers  = 0;
 265   UShort_t ivz       = 0;
 266   Double_t vz        = 0;
 267   Double_t cent      = -1;
 268   UShort_t nClusters = 0;
 269   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 270                                                ivz, vz, cent, nClusters);
 271
 272   if (found & AliFMDEventInspector::kNoEvent)    return;
 273   if (found & AliFMDEventInspector::kNoTriggers) return;
 274
 275   // Set trigger bits, and mark this event for storage
 276   fAODFMD.SetTriggerBits(triggers);
 277   fAODFMD.SetSNN(fEventInspector.GetEnergy());
 278   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 279   fAODFMD.SetCentrality(cent);
 280   fAODFMD.SetNClusters(nClusters);
 281   MarkEventForStore();
 282
 283   if (found & AliFMDEventInspector::kNoSPD)      return;
 284   if (found & AliFMDEventInspector::kNoFMD)      return;
 285   if (found & AliFMDEventInspector::kNoVertex)   return;
 286
 287   if (triggers & AliAODForwardMult::kPileUp) return;
 288
 289   fAODFMD.SetIpZ(vz);
 290
 291   if (found & AliFMDEventInspector::kBadVertex) return;
 292
 293   // We we do not want to use low flux specific code, we disable it here.
 294   if (!fEnableLowFlux) lowFlux = false;
 295
 296   // Get FMD data
 297   AliESDFMD* esdFMD = esd->GetFMDData();
 298   //  // Apply the sharing filter (or hit merging or clustering if you like)
 299   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 300     AliWarning("Sharing filter failed!");
 301     return;
 302   }
 303
 304   // Calculate the inclusive charged particle density
 305   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent,vz)) {
 306     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
 307     AliWarning("Density calculator failed!");
 308     return;
 309   }
 310
 311   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
 312     if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos))
 313       AliWarning("Eventplane finder failed!");
 314   }
 315
 316   // Do the secondary and other corrections.
 317   if (!fCorrections.Correct(fHistos, ivz)) {
 318     AliWarning("Corrections failed");
 319     return;
 320   }
 321
 322   if (!fHistCollector.Collect(fHistos, fRingSums,
 323                               ivz, fAODFMD.GetHistogram())) {
 324     AliWarning("Histogram collector failed");
 325     return;
 326   }
 327
 328   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 329     fHData->Add(&(fAODFMD.GetHistogram()));
 330
 331   PostData(1, fList);
 332 }
 333
 334 //____________________________________________________________________
 335 void
 336 AliForwardMultiplicityTask::Terminate(Option_t*)
 337 {
 338   //
 339   // End of job
 340   //
 341   // Parameters:
 342   //    option Not used
 343   //
 344   DGUARD(fDebug,1,"Processing the merged results");
 345
 346   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 347   if (!list) {
 348     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 349     if (GetOutputData(1)) GetOutputData(1)->Print();
 350     return;
 351   }
 352
 353   // Get our histograms from the container
 354   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 355   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 356   TH1I* hTriggers    = 0;
 357   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 358                                        hEventsTrVtx, hTriggers)) {
 359     AliError(Form("Didn't get histograms from event selector "
 360                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 361                   hEventsTr, hEventsTrVtx));
 362     list->ls();
 363     return;
 364   }
 365
 366   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 367   if (!hData) {
 368     AliError(Form("Couldn't get our summed histogram from output "
 369                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 370     list->ls();
 371     return;
 372   }
 373
 374   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 375   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 376   TH1D* norm   = hData->ProjectionX("norm",   0,  0,  "");
 377   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 378   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 379   dNdeta->Divide(norm);
 380   dNdeta->SetStats(0);
 381   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 382                 "width");
 383   list->Add(dNdeta);
 384   list->Add(norm);
 385
 386   MakeRingdNdeta(list, "ringSums", list, "ringResults");
 387
 388   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
 389   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 390   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 391 }
 392
 393 //
 394 // EOF
 395 //