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