PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx

   1 //
   2 // Calculate the multiplicity in the forward regions event-by-event
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //   - Kinematics
   7 //   - Track references
   8 //
   9 // Outputs:
  10 //   - AliAODForwardMult
  11 //
  12 // Histograms
  13 //
  14 // Corrections used
  15 //
  16 #include "AliForwardMCMultiplicityTask.h"
  17 #include "AliTriggerAnalysis.h"
  18 #include "AliPhysicsSelection.h"
  19 #include "AliLog.h"
  20 #include "AliESDEvent.h"
  21 #include "AliAODHandler.h"
  22 #include "AliMultiplicity.h"
  23 #include "AliInputEventHandler.h"
  24 #include "AliForwardCorrectionManager.h"
  25 #include "AliAnalysisManager.h"
  26 #include <TH1.h>
  27 #include <TDirectory.h>
  28 #include <TTree.h>
  29 #include <TROOT.h>
  30 #include <iostream>
  31 #include <iomanip>
  32
  33 //====================================================================
  34 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
  35   : AliForwardMultiplicityBase(),
  36     fHData(0),
  37     fESDFMD(),
  38     fHistos(),
  39     fAODFMD(),
  40     fMCESDFMD(),
  41     fMCHistos(),
  42     fMCAODFMD(),
  43     fPrimary(0),
  44     fEventInspector(),
  45     fEnergyFitter(),
  46     fSharingFilter(),
  47     fDensityCalculator(),
  48     fCorrections(),
  49     fHistCollector(),
  50     fList(0)
  51 {
  52   //
  53   // Constructor
  54   //
  55 }
  56
  57 //____________________________________________________________________
  58 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  59   : AliForwardMultiplicityBase(name),
  60     fHData(0),
  61     fESDFMD(),
  62     fHistos(),
  63     fAODFMD(kFALSE),
  64     fMCESDFMD(),
  65     fMCHistos(),
  66     fMCAODFMD(kTRUE),
  67     fPrimary(0),
  68     fEventInspector("event"),
  69     fEnergyFitter("energy"),
  70     fSharingFilter("sharing"),
  71     fDensityCalculator("density"),
  72     fCorrections("corrections"),
  73     fHistCollector("collector"),
  74     fList(0)
  75 {
  76   //
  77   // Constructor
  78   //
  79   // Parameters:
  80   //    name Name of task
  81   //
  82   DefineOutput(1, TList::Class());
  83 }
  84
  85 //____________________________________________________________________
  86 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
  87   : AliForwardMultiplicityBase(o),
  88     fHData(o.fHData),
  89     fESDFMD(o.fESDFMD),
  90     fHistos(o.fHistos),
  91     fAODFMD(o.fAODFMD),
  92     fMCESDFMD(o.fMCESDFMD),
  93     fMCHistos(o.fMCHistos),
  94     fMCAODFMD(o.fMCAODFMD),
  95     fPrimary(o.fPrimary),
  96     fEventInspector(o.fEventInspector),
  97     fEnergyFitter(o.fEnergyFitter),
  98     fSharingFilter(o.fSharingFilter),
  99     fDensityCalculator(o.fDensityCalculator),
 100     fCorrections(o.fCorrections),
 101     fHistCollector(o.fHistCollector),
 102     fList(o.fList)
 103 {
 104   //
 105   // Copy constructor
 106   //
 107   // Parameters:
 108   //    o Object to copy from
 109   //
 110   DefineOutput(1, TList::Class());
 111 }
 112
 113 //____________________________________________________________________
 114 AliForwardMCMultiplicityTask&
 115 AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
 116 {
 117   //
 118   // Assignment operator
 119   //
 120   // Parameters:
 121   //    o Object to assign from
 122   //
 123   // Return:
 124   //    Reference to this object
 125   //
 126   AliForwardMultiplicityBase::operator=(o);
 127
 128   fHData             = o.fHData;
 129   fEventInspector    = o.fEventInspector;
 130   fEnergyFitter      = o.fEnergyFitter;
 131   fSharingFilter     = o.fSharingFilter;
 132   fDensityCalculator = o.fDensityCalculator;
 133   fCorrections       = o.fCorrections;
 134   fHistCollector     = o.fHistCollector;
 135   fHistos            = o.fHistos;
 136   fAODFMD            = o.fAODFMD;
 137   fMCHistos          = o.fMCHistos;
 138   fMCAODFMD          = o.fMCAODFMD;
 139   fPrimary           = o.fPrimary;
 140   fList              = o.fList;
 141
 142   return *this;
 143 }
 144
 145 //____________________________________________________________________
 146 void
 147 AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
 148 {
 149   //
 150   // Set debug level
 151   //
 152   // Parameters:
 153   //    dbg debug level
 154   //
 155   fEventInspector.SetDebug(dbg);
 156   fEnergyFitter.SetDebug(dbg);
 157   fSharingFilter.SetDebug(dbg);
 158   fDensityCalculator.SetDebug(dbg);
 159   fCorrections.SetDebug(dbg);
 160   fHistCollector.SetDebug(dbg);
 161 }
 162 //____________________________________________________________________
 163 void
 164 AliForwardMCMultiplicityTask::InitializeSubs()
 165 {
 166   //
 167   // Initialise the sub objects and stuff.  Called on first event
 168   //
 169   //
 170   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 171   fcm.Init(fEventInspector.GetCollisionSystem(),
 172            fEventInspector.GetEnergy(),
 173            fEventInspector.GetField(),
 174            true); // Last true is for MC flag
 175   // Check that we have the energy loss fits, needed by
 176   //   AliFMDSharingFilter
 177   //   AliFMDDensityCalculator
 178   if (!fcm.GetELossFit()) {
 179     AliFatal(Form("No energy loss fits"));
 180     return;
 181   }
 182   // Check that we have the double hit correction - (optionally) used by
 183   //  AliFMDDensityCalculator
 184   if (!fcm.GetDoubleHit()) {
 185     AliWarning("No double hit corrections");
 186   }
 187   // Check that we have the secondary maps, needed by
 188   //   AliFMDCorrections
 189   //   AliFMDHistCollector
 190   if (!fcm.GetSecondaryMap()) {
 191     AliFatal("No secondary corrections");
 192     return;
 193   }
 194   // Check that we have the vertex bias correction, needed by
 195   //   AliFMDCorrections
 196   if (!fcm.GetVertexBias()) {
 197     AliFatal("No event vertex bias corrections");
 198     return;
 199   }
 200   // Check that we have the merging efficiencies, optionally used by
 201   //   AliFMDCorrections
 202   if (!fcm.GetMergingEfficiency()) {
 203     AliWarning("No merging efficiencies");
 204   }
 205
 206
 207   const TAxis* pe = fcm.GetEtaAxis();
 208   const TAxis* pv = fcm.GetVertexAxis();
 209   if (!pe) AliFatal("No eta axis defined");
 210   if (!pv) AliFatal("No vertex axis defined");
 211
 212   fHistos.Init(*pe);
 213   fAODFMD.Init(*pe);
 214   fMCHistos.Init(*pe);
 215   fMCAODFMD.Init(*pe);
 216
 217   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 218   fHData->SetStats(0);
 219   fHData->SetDirectory(0);
 220
 221   fList->Add(fHData);
 222
 223   fEnergyFitter.Init(*pe);
 224   fEventInspector.Init(*pv);
 225   fDensityCalculator.Init(*pe);
 226   fCorrections.Init(*pe);
 227   fHistCollector.Init(*pv);
 228
 229   this->Print();
 230 }
 231
 232 //____________________________________________________________________
 233 void
 234 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
 235 {
 236   //
 237   // Create output objects
 238   //
 239   //
 240   fList = new TList;
 241
 242   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 243   AliAODHandler*      ah =
 244     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 245   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 246
 247
 248   TObject* obj = &fAODFMD;
 249   ah->AddBranch("AliAODForwardMult", &obj);
 250
 251   TObject* mcobj = &fMCAODFMD;
 252   ah->AddBranch("AliAODForwardMult", &mcobj);
 253
 254   fPrimary = new TH2D("primary", "MC Primaries",
 255                       200, -4, 6, 20, 0, 2*TMath::Pi());
 256   fPrimary->SetXTitle("#eta");
 257   fPrimary->SetYTitle("#varphi [radians]");
 258   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 259   fPrimary->Sumw2();
 260   fPrimary->SetStats(0);
 261   fPrimary->SetDirectory(0);
 262   ah->AddBranch("TH2D", &fPrimary);
 263
 264   fEventInspector.DefineOutput(fList);
 265   fEnergyFitter.DefineOutput(fList);
 266   fSharingFilter.DefineOutput(fList);
 267   fDensityCalculator.DefineOutput(fList);
 268   fCorrections.DefineOutput(fList);
 269
 270   PostData(1, fList);
 271 }
 272 //____________________________________________________________________
 273 void
 274 AliForwardMCMultiplicityTask::UserExec(Option_t*)
 275 {
 276   //
 277   // Process each event
 278   //
 279   // Parameters:
 280   //    option Not used
 281   //
 282
 283   // Get the input data
 284   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 285   if (!esd) {
 286     AliWarning("No ESD event found for input event");
 287     return;
 288   }
 289
 290   // On the first event, initialize the parameters
 291   if (fFirstEvent && esd->GetESDRun()) {
 292     fEventInspector.ReadRunDetails(esd);
 293
 294     AliInfo(Form("Initializing with parameters from the ESD:\n"
 295                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 296                  "         AliESDEvent::GetBeamType()     ->%s\n"
 297                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 298                  "         AliESDEvent::GetMagneticField()->%f\n"
 299                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 300                  esd->GetBeamEnergy(),
 301                  esd->GetBeamType(),
 302                  esd->GetCurrentL3(),
 303                  esd->GetMagneticField(),
 304                  esd->GetRunNumber()));
 305
 306     fFirstEvent = false;
 307
 308     InitializeSubs();
 309   }
 310   // Clear stuff
 311   fHistos.Clear();
 312   fESDFMD.Clear();
 313   fAODFMD.Clear();
 314   fMCHistos.Clear();
 315   fMCESDFMD.Clear();
 316   fMCAODFMD.Clear();
 317   fPrimary->Reset();
 318
 319   Bool_t   lowFlux  = kFALSE;
 320   UInt_t   triggers = 0;
 321   UShort_t ivz      = 0;
 322   Double_t vz       = 0;
 323   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
 324   if (found & AliFMDEventInspector::kNoEvent)    return;
 325   if (found & AliFMDEventInspector::kNoTriggers) return;
 326
 327   // Set trigger bits, and mark this event for storage
 328   fAODFMD.SetTriggerBits(triggers);
 329   fMCAODFMD.SetTriggerBits(triggers);
 330   MarkEventForStore();
 331
 332   if (found & AliFMDEventInspector::kNoSPD)     return;
 333   if (found & AliFMDEventInspector::kNoFMD)     return;
 334   if (found & AliFMDEventInspector::kNoVertex)  return;
 335   fAODFMD.SetIpZ(vz);
 336   fMCAODFMD.SetIpZ(vz);
 337
 338   if (found & AliFMDEventInspector::kBadVertex) return;
 339
 340   // We we do not want to use low flux specific code, we disable it here.
 341   if (!fEnableLowFlux) lowFlux = false;
 342
 343   // Get FMD data
 344   AliESDFMD*  esdFMD  = esd->GetFMDData();
 345   AliMCEvent* mcEvent = MCEvent();
 346   // Apply the sharing filter (or hit merging or clustering if you like)
 347   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
 348     AliWarning("Sharing filter failed!");
 349     return;
 350   }
 351   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
 352     AliWarning("MC Sharing filter failed!");
 353     return;
 354   }
 355   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 356
 357   // Do the energy stuff
 358   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
 359     AliWarning("Energy fitter failed");
 360     return;
 361   }
 362
 363   // Calculate the inclusive charged particle density
 364   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 365     AliWarning("Density calculator failed!");
 366     return;
 367   }
 368   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
 369     AliWarning("MC Density calculator failed!");
 370     return;
 371   }
 372   fDensityCalculator.CompareResults(fHistos, fMCHistos);
 373
 374   // Do the secondary and other corrections.
 375   if (!fCorrections.Correct(fHistos, ivz)) {
 376     AliWarning("Corrections failed");
 377     return;
 378   }
 379   if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
 380     AliWarning("MC Corrections failed");
 381     return;
 382   }
 383   fCorrections.CompareResults(fHistos, fMCHistos);
 384
 385   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 386     AliWarning("Histogram collector failed");
 387     return;
 388   }
 389   if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
 390     AliWarning("MC Histogram collector failed");
 391     return;
 392   }
 393
 394   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 395     fHData->Add(&(fAODFMD.GetHistogram()));
 396
 397   PostData(1, fList);
 398 }
 399
 400 //____________________________________________________________________
 401 void
 402 AliForwardMCMultiplicityTask::Terminate(Option_t*)
 403 {
 404   //
 405   // End of job
 406   //
 407   // Parameters:
 408   //    option Not used
 409   //
 410   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 411   if (!list) {
 412     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 413     if (GetOutputData(1)) GetOutputData(1)->Print();
 414     return;
 415   }
 416
 417   // Get our histograms from the container
 418   TH1I* hEventsTr    = 0;
 419   TH1I* hEventsTrVtx = 0;
 420   TH1I* hTriggers    = 0;
 421   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 422                                        hEventsTrVtx, hTriggers)) {
 423     AliError(Form("Didn't get histograms from event selector "
 424                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 425                   hEventsTr, hEventsTrVtx));
 426     list->ls();
 427     return;
 428   }
 429
 430   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 431   if (!hData) {
 432     AliError(Form("Couldn't get our summed histogram from output "
 433                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 434     list->ls();
 435     return;
 436   }
 437
 438   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 439   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 440   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
 441   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 442   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 443   dNdeta->Divide(norm);
 444   dNdeta->SetStats(0);
 445   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 446                 "width");
 447   list->Add(dNdeta);
 448   list->Add(norm);
 449
 450   fEnergyFitter.Fit(list);
 451   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 452   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 453   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 454 }
 455
 456 //____________________________________________________________________
 457 void
 458 AliForwardMCMultiplicityTask::Print(Option_t* option) const
 459 {
 460   //
 461   // Print information
 462   //
 463   // Parameters:
 464   //    option Not used
 465   //
 466   AliForwardMultiplicityBase::Print(option);
 467   gROOT->IncreaseDirLevel();
 468   fEventInspector   .Print(option);
 469   fEnergyFitter     .Print(option);
 470   fSharingFilter    .Print(option);
 471   fDensityCalculator.Print(option);
 472   fCorrections      .Print(option);
 473   fHistCollector    .Print(option);
 474   gROOT->DecreaseDirLevel();
 475 }
 476
 477 //
 478 // EOF
 479 //