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