PWGLF/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     fAODEP(),
  39     fMCESDFMD(),
  40     fMCHistos(),
  41     fMCAODFMD(),
  42     fRingSums(),
  43     fMCRingSums(),
  44     fPrimary(0),
  45     fEventInspector(),
  46     fSharingFilter(),
  47     fDensityCalculator(),
  48     fCorrections(),
  49     fHistCollector(),
  50     fEventPlaneFinder(),
  51     fList(0)
  52 {
  53   //
  54   // Constructor
  55   //
  56 }
  57
  58 //____________________________________________________________________
  59 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  60   : AliForwardMultiplicityBase(name),
  61     fHData(0),
  62     fESDFMD(),
  63     fHistos(),
  64     fAODFMD(kFALSE),
  65     fAODEP(kFALSE),
  66     fMCESDFMD(),
  67     fMCHistos(),
  68     fMCAODFMD(kTRUE),
  69     fRingSums(),
  70     fMCRingSums(),
  71     fPrimary(0),
  72     fEventInspector("event"),
  73     fSharingFilter("sharing"),
  74     fDensityCalculator("density"),
  75     fCorrections("corrections"),
  76     fHistCollector("collector"),
  77     fEventPlaneFinder("eventplane"),
  78     fList(0)
  79 {
  80   //
  81   // Constructor
  82   //
  83   // Parameters:
  84   //    name Name of task
  85   //
  86   DefineOutput(1, TList::Class());
  87 }
  88
  89 //____________________________________________________________________
  90 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
  91   : AliForwardMultiplicityBase(o),
  92     fHData(o.fHData),
  93     fESDFMD(o.fESDFMD),
  94     fHistos(o.fHistos),
  95     fAODFMD(o.fAODFMD),
  96     fAODEP(o.fAODEP),
  97     fMCESDFMD(o.fMCESDFMD),
  98     fMCHistos(o.fMCHistos),
  99     fMCAODFMD(o.fMCAODFMD),
 100     fRingSums(o.fRingSums),
 101     fMCRingSums(o.fMCRingSums),
 102     fPrimary(o.fPrimary),
 103     fEventInspector(o.fEventInspector),
 104     fSharingFilter(o.fSharingFilter),
 105     fDensityCalculator(o.fDensityCalculator),
 106     fCorrections(o.fCorrections),
 107     fHistCollector(o.fHistCollector),
 108     fEventPlaneFinder(o.fEventPlaneFinder),
 109     fList(o.fList)
 110 {
 111   //
 112   // Copy constructor
 113   //
 114   // Parameters:
 115   //    o Object to copy from
 116   //
 117   DefineOutput(1, TList::Class());
 118 }
 119
 120 //____________________________________________________________________
 121 AliForwardMCMultiplicityTask&
 122 AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
 123 {
 124   //
 125   // Assignment operator
 126   //
 127   // Parameters:
 128   //    o Object to assign from
 129   //
 130   // Return:
 131   //    Reference to this object
 132   //
 133   if (&o == this) return *this;
 134   AliForwardMultiplicityBase::operator=(o);
 135
 136   fHData             = o.fHData;
 137   fEventInspector    = o.fEventInspector;
 138   fSharingFilter     = o.fSharingFilter;
 139   fDensityCalculator = o.fDensityCalculator;
 140   fCorrections       = o.fCorrections;
 141   fHistCollector     = o.fHistCollector;
 142   fEventPlaneFinder  = o.fEventPlaneFinder;
 143   fHistos            = o.fHistos;
 144   fAODFMD            = o.fAODFMD;
 145   fAODEP             = o.fAODEP;
 146   fMCHistos          = o.fMCHistos;
 147   fMCAODFMD          = o.fMCAODFMD;
 148   fRingSums          = o.fRingSums;
 149   fMCRingSums        = o.fMCRingSums;
 150   fPrimary           = o.fPrimary;
 151   fList              = o.fList;
 152
 153   return *this;
 154 }
 155
 156 //____________________________________________________________________
 157 void
 158 AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
 159 {
 160   //
 161   // Set debug level
 162   //
 163   // Parameters:
 164   //    dbg debug level
 165   //
 166   fEventInspector.SetDebug(dbg);
 167   fSharingFilter.SetDebug(dbg);
 168   fDensityCalculator.SetDebug(dbg);
 169   fCorrections.SetDebug(dbg);
 170   fHistCollector.SetDebug(dbg);
 171   fEventPlaneFinder.SetDebug(dbg);
 172 }
 173 //____________________________________________________________________
 174 void
 175 AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
 176 {
 177   fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
 178   fCorrections.SetSecondaryForMC(!use);
 179 }
 180
 181 //____________________________________________________________________
 182 void
 183 AliForwardMCMultiplicityTask::InitializeSubs()
 184 {
 185   //
 186   // Initialise the sub objects and stuff.  Called on first event
 187   //
 188   //
 189   const TAxis* pe = 0;
 190   const TAxis* pv = 0;
 191
 192   if (!ReadCorrections(pe,pv,true)) return;
 193
 194   fHistos.Init(*pe);
 195   fAODFMD.Init(*pe);
 196   fAODEP.Init(*pe);
 197   fMCHistos.Init(*pe);
 198   fMCAODFMD.Init(*pe);
 199   fRingSums.Init(*pe);
 200   fMCRingSums.Init(*pe);
 201
 202   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 203   fHData->SetStats(0);
 204   fHData->SetDirectory(0);
 205
 206   fList->Add(fHData);
 207
 208   TList* rings = new TList;
 209   rings->SetName("ringSums");
 210   rings->SetOwner();
 211   fList->Add(rings);
 212
 213   rings->Add(fRingSums.Get(1, 'I'));
 214   rings->Add(fRingSums.Get(2, 'I'));
 215   rings->Add(fRingSums.Get(2, 'O'));
 216   rings->Add(fRingSums.Get(3, 'I'));
 217   rings->Add(fRingSums.Get(3, 'O'));
 218   fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 219   fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 220   fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 221   fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 222   fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 223
 224   TList* mcRings = new TList;
 225   mcRings->SetName("mcRingSums");
 226   mcRings->SetOwner();
 227   fList->Add(mcRings);
 228
 229   mcRings->Add(fMCRingSums.Get(1, 'I'));
 230   mcRings->Add(fMCRingSums.Get(2, 'I'));
 231   mcRings->Add(fMCRingSums.Get(2, 'O'));
 232   mcRings->Add(fMCRingSums.Get(3, 'I'));
 233   mcRings->Add(fMCRingSums.Get(3, 'O'));
 234   fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 235   fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 236   fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 237   fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 238   fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 239
 240
 241
 242   fEventInspector.Init(*pv);
 243   fSharingFilter.Init();
 244   fDensityCalculator.Init(*pe);
 245   fCorrections.Init(*pe);
 246   fHistCollector.Init(*pv,*pe);
 247   fEventPlaneFinder.Init(*pe);
 248
 249   this->Print();
 250 }
 251
 252 //____________________________________________________________________
 253 void
 254 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
 255 {
 256   //
 257   // Create output objects
 258   //
 259   //
 260   fList = new TList;
 261   fList->SetOwner();
 262
 263   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 264   AliAODHandler*      ah =
 265     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 266   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 267
 268
 269   TObject* obj = &fAODFMD;
 270   ah->AddBranch("AliAODForwardMult", &obj);
 271
 272   TObject* mcobj = &fMCAODFMD;
 273   ah->AddBranch("AliAODForwardMult", &mcobj);
 274
 275   TObject* epobj = &fAODFMD;
 276   ah->AddBranch("AliAODForwardEP", &epobj);
 277
 278   fPrimary = new TH2D("primary", "MC Primaries",
 279                       200, -4, 6, 20, 0, 2*TMath::Pi());
 280   fPrimary->SetXTitle("#eta");
 281   fPrimary->SetYTitle("#varphi [radians]");
 282   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 283   fPrimary->Sumw2();
 284   fPrimary->SetStats(0);
 285   fPrimary->SetDirectory(0);
 286   ah->AddBranch("TH2D", &fPrimary);
 287
 288   fEventInspector.DefineOutput(fList);
 289   fSharingFilter.DefineOutput(fList);
 290   fDensityCalculator.DefineOutput(fList);
 291   fCorrections.DefineOutput(fList);
 292   fHistCollector.DefineOutput(fList);
 293   fEventPlaneFinder.DefineOutput(fList);
 294
 295   PostData(1, fList);
 296 }
 297 //____________________________________________________________________
 298 void
 299 AliForwardMCMultiplicityTask::UserExec(Option_t*)
 300 {
 301   //
 302   // Process each event
 303   //
 304   // Parameters:
 305   //    option Not used
 306   //
 307
 308   // Read production details
 309   if (fFirstEvent)
 310     fEventInspector.ReadProductionDetails(MCEvent());
 311
 312   // Get the input data
 313   AliESDEvent* esd     = GetESDEvent();
 314   AliMCEvent*  mcEvent = MCEvent();
 315
 316   // Clear stuff
 317   fHistos.Clear();
 318   fESDFMD.Clear();
 319   fAODFMD.Clear();
 320   fAODEP.Clear();
 321   fMCHistos.Clear();
 322   fMCESDFMD.Clear();
 323   fMCAODFMD.Clear();
 324   fPrimary->Reset();
 325
 326   Bool_t   lowFlux   = kFALSE;
 327   UInt_t   triggers  = 0;
 328   UShort_t ivz       = 0;
 329   Double_t vz        = 0;
 330   Double_t cent      = -1;
 331   UShort_t nClusters = 0;
 332   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 333                                                ivz, vz, cent, nClusters);
 334   UShort_t ivzMC    = 0;
 335   Double_t vzMC     = 0;
 336   Double_t phiR     = 0;
 337   Double_t b        = 0;
 338   Int_t    npart    = 0;
 339   Int_t    nbin     = 0;
 340   // UInt_t   foundMC  =
 341   fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
 342                             npart, nbin, phiR);
 343   fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
 344
 345   //Store all events
 346   MarkEventForStore();
 347
 348   Bool_t isAccepted = true;
 349   if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
 350   if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
 351   //MarkEventForStore();
 352   // Always set the B trigger - each MC event _is_ a collision
 353   triggers |= AliAODForwardMult::kB;
 354   // Set trigger bits, and mark this event for storage
 355   fAODFMD.SetTriggerBits(triggers);
 356   fAODFMD.SetSNN(fEventInspector.GetEnergy());
 357   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 358   fAODFMD.SetCentrality(cent);
 359   fAODFMD.SetNClusters(nClusters);
 360
 361   fMCAODFMD.SetTriggerBits(triggers);
 362   fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
 363   fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 364   fMCAODFMD.SetCentrality(cent);
 365   fMCAODFMD.SetNClusters(nClusters);
 366
 367   //All events should be stored - HHD
 368   //if (isAccepted) MarkEventForStore();
 369
 370   // Disable this check on SPD - will bias data
 371   // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
 372   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
 373   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
 374
 375   if (isAccepted) {
 376     fAODFMD.SetIpZ(vz);
 377     fMCAODFMD.SetIpZ(vz);
 378   }
 379   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 380
 381   // We we do not want to use low flux specific code, we disable it here.
 382   if (!fEnableLowFlux) lowFlux = false;
 383
 384
 385
 386   // Get FMD data
 387   AliESDFMD*  esdFMD  = esd->GetFMDData();
 388
 389   // Apply the sharing filter (or hit merging or clustering if you like)
 390   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 391     AliWarning("Sharing filter failed!");
 392     return;
 393   }
 394   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
 395     AliWarning("MC Sharing filter failed!");
 396     return;
 397   }
 398   if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
 399   // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
 400
 401   //MarkEventForStore();
 402   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 403
 404   // Calculate the inclusive charged particle density
 405   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
 406     AliWarning("Density calculator failed!");
 407     return;
 408   }
 409   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
 410     AliWarning("MC Density calculator failed!");
 411     return;
 412   }
 413   fDensityCalculator.CompareResults(fHistos, fMCHistos);
 414
 415   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
 416     if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
 417       AliWarning("Eventplane finder failed!");
 418   }
 419
 420   // Do the secondary and other corrections.
 421   if (!fCorrections.Correct(fHistos, ivz)) {
 422     AliWarning("Corrections failed");
 423     return;
 424   }
 425   if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
 426     AliWarning("MC Corrections failed");
 427     return;
 428   }
 429   fCorrections.CompareResults(fHistos, fMCHistos);
 430
 431   if (!fHistCollector.Collect(fHistos, fRingSums,
 432                               ivz, fAODFMD.GetHistogram())) {
 433     AliWarning("Histogram collector failed");
 434     return;
 435   }
 436   if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
 437                               ivz, fMCAODFMD.GetHistogram())) {
 438     AliWarning("MC Histogram collector failed");
 439     return;
 440   }
 441
 442   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 443     fHData->Add(&(fAODFMD.GetHistogram()));
 444
 445   PostData(1, fList);
 446 }
 447
 448 //____________________________________________________________________
 449 void
 450 AliForwardMCMultiplicityTask::Terminate(Option_t*)
 451 {
 452   //
 453   // End of job
 454   //
 455   // Parameters:
 456   //    option Not used
 457   //
 458   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 459   if (!list) {
 460     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 461     if (GetOutputData(1)) GetOutputData(1)->Print();
 462     return;
 463   }
 464
 465   // Get our histograms from the container
 466   TH1I* hEventsTr    = 0;
 467   TH1I* hEventsTrVtx = 0;
 468   TH1I* hTriggers    = 0;
 469   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 470                                        hEventsTrVtx, hTriggers)) {
 471     AliError(Form("Didn't get histograms from event selector "
 472                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 473                   hEventsTr, hEventsTrVtx));
 474     list->ls();
 475     return;
 476   }
 477
 478   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 479   if (!hData) {
 480     AliError(Form("Couldn't get our summed histogram from output "
 481                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 482     list->ls();
 483     return;
 484   }
 485
 486   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 487   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 488   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
 489   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 490   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 491   dNdeta->Divide(norm);
 492   dNdeta->SetStats(0);
 493   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 494                 "width");
 495   list->Add(dNdeta);
 496   list->Add(norm);
 497
 498   MakeRingdNdeta(list, "ringSums", list, "ringResults");
 499   MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
 500
 501   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
 502   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 503   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 504 }
 505
 506
 507 //
 508 // EOF
 509 //