PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx

   1 //====================================================================
   2 //
   3 // Base class for classes that calculate the multiplicity in the
   4 // forward regions event-by-event
   5 //
   6 // Inputs:
   7 //   - AliESDEvent
   8 //
   9 // Outputs:
  10 //   - AliAODForwardMult
  11 //
  12 // Histograms
  13 //
  14 // Corrections used
  15 #include "AliForwardMultiplicityBase.h"
  16 #include "AliForwardCorrectionManager.h"
  17 #include "AliForwardUtil.h"
  18 #include "AliFMDCorrELossFit.h"
  19 #include "AliLog.h"
  20 #include "AliAODHandler.h"
  21 #include "AliInputEventHandler.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliFMDEventInspector.h"
  24 #include "AliFMDSharingFilter.h"
  25 #include "AliFMDDensityCalculator.h"
  26 #include "AliFMDCorrector.h"
  27 #include "AliFMDHistCollector.h"
  28 #include "AliFMDEventPlaneFinder.h"
  29 #include "AliESDEvent.h"
  30 #include <TROOT.h>
  31 #include <TSystem.h>
  32 #include <TAxis.h>
  33 #include <TProfile.h>
  34 #include <THStack.h>
  35 #include <iostream>
  36 #include <iomanip>
  37
  38 //====================================================================
  39 AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name)
  40   : AliBaseESDTask(name, "AliForwardMultiplicityBase",
  41                    &(AliForwardCorrectionManager::Instance())),
  42     fEnableLowFlux(false),
  43     fStorePerRing(false),
  44     fHData(0),
  45     fHistos(),
  46     fAODFMD(false),
  47     fAODEP(false),
  48     fRingSums(),
  49     fDoTiming(false),
  50     fHTiming(0)
  51 {
  52   DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name);
  53 }
  54
  55
  56 //____________________________________________________________________
  57 void
  58 AliForwardMultiplicityBase::SetDebug(Int_t dbg)
  59 {
  60   //
  61   // Set debug level
  62   //
  63   // Parameters:
  64   //    dbg debug level
  65   //
  66   AliBaseESDTask::        SetDebug(dbg);
  67   GetSharingFilter()    .SetDebug(dbg);
  68   GetDensityCalculator().SetDebug(dbg);
  69   GetCorrections()      .SetDebug(dbg);
  70   GetHistCollector()    .SetDebug(dbg);
  71   GetEventPlaneFinder() .SetDebug(dbg);
  72 }
  73
  74 //____________________________________________________________________
  75 Bool_t
  76 AliForwardMultiplicityBase::Book()
  77 {
  78   //
  79   // Create output objects
  80   //
  81   //
  82   DGUARD(fDebug,1,"Create user ouput");
  83   UInt_t what = AliForwardCorrectionManager::kAll;
  84   if (!fEnableLowFlux)
  85     what ^= AliForwardCorrectionManager::kDoubleHit;
  86   if (!GetCorrections().IsUseVertexBias())
  87     what ^= AliForwardCorrectionManager::kVertexBias;
  88   if (!GetCorrections().IsUseAcceptance())
  89     what ^= AliForwardCorrectionManager::kAcceptance;
  90   if (!GetCorrections().IsUseMergingEfficiency())
  91     what ^= AliForwardCorrectionManager::kMergingEfficiency;
  92   fNeededCorrections = what;
  93
  94   GetSharingFilter()    .CreateOutputObjects(fList);
  95   GetDensityCalculator().CreateOutputObjects(fList);
  96   GetCorrections()      .CreateOutputObjects(fList);
  97   GetHistCollector()    .CreateOutputObjects(fList);
  98   GetEventPlaneFinder() .CreateOutputObjects(fList);
  99
 100   TAxis tmp(1, 0, 1);
 101   fHistos.Init(tmp);
 102
 103   if (fDebug > 1) fDoTiming = true;
 104   if (fDoTiming) {
 105     fHTiming = new TProfile("timing", "Timing of task",
 106                             kTimingTotal, 0.5, kTimingTotal+.5);
 107     fHTiming->SetDirectory(0);
 108     fHTiming->SetFillColor(kRed+1);
 109     fHTiming->SetFillStyle(3001);
 110     fHTiming->SetMarkerStyle(20);
 111     fHTiming->SetMarkerColor(kBlack);
 112     fHTiming->SetLineColor(kBlack);
 113     fHTiming->SetXTitle("Part");
 114     fHTiming->SetYTitle("#LTt_{part}#GT [CPU]");
 115     fHTiming->SetStats(0);
 116     TAxis* xaxis = fHTiming->GetXaxis();
 117     xaxis->SetBinLabel(kTimingEventInspector,
 118                        GetEventInspector()   .GetName());
 119     xaxis->SetBinLabel(kTimingSharingFilter,
 120                        GetSharingFilter()    .GetName());
 121     xaxis->SetBinLabel(kTimingDensityCalculator,
 122                        GetDensityCalculator().GetName());
 123     xaxis->SetBinLabel(kTimingCorrections,
 124                        GetCorrections()      .GetName());
 125     xaxis->SetBinLabel(kTimingHistCollector,
 126                        GetHistCollector()    .GetName());
 127     xaxis->SetBinLabel(kTimingEventPlaneFinder,
 128                        GetEventPlaneFinder() .GetName());
 129     xaxis->SetBinLabel(kTimingTotal, "Total");
 130     fList->Add(fHTiming);
 131   }
 132   return true;
 133 }
 134 //____________________________________________________________________
 135 void
 136 AliForwardMultiplicityBase::CreateBranches(AliAODHandler* ah)
 137 {
 138   TObject* obj   = &fAODFMD; ah->AddBranch("AliAODForwardMult", &obj);
 139   TObject* epobj = &fAODEP;  ah->AddBranch("AliAODForwardEP", &epobj);
 140
 141   if (!fStorePerRing) return;
 142
 143   AliWarning("Per-ring histograms in AOD\n"
 144              "*********************************************************\n"
 145              "* For each event 5 additional 2D histogram are stored   *\n"
 146              "* in separate branches of the AODs.  This will increase *\n"
 147              "* the size of the AODs - proceed with caution           *\n"
 148              "*********************************************************");
 149   TObject* hists[] = { fHistos.fFMD1i,
 150                        fHistos.fFMD2i, fHistos.fFMD2o,
 151                        fHistos.fFMD3i, fHistos.fFMD3o };
 152   for (Int_t i = 0; i < 5; i++) {
 153     ah->AddBranch("TH2D", &(hists[i]));
 154   }
 155 }
 156
 157 //____________________________________________________________________
 158 Bool_t
 159 AliForwardMultiplicityBase::PreData(const TAxis& vertex, const TAxis& eta)
 160 {
 161   //
 162   // Initialise the sub objects and stuff.  Called on first event
 163   //
 164   //
 165   DGUARD(fDebug,1,"Initialize sub-algorithms");
 166
 167   // Force this here so we select the proper quality
 168   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 169
 170   UInt_t what = fNeededCorrections;
 171   // Check that we have the energy loss fits, needed by
 172   //   AliFMDSharingFilter
 173   //   AliFMDDensityCalculator
 174   if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit())
 175     AliFatal(Form("No energy loss fits"));
 176
 177   // Check that we have the double hit correction - (optionally) used by
 178   //  AliFMDDensityCalculator
 179   if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit())
 180     AliFatal("No double hit corrections");
 181
 182   // Check that we have the secondary maps, needed by
 183   //   AliFMDCorrector
 184   //   AliFMDHistCollector
 185   if (what & AliForwardCorrectionManager::kSecondaryMap &&
 186       !fcm.GetSecondaryMap())
 187     AliFatal("No secondary corrections");
 188
 189   // Check that we have the vertex bias correction, needed by
 190   //   AliFMDCorrector
 191   if (what & AliForwardCorrectionManager::kVertexBias && !fcm.GetVertexBias())
 192     AliFatal("No event vertex bias corrections");
 193
 194   // Check that we have the merging efficiencies, optionally used by
 195   //   AliFMDCorrector
 196   if (what & AliForwardCorrectionManager::kMergingEfficiency &&
 197       !fcm.GetMergingEfficiency())
 198     AliFatal("No merging efficiencies");
 199
 200   // Check that we have the acceptance correction, needed by
 201   //   AliFMDCorrector
 202   if (what & AliForwardCorrectionManager::kAcceptance && !fcm.GetAcceptance())
 203     AliFatal("No acceptance corrections");
 204
 205   const AliFMDCorrELossFit* fits = fcm.GetELossFit();
 206   fits->CacheBins(GetDensityCalculator().GetMinQuality());
 207
 208   InitMembers(eta,vertex);
 209
 210   GetSharingFilter()    .SetupForData(eta);
 211   GetDensityCalculator().SetupForData(eta);
 212   GetCorrections()      .SetupForData(eta);
 213   GetHistCollector()    .SetupForData(vertex,eta);
 214
 215   GetEventPlaneFinder() .SetRunNumber(GetEventInspector().GetRunNumber());
 216   GetEventPlaneFinder() .SetupForData(eta);
 217
 218   fAODFMD.SetBit(AliAODForwardMult::kSecondary,
 219                  GetCorrections().IsUseSecondaryMap());
 220   fAODFMD.SetBit(AliAODForwardMult::kVertexBias,
 221                  GetCorrections().IsUseVertexBias());
 222   fAODFMD.SetBit(AliAODForwardMult::kAcceptance,
 223                  GetCorrections().IsUseAcceptance());
 224   fAODFMD.SetBit(AliAODForwardMult::kMergingEfficiency,
 225                  GetCorrections().IsUseMergingEfficiency());
 226   fAODFMD.SetBit(AliAODForwardMult::kSum,
 227                  GetHistCollector().GetMergeMethod() ==
 228                  AliFMDHistCollector::kSum);
 229   return true;
 230 }
 231
 232 //____________________________________________________________________
 233 void
 234 AliForwardMultiplicityBase::InitMembers(const TAxis& eta, const TAxis& /*pv*/)
 235 {
 236   fHistos.ReInit(eta);
 237   fAODFMD.Init(eta);
 238   fAODEP.Init(eta);
 239   fRingSums.Init(eta);
 240
 241   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 242   fHData->SetStats(0);
 243   fHData->SetDirectory(0);
 244   fList->Add(fHData);
 245
 246   TList* rings = new TList;
 247   rings->SetName("ringSums");
 248   rings->SetOwner();
 249   fList->Add(rings);
 250
 251   rings->Add(fRingSums.Get(1, 'I'));
 252   rings->Add(fRingSums.Get(2, 'I'));
 253   rings->Add(fRingSums.Get(2, 'O'));
 254   rings->Add(fRingSums.Get(3, 'I'));
 255   rings->Add(fRingSums.Get(3, 'O'));
 256   fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
 257   fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
 258   fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
 259   fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
 260   fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
 261 }
 262 //____________________________________________________________________
 263 Bool_t
 264 AliForwardMultiplicityBase::Finalize()
 265 {
 266   //
 267   // End of job
 268   //
 269   // Parameters:
 270   //    option Not used
 271   //
 272   DGUARD(fDebug,1,"Processing the merged results");
 273
 274   TList* list   = fList;
 275   TList* output = fResults;
 276   Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
 277   MakeSimpledNdeta(list, output, nTr, nTrVtx, nAcc);
 278
 279   EstimatedNdeta(list, output);
 280
 281   GetSharingFilter()    .Terminate(list,output,Int_t(nTr));
 282   GetDensityCalculator().Terminate(list,output,Int_t(nTrVtx));
 283   GetCorrections()      .Terminate(list,output,Int_t(nTrVtx));
 284
 285   TProfile* timing = static_cast<TProfile*>(list->FindObject("timing"));
 286   if (timing) {
 287     TProfile* p = static_cast<TProfile*>(timing->Clone());
 288     p->SetDirectory(0);
 289     p->Scale(100. / p->GetBinContent(p->GetNbinsX()));
 290     p->SetYTitle("#LTt_{part}#GT/#LTt_{total}#GT [%]");
 291     p->SetTitle("Relative timing of task");
 292     output->Add(p);
 293   }
 294   return true;
 295 }
 296
 297 //____________________________________________________________________
 298 void
 299 AliForwardMultiplicityBase::EstimatedNdeta(const TList* input,
 300                                            TList*       output) const
 301 {
 302   MakeRingdNdeta(input, "ringSums", output, "ringResults");
 303 }
 304
 305 //____________________________________________________________________
 306 Bool_t
 307 AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input,
 308                                              TList*       output,
 309                                              Double_t&    nTr,
 310                                              Double_t&    nTrVtx,
 311                                              Double_t&    nAcc)
 312 {
 313   // Get our histograms from the container
 314   TH1I* hEventsTr    = 0;
 315   TH1I* hEventsTrVtx = 0;
 316   TH1I* hEventsAcc   = 0;
 317   TH1I* hTriggers    = 0;
 318   if (!GetEventInspector().FetchHistograms(input,
 319                                            hEventsTr,
 320                                            hEventsTrVtx,
 321                                            hEventsAcc,
 322                                            hTriggers)) {
 323     AliError(Form("Didn't get histograms from event selector "
 324                   "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
 325                   hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
 326     input->ls();
 327     return false;
 328   }
 329   nTr             = hEventsTr->Integral();
 330   nTrVtx          = hEventsTrVtx->Integral();
 331   nAcc            = hEventsAcc->Integral();
 332   Double_t vtxEff = nTrVtx / nTr;
 333   TH2D*   hData   = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
 334   if (!hData) {
 335     AliError(Form("Couldn't get our summed histogram from output "
 336                   "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
 337     input->ls();
 338     return false;
 339   }
 340
 341   Int_t nY      = hData->GetNbinsY();
 342   TH1D* dNdeta  = hData->ProjectionX("dNdeta",  1,     nY, "e");
 343   TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1,     nY, "e");
 344   TH1D* norm    = hData->ProjectionX("norm",    0,     0,  "");
 345   TH1D* phi     = hData->ProjectionX("phi",     nY+1,  nY+1,  "");
 346   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 347   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 348   dNdeta->SetMarkerColor(kRed+1);
 349   dNdeta->SetMarkerStyle(20);
 350   dNdeta->SetDirectory(0);
 351
 352   dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
 353   dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 354   dNdeta_->SetMarkerColor(kMagenta+1);
 355   dNdeta_->SetMarkerStyle(21);
 356   dNdeta_->SetDirectory(0);
 357
 358   norm->SetTitle("Normalization to  #eta coverage");
 359   norm->SetYTitle("#eta coverage");
 360   norm->SetLineColor(kBlue+1);
 361   norm->SetMarkerColor(kBlue+1);
 362   norm->SetMarkerStyle(21);
 363   norm->SetFillColor(kBlue+1);
 364   norm->SetFillStyle(3005);
 365   norm->SetDirectory(0);
 366
 367   phi->SetTitle("Normalization to  #phi acceptance");
 368   phi->SetYTitle("#phi acceptance");
 369   phi->SetLineColor(kGreen+1);
 370   phi->SetMarkerColor(kGreen+1);
 371   phi->SetMarkerStyle(20);
 372   phi->SetFillColor(kGreen+1);
 373   phi->SetFillStyle(3004);
 374   // phi->Scale(1. / nAcc);
 375   phi->SetDirectory(0);
 376
 377   // dNdeta->Divide(norm);
 378   dNdeta->Divide(phi);
 379   dNdeta->SetStats(0);
 380   dNdeta->Scale(vtxEff, "width");
 381
 382   dNdeta_->Divide(norm);
 383   dNdeta_->SetStats(0);
 384   dNdeta_->Scale(vtxEff, "width");
 385
 386   output->Add(dNdeta);
 387   output->Add(dNdeta_);
 388   output->Add(norm);
 389   output->Add(phi);
 390
 391   return true;
 392 }
 393
 394
 395 //____________________________________________________________________
 396 void
 397 AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
 398                                            const char*  inName,
 399                                            TList*       output,
 400                                            const char*  outName,
 401                                            Int_t        style) const
 402 {
 403   // Make dN/deta for each ring found in the input list.
 404   //
 405   // A stack of all the dN/deta is also made for easy drawing.
 406   //
 407   // Note, that the distributions are normalised to the number of
 408   // observed events only - they should be corrected for
 409   DGUARD(fDebug,3,"Make first-shot ring dN/deta");
 410
 411   if (!input) return;
 412   TList* list = static_cast<TList*>(input->FindObject(inName));
 413   if (!list) {
 414     AliWarning(Form("No list %s found in %s", inName, input->GetName()));
 415     return;
 416   }
 417
 418   TList* out = new TList;
 419   out->SetName(outName);
 420   out->SetOwner();
 421   output->Add(out);
 422
 423   THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
 424   const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 425   const char** ptr         = names;
 426
 427   while (*ptr) {
 428     TList* thisList = new TList;
 429     thisList->SetOwner();
 430     thisList->SetName(*ptr);
 431     out->Add(thisList);
 432
 433     TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
 434     if (!h) {
 435       AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
 436       ptr++;
 437       continue;
 438     }
 439     TH2D* sumPhi = static_cast<TH2D*>(h->Clone("sum_phi"));
 440     sumPhi->SetDirectory(0);
 441     thisList->Add(sumPhi);
 442
 443     TH2D* sumEta = static_cast<TH2D*>(h->Clone("sum_eta"));
 444     sumEta->SetDirectory(0);
 445     thisList->Add(sumEta);
 446
 447     Int_t nY   = sumEta->GetNbinsY();
 448     TH1D* etaCov =static_cast<TH1D*>(h->ProjectionX("etaCov", 0,    0,    ""));
 449     TH1D* phiAcc =static_cast<TH1D*>(h->ProjectionX("phiAcc", nY+1, nY+1, ""));
 450
 451     etaCov->SetTitle("Normalization to #eta coverage");
 452     etaCov->SetYTitle("#eta coverage");
 453     etaCov->SetMarkerColor(kBlue+1);
 454     etaCov->SetFillColor(kBlue+1);
 455     etaCov->SetFillStyle(3005);
 456     etaCov->SetDirectory(0);
 457
 458     phiAcc->SetTitle("Normalization to #phi acceptance");
 459     phiAcc->SetYTitle("#phi acceptance");
 460     phiAcc->SetMarkerColor(kGreen+1);
 461     phiAcc->SetFillColor(kGreen+1);
 462     phiAcc->SetFillStyle(3004);
 463     // phiAcc->Scale(1. / nAcc);
 464     phiAcc->SetDirectory(0);
 465
 466     // Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1);
 467     for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) {
 468       for (Int_t j = 1; j <= nY; j++) {
 469         Double_t c = sumEta->GetBinContent(i, j);
 470         Double_t e = sumEta->GetBinError(i, j);
 471         Double_t a = etaCov->GetBinContent(i);
 472         Double_t p = phiAcc->GetBinContent(i);
 473         // Double_t t = p; // * a
 474         sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a);
 475         sumEta->SetBinError(  i, j, a <= 0 ? 0 : e / a);
 476         sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p);
 477         sumPhi->SetBinError(  i, j, p <= 0 ? 0 : e / p);
 478       }
 479     }
 480     // etaCov->Scale(s);
 481     // phiAcc->Scale(s);
 482
 483     TH1D* resPhi  =static_cast<TH1D*>(sumPhi->ProjectionX("dndeta_phi",
 484                                                           1,nY,"e"));
 485     resPhi->SetMarkerStyle(style);
 486     resPhi->SetDirectory(0);
 487     resPhi->Scale(1, "width");
 488
 489     TH1D* resEta  =static_cast<TH1D*>(sumEta->ProjectionX("dndeta_eta",
 490                                                           1,nY,"e"));
 491     resEta->SetMarkerStyle(style);
 492     resEta->SetDirectory(0);
 493     resEta->Scale(1, "width");
 494
 495     thisList->Add(resEta);
 496     thisList->Add(etaCov);
 497     thisList->Add(resPhi);
 498     thisList->Add(phiAcc);
 499     dndetaRings->Add(resPhi);
 500     ptr++;
 501   }
 502   out->Add(dndetaRings);
 503 }
 504 #define PFB(N,FLAG)                             \
 505   do {                                                                  \
 506     AliForwardUtil::PrintName(N);                                       \
 507     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 508   } while(false)
 509 //____________________________________________________________________
 510 void
 511 AliForwardMultiplicityBase::Print(Option_t* option) const
 512 {
 513   //
 514   // Print information
 515   //
 516   // Parameters:
 517   //    option Not used
 518   //
 519   AliBaseESDTask::Print(option);
 520   gROOT->IncreaseDirLevel();
 521   PFB("Enable low flux code", fEnableLowFlux);
 522   PFB("Store per-ring hists", fStorePerRing);
 523   PFB("Make timing histogram", fDoTiming);
 524   // gROOT->IncreaseDirLevel();
 525   GetSharingFilter()    .Print(option);
 526   GetDensityCalculator().Print(option);
 527   GetCorrections()      .Print(option);
 528   GetHistCollector()    .Print(option);
 529   GetEventPlaneFinder() .Print(option);
 530   // gROOT->DecreaseDirLevel();
 531   gROOT->DecreaseDirLevel();
 532 }
 533
 534
 535 //
 536 // EOF
 537 //