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   UInt_t what = AliForwardCorrectionManager::kAll;
 171   if (!fEnableLowFlux)
 172     what ^= AliForwardCorrectionManager::kDoubleHit;
 173   if (!fCorrections.IsUseMergingEfficiency())
 174     what ^= AliForwardCorrectionManager::kMergingEfficiency;
 175
 176   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 177   fcm.Init(fEventInspector.GetCollisionSystem(),
 178            fEventInspector.GetEnergy(),
 179            fEventInspector.GetField(),
 180            true,
 181            what);
 182   if (!CheckCorrections(what)) return;
 183
 184   const TAxis* pe = fcm.GetEtaAxis();
 185   const TAxis* pv = fcm.GetVertexAxis();
 186   if (!pe) AliFatal("No eta axis defined");
 187   if (!pv) AliFatal("No vertex axis defined");
 188
 189   fHistos.Init(*pe);
 190   fAODFMD.Init(*pe);
 191   fMCHistos.Init(*pe);
 192   fMCAODFMD.Init(*pe);
 193
 194   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 195   fHData->SetStats(0);
 196   fHData->SetDirectory(0);
 197
 198   fList->Add(fHData);
 199
 200   fEnergyFitter.Init(*pe);
 201   fEventInspector.Init(*pv);
 202   fDensityCalculator.Init(*pe);
 203   fCorrections.Init(*pe);
 204   fHistCollector.Init(*pv);
 205
 206   this->Print();
 207 }
 208
 209 //____________________________________________________________________
 210 void
 211 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
 212 {
 213   //
 214   // Create output objects
 215   //
 216   //
 217   fList = new TList;
 218
 219   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 220   AliAODHandler*      ah =
 221     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 222   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 223
 224
 225   TObject* obj = &fAODFMD;
 226   ah->AddBranch("AliAODForwardMult", &obj);
 227
 228   TObject* mcobj = &fMCAODFMD;
 229   ah->AddBranch("AliAODForwardMult", &mcobj);
 230
 231   fPrimary = new TH2D("primary", "MC Primaries",
 232                       200, -4, 6, 20, 0, 2*TMath::Pi());
 233   fPrimary->SetXTitle("#eta");
 234   fPrimary->SetYTitle("#varphi [radians]");
 235   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 236   fPrimary->Sumw2();
 237   fPrimary->SetStats(0);
 238   fPrimary->SetDirectory(0);
 239   ah->AddBranch("TH2D", &fPrimary);
 240
 241   fEventInspector.DefineOutput(fList);
 242   fEnergyFitter.DefineOutput(fList);
 243   fSharingFilter.DefineOutput(fList);
 244   fDensityCalculator.DefineOutput(fList);
 245   fCorrections.DefineOutput(fList);
 246
 247   PostData(1, fList);
 248 }
 249 //____________________________________________________________________
 250 void
 251 AliForwardMCMultiplicityTask::UserExec(Option_t*)
 252 {
 253   //
 254   // Process each event
 255   //
 256   // Parameters:
 257   //    option Not used
 258   //
 259
 260   // Get the input data
 261   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 262   if (!esd) {
 263     AliWarning("No ESD event found for input event");
 264     return;
 265   }
 266
 267   // On the first event, initialize the parameters
 268   if (fFirstEvent && esd->GetESDRun()) {
 269     fEventInspector.ReadRunDetails(esd);
 270
 271     AliInfo(Form("Initializing with parameters from the ESD:\n"
 272                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 273                  "         AliESDEvent::GetBeamType()     ->%s\n"
 274                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 275                  "         AliESDEvent::GetMagneticField()->%f\n"
 276                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 277                  esd->GetBeamEnergy(),
 278                  esd->GetBeamType(),
 279                  esd->GetCurrentL3(),
 280                  esd->GetMagneticField(),
 281                  esd->GetRunNumber()));
 282
 283     fFirstEvent = false;
 284
 285     InitializeSubs();
 286   }
 287   // Clear stuff
 288   fHistos.Clear();
 289   fESDFMD.Clear();
 290   fAODFMD.Clear();
 291   fMCHistos.Clear();
 292   fMCESDFMD.Clear();
 293   fMCAODFMD.Clear();
 294   fPrimary->Reset();
 295
 296   Bool_t   lowFlux  = kFALSE;
 297   UInt_t   triggers = 0;
 298   UShort_t ivz      = 0;
 299   Double_t vz       = 0;
 300   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
 301   //Store all events
 302   MarkEventForStore();
 303
 304   Bool_t isAccepted = true;
 305   if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
 306   if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
 307
 308   // Set trigger bits, and mark this event for storage
 309   fAODFMD.SetTriggerBits(triggers);
 310   fMCAODFMD.SetTriggerBits(triggers);
 311
 312   //All events should be stored - HHD
 313   //MarkEventForStore();
 314
 315   if (found & AliFMDEventInspector::kNoSPD)     isAccepted = false; // return;
 316   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
 317   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
 318
 319   if (isAccepted) {
 320     fAODFMD.SetIpZ(vz);
 321     fMCAODFMD.SetIpZ(vz);
 322   }
 323   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 324
 325   // We we do not want to use low flux specific code, we disable it here.
 326   if (!fEnableLowFlux) lowFlux = false;
 327
 328   // Get FMD data
 329   AliESDFMD*  esdFMD  = esd->GetFMDData();
 330   AliMCEvent* mcEvent = MCEvent();
 331   // Apply the sharing filter (or hit merging or clustering if you like)
 332   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
 333     AliWarning("Sharing filter failed!");
 334     return;
 335   }
 336   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
 337     AliWarning("MC Sharing filter failed!");
 338     return;
 339   }
 340   if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
 341   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 342
 343   // Do the energy stuff
 344   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
 345     AliWarning("Energy fitter failed");
 346     return;
 347   }
 348
 349   // Calculate the inclusive charged particle density
 350   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 351     AliWarning("Density calculator failed!");
 352     return;
 353   }
 354   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
 355     AliWarning("MC Density calculator failed!");
 356     return;
 357   }
 358   fDensityCalculator.CompareResults(fHistos, fMCHistos);
 359
 360   // Do the secondary and other corrections.
 361   if (!fCorrections.Correct(fHistos, ivz)) {
 362     AliWarning("Corrections failed");
 363     return;
 364   }
 365   if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
 366     AliWarning("MC Corrections failed");
 367     return;
 368   }
 369   fCorrections.CompareResults(fHistos, fMCHistos);
 370
 371   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 372     AliWarning("Histogram collector failed");
 373     return;
 374   }
 375   if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
 376     AliWarning("MC Histogram collector failed");
 377     return;
 378   }
 379
 380   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 381     fHData->Add(&(fAODFMD.GetHistogram()));
 382
 383   PostData(1, fList);
 384 }
 385
 386 //____________________________________________________________________
 387 void
 388 AliForwardMCMultiplicityTask::Terminate(Option_t*)
 389 {
 390   //
 391   // End of job
 392   //
 393   // Parameters:
 394   //    option Not used
 395   //
 396   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 397   if (!list) {
 398     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 399     if (GetOutputData(1)) GetOutputData(1)->Print();
 400     return;
 401   }
 402
 403   // Get our histograms from the container
 404   TH1I* hEventsTr    = 0;
 405   TH1I* hEventsTrVtx = 0;
 406   TH1I* hTriggers    = 0;
 407   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 408                                        hEventsTrVtx, hTriggers)) {
 409     AliError(Form("Didn't get histograms from event selector "
 410                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 411                   hEventsTr, hEventsTrVtx));
 412     list->ls();
 413     return;
 414   }
 415
 416   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 417   if (!hData) {
 418     AliError(Form("Couldn't get our summed histogram from output "
 419                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 420     list->ls();
 421     return;
 422   }
 423
 424   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 425   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 426   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
 427   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 428   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 429   dNdeta->Divide(norm);
 430   dNdeta->SetStats(0);
 431   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 432                 "width");
 433   list->Add(dNdeta);
 434   list->Add(norm);
 435
 436   fEnergyFitter.Fit(list);
 437   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 438   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 439   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 440 }
 441
 442
 443 //
 444 // EOF
 445 //