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     fESDFMD(),
  35     fMCESDFMD(),
  36     fMCHistos(),
  37     fMCAODFMD(),
  38     fMCRingSums(),
  39     fPrimary(0),
  40     fEventInspector(),
  41     fSharingFilter(),
  42     fDensityCalculator(),
  43     fCorrections(),
  44     fHistCollector(),
  45     fEventPlaneFinder()
  46 {
  47   //
  48   // Constructor
  49   //
  50 }
  51
  52 //____________________________________________________________________
  53 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  54   : AliForwardMultiplicityBase(name),
  55     fESDFMD(),
  56     fMCESDFMD(),
  57     fMCHistos(),
  58     fMCAODFMD(kTRUE),
  59     fMCRingSums(),
  60     fPrimary(0),
  61     fEventInspector("event"),
  62     fSharingFilter("sharing"),
  63     fDensityCalculator("density"),
  64     fCorrections("corrections"),
  65     fHistCollector("collector"),
  66     fEventPlaneFinder("eventplane")
  67 {
  68   //
  69   // Constructor
  70   //
  71   // Parameters:
  72   //    name Name of task
  73   //
  74 }
  75
  76 //____________________________________________________________________
  77 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
  78   : AliForwardMultiplicityBase(o),
  79     fESDFMD(o.fESDFMD),
  80     fMCESDFMD(o.fMCESDFMD),
  81     fMCHistos(o.fMCHistos),
  82     fMCAODFMD(o.fMCAODFMD),
  83     fMCRingSums(o.fMCRingSums),
  84     fPrimary(o.fPrimary),
  85     fEventInspector(o.fEventInspector),
  86     fSharingFilter(o.fSharingFilter),
  87     fDensityCalculator(o.fDensityCalculator),
  88     fCorrections(o.fCorrections),
  89     fHistCollector(o.fHistCollector),
  90     fEventPlaneFinder(o.fEventPlaneFinder)
  91 {
  92   //
  93   // Copy constructor
  94   //
  95   // Parameters:
  96   //    o Object to copy from
  97   //
  98 }
  99
 100 //____________________________________________________________________
 101 AliForwardMCMultiplicityTask&
 102 AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
 103 {
 104   //
 105   // Assignment operator
 106   //
 107   // Parameters:
 108   //    o Object to assign from
 109   //
 110   // Return:
 111   //    Reference to this object
 112   //
 113   if (&o == this) return *this;
 114   AliForwardMultiplicityBase::operator=(o);
 115
 116   fEventInspector    = o.fEventInspector;
 117   fSharingFilter     = o.fSharingFilter;
 118   fDensityCalculator = o.fDensityCalculator;
 119   fCorrections       = o.fCorrections;
 120   fHistCollector     = o.fHistCollector;
 121   fEventPlaneFinder  = o.fEventPlaneFinder;
 122   fMCHistos          = o.fMCHistos;
 123   fMCAODFMD          = o.fMCAODFMD;
 124   fMCRingSums        = o.fMCRingSums;
 125   fPrimary           = o.fPrimary;
 126   return *this;
 127 }
 128
 129 //____________________________________________________________________
 130 void
 131 AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
 132 {
 133   fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
 134   fCorrections.SetSecondaryForMC(!use);
 135 }
 136
 137 //____________________________________________________________________
 138 void
 139 AliForwardMCMultiplicityTask::CreateBranches(AliAODHandler* ah)
 140 {
 141   //
 142   // Create output objects
 143   //
 144   //
 145   AliForwardMultiplicityBase::CreateBranches(ah);
 146
 147   TObject* mcobj = &fMCAODFMD;
 148   ah->AddBranch("AliAODForwardMult", &mcobj);
 149
 150   fPrimary = new TH2D("primary", "MC Primaries", 1,0,1,20,0,TMath::TwoPi());
 151   fPrimary->SetXTitle("#eta");
 152   fPrimary->SetYTitle("#varphi [radians]");
 153   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 154   fPrimary->Sumw2();
 155   fPrimary->SetStats(0);
 156   fPrimary->SetDirectory(0);
 157
 158   ah->AddBranch("TH2D", &fPrimary);
 159 }
 160
 161
 162 //____________________________________________________________________
 163 void
 164 AliForwardMCMultiplicityTask::InitMembers(const TAxis* pe, const TAxis* pv)
 165 {
 166   //
 167   // Initialise the sub objects and stuff.  Called on first event
 168   //
 169   //
 170   AliForwardMultiplicityBase::InitMembers(pe, pv);
 171
 172   fMCHistos.Init(*pe);
 173   fMCAODFMD.Init(*pe);
 174   fMCRingSums.Init(*pe);
 175
 176   AliForwardUtil::Histos::RebinEta(fPrimary, *pe);
 177
 178   TList* mcRings = new TList;
 179   mcRings->SetName("mcRingSums");
 180   mcRings->SetOwner();
 181   fList->Add(mcRings);
 182
 183   mcRings->Add(fMCRingSums.Get(1, 'I'));
 184   mcRings->Add(fMCRingSums.Get(2, 'I'));
 185   mcRings->Add(fMCRingSums.Get(2, 'O'));
 186   mcRings->Add(fMCRingSums.Get(3, 'I'));
 187   mcRings->Add(fMCRingSums.Get(3, 'O'));
 188   fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 189   fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 190   fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 191   fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 192   fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 193 }
 194
 195 //____________________________________________________________________
 196 void
 197 AliForwardMCMultiplicityTask::UserExec(Option_t*)
 198 {
 199   //
 200   // Process each event
 201   //
 202   // Parameters:
 203   //    option Not used
 204   //
 205
 206   // Read production details
 207   if (fFirstEvent)
 208     fEventInspector.ReadProductionDetails(MCEvent());
 209
 210   // Get the input data
 211   AliESDEvent* esd     = GetESDEvent();
 212   AliMCEvent*  mcEvent = MCEvent();
 213   if (!esd || !mcEvent) return;
 214
 215   // Clear stuff
 216   fHistos.Clear();
 217   fESDFMD.Clear();
 218   fAODFMD.Clear();
 219   fAODEP.Clear();
 220   fMCHistos.Clear();
 221   fMCESDFMD.Clear();
 222   fMCAODFMD.Clear();
 223   fPrimary->Reset();
 224
 225   Bool_t   lowFlux   = kFALSE;
 226   UInt_t   triggers  = 0;
 227   UShort_t ivz       = 0;
 228   TVector3 ip(1024, 1024, 0);
 229   Double_t cent      = -1;
 230   UShort_t nClusters = 0;
 231   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 232                                                ivz, ip, cent, nClusters);
 233   UShort_t ivzMC    = 0;
 234   Double_t vzMC     = 0;
 235   Double_t phiR     = 0;
 236   Double_t b        = 0;
 237   Double_t cMC      = 0;
 238   Int_t    npart    = 0;
 239   Int_t    nbin     = 0;
 240   // UInt_t   foundMC  =
 241   fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
 242                             npart, nbin, phiR);
 243   fEventInspector.CompareResults(ip.Z(), vzMC, cent, cMC, b, npart, nbin);
 244
 245   //Store all events
 246   MarkEventForStore();
 247
 248   Bool_t isAccepted = true;
 249   if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
 250   if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
 251   //MarkEventForStore();
 252   // Always set the B trigger - each MC event _is_ a collision
 253   triggers |= AliAODForwardMult::kB;
 254   // Set trigger bits, and mark this event for storage
 255   fAODFMD.SetTriggerBits(triggers);
 256   fAODFMD.SetSNN(fEventInspector.GetEnergy());
 257   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 258   fAODFMD.SetCentrality(cent);
 259   fAODFMD.SetNClusters(nClusters);
 260
 261   fMCAODFMD.SetTriggerBits(triggers);
 262   fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
 263   fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 264   fMCAODFMD.SetCentrality(cent);
 265   fMCAODFMD.SetNClusters(nClusters);
 266
 267   //All events should be stored - HHD
 268   //if (isAccepted) MarkEventForStore();
 269
 270   // Disable this check on SPD - will bias data
 271   // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
 272   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
 273   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
 274
 275   if (isAccepted) {
 276     fAODFMD.SetIpZ(ip.Z());
 277     fMCAODFMD.SetIpZ(ip.Z());
 278   }
 279   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 280
 281   // We we do not want to use low flux specific code, we disable it here.
 282   if (!fEnableLowFlux) lowFlux = false;
 283
 284
 285
 286   // Get FMD data
 287   AliESDFMD*  esdFMD  = esd->GetFMDData();
 288
 289   // Apply the sharing filter (or hit merging or clustering if you like)
 290   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())) {
 291     AliWarning("Sharing filter failed!");
 292     return;
 293   }
 294   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),fMCESDFMD,fPrimary)) {
 295     AliWarning("MC Sharing filter failed!");
 296     return;
 297   }
 298
 299   // Store some MC parameters in corners of histogram :-)
 300   fPrimary->SetBinContent(0,                      0,                    vzMC);
 301   fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,0,                    phiR);
 302   fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,fPrimary->GetNbinsY(),cMC);
 303
 304
 305   if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
 306   // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
 307
 308   //MarkEventForStore();
 309   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 310
 311   // Calculate the inclusive charged particle density
 312   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
 313     AliWarning("Density calculator failed!");
 314     return;
 315   }
 316   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
 317     AliWarning("MC Density calculator failed!");
 318     return;
 319   }
 320   fDensityCalculator.CompareResults(fHistos, fMCHistos);
 321
 322   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
 323     if (!fEventPlaneFinder.FindEventplane(esd, fAODEP,
 324                                           &(fAODFMD.GetHistogram()) , &fHistos))
 325       AliWarning("Eventplane finder failed!");
 326   }
 327
 328   // Do the secondary and other corrections.
 329   if (!fCorrections.Correct(fHistos, ivz)) {
 330     AliWarning("Corrections failed");
 331     return;
 332   }
 333   if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
 334     AliWarning("MC Corrections failed");
 335     return;
 336   }
 337   fCorrections.CompareResults(fHistos, fMCHistos);
 338
 339   if (!fHistCollector.Collect(fHistos, fRingSums,
 340                               ivz, fAODFMD.GetHistogram(),
 341                               fAODFMD.GetCentrality())) {
 342     AliWarning("Histogram collector failed");
 343     return;
 344   }
 345   if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
 346                               ivz, fMCAODFMD.GetHistogram())) {
 347     AliWarning("MC Histogram collector failed");
 348     return;
 349   }
 350   // Copy underflow bins to overflow bins - always full phi coverage
 351   TH2&  hMC  = fMCAODFMD.GetHistogram();
 352   Int_t nEta = hMC.GetNbinsX();
 353   Int_t nY   = hMC.GetNbinsY();
 354   for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 355     hMC.SetBinContent(iEta, nY+1, hMC.GetBinContent(iEta, 0));
 356   }
 357
 358   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 359     fHData->Add(&(fAODFMD.GetHistogram()));
 360
 361   PostData(1, fList);
 362 }
 363
 364 //____________________________________________________________________
 365 void
 366 AliForwardMCMultiplicityTask::EstimatedNdeta(const TList* input,
 367                                              TList*       output) const
 368 {
 369   AliForwardMultiplicityBase::EstimatedNdeta(input, output);
 370   MakeRingdNdeta(input, "mcRingSums", output, "mcRingResults", 24);
 371 }
 372
 373
 374
 375 //
 376 // EOF
 377 //