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 #define MCAOD_SLOT 4
  31 #define PRIMARY_SLOT 5
  32 #ifdef POST_AOD
  33 # define DEFINE(N,C) DefineOutput(N,C)
  34 # define POST(N,O)   PostData(N,O)
  35 #else
  36 # define DEFINE(N,C) do { } while(false)
  37 # define POST(N,O)   do { } while(false)
  38 #endif
  39 #ifndef ENABLE_TIMING
  40 # define MAKE_SW(NAME) do {} while(false)
  41 # define START_SW(NAME) do {} while(false)
  42 # define FILL_SW(NAME,WHICH) do {} while(false)
  43 #else
  44 # define MAKE_SW(NAME) TStopwatch NAME
  45 # define START_SW(NAME) if (fDoTiming) NAME.Start(true)
  46 # define FILL_SW(NAME,WHICH)                            \
  47   if (fDoTiming) fHTiming->Fill(WHICH,NAME.CpuTime())
  48 #endif
  49
  50 //====================================================================
  51 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
  52   : AliForwardMultiplicityBase(),
  53     fESDFMD(),
  54     fMCESDFMD(),
  55     fMCHistos(),
  56     fMCAODFMD(),
  57     fMCRingSums(),
  58     fPrimary(0),
  59     fEventInspector(),
  60     fESDFixer(),
  61     fSharingFilter(),
  62     fDensityCalculator(),
  63     fCorrections(),
  64     fHistCollector(),
  65     fEventPlaneFinder()
  66 {
  67   //
  68   // Constructor
  69   //
  70 }
  71
  72 //____________________________________________________________________
  73 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  74   : AliForwardMultiplicityBase(name),
  75     fESDFMD(),
  76     fMCESDFMD(),
  77     fMCHistos(),
  78     fMCAODFMD(kTRUE),
  79     fMCRingSums(),
  80     fPrimary(0),
  81     fEventInspector("event"),
  82     fESDFixer("esdFizer"),
  83     fSharingFilter("sharing"),
  84     fDensityCalculator("density"),
  85     fCorrections("corrections"),
  86     fHistCollector("collector"),
  87     fEventPlaneFinder("eventplane")
  88 {
  89   //
  90   // Constructor
  91   //
  92   // Parameters:
  93   //    name Name of task
  94   //
  95   fPrimary = new TH2D("primary", "MC Primaries", 1,0,1,20,0,TMath::TwoPi());
  96   fPrimary->SetXTitle("#eta");
  97   fPrimary->SetYTitle("#varphi [radians]");
  98   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
  99   fPrimary->Sumw2();
 100   fPrimary->SetStats(0);
 101   fPrimary->SetDirectory(0);
 102   DEFINE(MCAOD_SLOT,AliAODForwardMult::Class());
 103   DEFINE(PRIM_SLOT, TH2D::Class());
 104 }
 105
 106 //____________________________________________________________________
 107 void
 108 AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
 109 {
 110   fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
 111   fCorrections.SetSecondaryForMC(!use);
 112 }
 113
 114 //____________________________________________________________________
 115 void
 116 AliForwardMCMultiplicityTask::CreateBranches(AliAODHandler* ah)
 117 {
 118   //
 119   // Create output objects
 120   //
 121   //
 122   AliForwardMultiplicityBase::CreateBranches(ah);
 123
 124   TObject* mcobj = &fMCAODFMD;
 125   ah->AddBranch("AliAODForwardMult", &mcobj);
 126   ah->AddBranch("TH2D", &fPrimary);
 127 }
 128
 129 //____________________________________________________________________
 130 Bool_t
 131 AliForwardMCMultiplicityTask::Book()
 132 {
 133   // We do this to explicitly disable the noise corrector for MC
 134   GetESDFixer().SetRecoNoiseFactor(5);
 135
 136   Bool_t ret = AliForwardMultiplicityBase::Book();
 137   POST(MCAOD_SLOT, &fMCAODFMD);
 138   POST(PRIM_SLOT,  fPrimary);
 139   return ret;
 140 }
 141
 142 //____________________________________________________________________
 143 void
 144 AliForwardMCMultiplicityTask::InitMembers(const TAxis& eta, const TAxis& vertex)
 145 {
 146   //
 147   // Initialise the sub objects and stuff.  Called on first event
 148   //
 149   //
 150   AliForwardMultiplicityBase::InitMembers(eta, vertex);
 151
 152   fMCHistos.Init(eta);
 153   fMCAODFMD.Init(eta);
 154   fMCRingSums.Init(eta);
 155
 156   AliForwardUtil::Histos::RebinEta(fPrimary, eta);
 157   DMSG(fDebug,0,"Primary histogram rebinned to %d,%f,%f eta axis %d,%f,%f",
 158        fPrimary->GetXaxis()->GetNbins(),
 159        fPrimary->GetXaxis()->GetXmin(),
 160        fPrimary->GetXaxis()->GetXmax(),
 161        eta.GetNbins(),
 162        eta.GetXmin(),
 163        eta.GetXmax());
 164
 165
 166   TList* mcRings = new TList;
 167   mcRings->SetName("mcRingSums");
 168   mcRings->SetOwner();
 169   fList->Add(mcRings);
 170
 171   mcRings->Add(fMCRingSums.Get(1, 'I'));
 172   mcRings->Add(fMCRingSums.Get(2, 'I'));
 173   mcRings->Add(fMCRingSums.Get(2, 'O'));
 174   mcRings->Add(fMCRingSums.Get(3, 'I'));
 175   mcRings->Add(fMCRingSums.Get(3, 'O'));
 176   fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 177   fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 178   fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 179   fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 180   fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 181 }
 182
 183 //____________________________________________________________________
 184 Bool_t
 185 AliForwardMCMultiplicityTask::PreEvent()
 186 {
 187  if (fFirstEvent)
 188     fEventInspector.ReadProductionDetails(MCEvent());
 189   // Clear stuff
 190   fHistos.Clear();
 191   fESDFMD.Clear();
 192   fAODFMD.Clear();
 193   fAODEP.Clear();
 194   fMCHistos.Clear();
 195   fMCESDFMD.Clear();
 196   fMCAODFMD.Clear();
 197   fPrimary->Reset();
 198   return true;
 199 }
 200 //____________________________________________________________________
 201 Bool_t
 202 AliForwardMCMultiplicityTask::Event(AliESDEvent& esd)
 203 {
 204   //
 205   // Process each event
 206   //
 207   // Parameters:
 208   //    option Not used
 209   //
 210   MAKE_SW(total);
 211   MAKE_SW(individual);
 212   START_SW(total);
 213   START_SW(individual);
 214
 215   // Read production details
 216
 217   // Get the input data
 218   AliMCEvent*  mcEvent = MCEvent();
 219   if (!mcEvent) return false;
 220
 221   Bool_t   lowFlux   = kFALSE;
 222   UInt_t   triggers  = 0;
 223   UShort_t ivz       = 0;
 224   TVector3 ip(1024, 1024, 0);
 225   Double_t cent      = -1;
 226   UShort_t nClusters = 0;
 227   UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux,
 228                                                ivz, ip, cent, nClusters);
 229   UShort_t ivzMC    = 0;
 230   Double_t vzMC     = 0;
 231   Double_t phiR     = 0;
 232   Double_t b        = 0;
 233   Double_t cMC      = 0;
 234   Int_t    npart    = 0;
 235   Int_t    nbin     = 0;
 236   // UInt_t   foundMC  =
 237   fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
 238                             npart, nbin, phiR);
 239   fEventInspector.CompareResults(ip.Z(), vzMC, cent, cMC, b, npart, nbin);
 240   FILL_SW(individual,kTimingEventInspector);
 241
 242   //Store all events
 243   MarkEventForStore();
 244
 245   Bool_t isAccepted = true;
 246   if (found & AliFMDEventInspector::kNoEvent)    {
 247     fHStatus->Fill(1);
 248     isAccepted = false;
 249     // return;
 250   }
 251   if (found & AliFMDEventInspector::kNoTriggers) {
 252     fHStatus->Fill(2);
 253     isAccepted = false;
 254     // return;
 255   }
 256   //MarkEventForStore();
 257   // Always set the B trigger - each MC event _is_ a collision
 258   triggers |= AliAODForwardMult::kB;
 259   // Set trigger bits, and mark this event for storage
 260   fAODFMD.SetTriggerBits(triggers);
 261   fAODFMD.SetSNN(fEventInspector.GetEnergy());
 262   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 263   fAODFMD.SetCentrality(cent);
 264   fAODFMD.SetNClusters(nClusters);
 265
 266   fMCAODFMD.SetTriggerBits(triggers);
 267   fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
 268   fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
 269   fMCAODFMD.SetCentrality(cent);
 270   fMCAODFMD.SetNClusters(nClusters);
 271
 272   // Disable this check on SPD - will bias data
 273   // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
 274   if (found & AliFMDEventInspector::kNoFMD)     {
 275     fHStatus->Fill(4);
 276     isAccepted = false;
 277     // return;
 278   }
 279   if (found & AliFMDEventInspector::kNoVertex)  {
 280     fHStatus->Fill(5);
 281     isAccepted = false;
 282     // return;
 283   }
 284
 285   if (isAccepted) {
 286     fAODFMD.SetIpZ(ip.Z());
 287     fMCAODFMD.SetIpZ(ip.Z());
 288   }
 289   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 290
 291   // We we do not want to use low flux specific code, we disable it here.
 292   if (!fEnableLowFlux) lowFlux = false;
 293
 294
 295
 296   // Get FMD data
 297   AliESDFMD*  esdFMD  = esd.GetFMDData();
 298
 299   // Fix up the the ESD
 300   GetESDFixer().Fix(*esdFMD, ip.Z());
 301
 302   // Apply the sharing filter (or hit merging or clustering if you like)
 303   START_SW(individual);
 304   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())){
 305     AliWarning("Sharing filter failed!");
 306     fHStatus->Fill(8);
 307     return false;
 308   }
 309   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),fMCESDFMD,fPrimary)){
 310     AliWarning("MC Sharing filter failed!");
 311     fHStatus->Fill(8);
 312     return false;
 313   }
 314
 315   // Store some MC parameters in corners of histogram :-)
 316   fPrimary->SetBinContent(0,                      0,                    vzMC);
 317   fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,0,                    phiR);
 318   fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,fPrimary->GetNbinsY(),cMC);
 319
 320
 321   if (!isAccepted) {
 322     // Exit on MC event w/o trigger, vertex, data - since there's no more
 323     // to be done for MC
 324     FILL_SW(individual,kTimingSharingFilter);
 325     return false;
 326   }
 327
 328   //MarkEventForStore();
 329   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 330   FILL_SW(individual,kTimingSharingFilter);
 331
 332
 333   // Calculate the inclusive charged particle density
 334   START_SW(individual);
 335   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
 336     AliWarning("Density calculator failed!");
 337     fHStatus->Fill(9);
 338     return false;
 339   }
 340   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
 341     AliWarning("MC Density calculator failed!");
 342     fHStatus->Fill(9);
 343     return false;
 344   }
 345   fDensityCalculator.CompareResults(fHistos, fMCHistos);
 346   FILL_SW(individual,kTimingDensityCalculator);
 347
 348   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
 349     START_SW(individual);
 350     if (!fEventPlaneFinder.FindEventplane(&esd, fAODEP,
 351                                           &(fAODFMD.GetHistogram()), &fHistos)){
 352       AliWarning("Eventplane finder failed!");
 353       fHStatus->Fill(10);
 354     }
 355     FILL_SW(individual,kTimingEventPlaneFinder);
 356   }
 357
 358   // Do the secondary and other corrections.
 359   START_SW(individual);
 360   if (!fCorrections.Correct(fHistos, ivz)) {
 361     AliWarning("Corrections failed");
 362     fHStatus->Fill(12);
 363     return false;
 364   }
 365   if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
 366     AliWarning("MC Corrections failed");
 367     fHStatus->Fill(12);
 368     return false;
 369   }
 370   fCorrections.CompareResults(fHistos, fMCHistos);
 371   FILL_SW(individual,kTimingCorrections);
 372
 373   // Collect our 'super' histogram
 374   Bool_t add = fAODFMD.IsTriggerBits(AliAODForwardMult::kInel);
 375   START_SW(individual);
 376   if (!fHistCollector.Collect(fHistos,
 377                               fRingSums,
 378                               ivz,
 379                               fAODFMD.GetHistogram(),
 380                               fAODFMD.GetCentrality(),
 381                               false,
 382                               add)) {
 383     AliWarning("Histogram collector failed");
 384     fHStatus->Fill(13);
 385     return false;
 386   }
 387   if (!fHistCollector.Collect(fMCHistos,
 388                               fMCRingSums,
 389                               ivz,
 390                               fMCAODFMD.GetHistogram(),
 391                               -1,
 392                               true,
 393                               add)) {
 394     AliWarning("MC Histogram collector failed");
 395     fHStatus->Fill(13);
 396     return false;
 397   }
 398   FILL_SW(individual,kTimingHistCollector);
 399 #if 0
 400   // Copy underflow bins to overflow bins - always full phi coverage
 401   TH2&  hMC  = fMCAODFMD.GetHistogram();
 402   Int_t nEta = hMC.GetNbinsX();
 403   Int_t nY   = hMC.GetNbinsY();
 404   for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 405     hMC.SetBinContent(iEta, nY+1, hMC.GetBinContent(iEta, 0));
 406   }
 407 #endif
 408
 409   if (add) {
 410     fHData->Add(&(fAODFMD.GetHistogram()));
 411     fHStatus->Fill(15);
 412   }
 413   else {
 414     fHStatus->Fill(14);
 415   }
 416   FILL_SW(total,kTimingTotal);
 417
 418   return true;
 419 }
 420
 421 //____________________________________________________________________
 422 Bool_t
 423 AliForwardMCMultiplicityTask::PostEvent()
 424 {
 425   Bool_t ret = AliForwardMultiplicityBase::PostEvent();
 426   POST(MCAOD_SLOT, &fMCAODFMD);
 427   POST(PRIMARY_SLOT, fPrimary);
 428   return ret;
 429 }
 430
 431 //____________________________________________________________________
 432 void
 433 AliForwardMCMultiplicityTask::EstimatedNdeta(const TList* input,
 434                                              TList*       output) const
 435 {
 436   AliForwardMultiplicityBase::EstimatedNdeta(input, output);
 437   MakeRingdNdeta(input, "mcRingSums", output, "mcRingResults", 24);
 438 }
 439
 440
 441
 442 //
 443 // EOF
 444 //