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