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 "AliLog.h"
  19 #include "AliAODHandler.h"
  20 #include "AliInputEventHandler.h"
  21 #include "AliAnalysisManager.h"
  22 #include "AliFMDEventInspector.h"
  23 #include "AliFMDSharingFilter.h"
  24 #include "AliFMDDensityCalculator.h"
  25 #include "AliFMDCorrector.h"
  26 #include "AliFMDHistCollector.h"
  27 #include "AliFMDEventPlaneFinder.h"
  28 #include "AliESDEvent.h"
  29 #include <TROOT.h>
  30 #include <TSystem.h>
  31 #include <TAxis.h>
  32 #include <THStack.h>
  33 #include <iostream>
  34 #include <iomanip>
  35
  36 //====================================================================
  37 AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name)
  38   : AliAnalysisTaskSE(name),
  39     fEnableLowFlux(false),
  40     fFirstEvent(true),
  41     fCorrManager(0)
  42 {
  43   DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name);
  44   // Set our persistent pointer
  45   fCorrManager = &AliForwardCorrectionManager::Instance();
  46   fBranchNames =
  47     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  48     "AliESDFMD.,SPDVertex.,PrimaryVertex.";
  49 }
  50
  51 //____________________________________________________________________
  52 AliForwardMultiplicityBase&
  53 AliForwardMultiplicityBase::operator=(const AliForwardMultiplicityBase& o)
  54 {
  55   DGUARD(fDebug,2,"Assignment to AliForwardMultiplicityBase");
  56   if (&o == this) return *this;
  57   fEnableLowFlux = o.fEnableLowFlux;
  58   fFirstEvent    = o.fFirstEvent;
  59   fCorrManager   = o.fCorrManager;
  60   return *this;
  61 }
  62 //____________________________________________________________________
  63 Bool_t
  64 AliForwardMultiplicityBase::Configure(const char* macro)
  65 {
  66   // --- Configure the task ------------------------------------------
  67   TString macroPath(gROOT->GetMacroPath());
  68   if (!macroPath.Contains("$(ALICE_ROOT)/PWGLF/FORWARD/analysis2")) {
  69     macroPath.Append(":$(ALICE_ROOT)/PWGLF/FORWARD/analysis2");
  70     gROOT->SetMacroPath(macroPath);
  71   }
  72   const char* config = gSystem->Which(gROOT->GetMacroPath(), macro);
  73   if (!config) {
  74     AliWarningF("%s not found in %s", macro, gROOT->GetMacroPath());
  75     return false;
  76   }
  77
  78   AliInfoF("Loading configuration of '%s' from %s",  ClassName(), config);
  79   gROOT->Macro(Form("%s((AliForwardMultiplicityBase*)%p)", config, this));
  80   delete config;
  81
  82  return true;
  83 }
  84
  85 //____________________________________________________________________
  86 Bool_t
  87 AliForwardMultiplicityBase::CheckCorrections(UInt_t what) const
  88 {
  89   //
  90   // Check if all needed corrections are there and accounted for.  If not,
  91   // do a Fatal exit
  92   //
  93   // Parameters:
  94   //    what Which corrections is needed
  95   //
  96   // Return:
  97   //    true if all present, false otherwise
  98   //
  99   DGUARD(fDebug,1,"Checking corrections 0x%x", what);
 100
 101   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 102   // Check that we have the energy loss fits, needed by
 103   //   AliFMDSharingFilter
 104   //   AliFMDDensityCalculator
 105   if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) {
 106     AliFatal(Form("No energy loss fits"));
 107     return false;
 108   }
 109   // Check that we have the double hit correction - (optionally) used by
 110   //  AliFMDDensityCalculator
 111   if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit()) {
 112     AliFatal("No double hit corrections");
 113     return false;
 114   }
 115   // Check that we have the secondary maps, needed by
 116   //   AliFMDCorrector
 117   //   AliFMDHistCollector
 118   if (what & AliForwardCorrectionManager::kSecondaryMap &&
 119       !fcm.GetSecondaryMap()) {
 120     AliFatal("No secondary corrections");
 121     return false;
 122   }
 123   // Check that we have the vertex bias correction, needed by
 124   //   AliFMDCorrector
 125   if (what & AliForwardCorrectionManager::kVertexBias &&
 126       !fcm.GetVertexBias()) {
 127     AliFatal("No event vertex bias corrections");
 128     return false;
 129   }
 130   // Check that we have the merging efficiencies, optionally used by
 131   //   AliFMDCorrector
 132   if (what & AliForwardCorrectionManager::kMergingEfficiency &&
 133       !fcm.GetMergingEfficiency()) {
 134     AliFatal("No merging efficiencies");
 135     return false;
 136   }
 137   // Check that we have the acceptance correction, needed by
 138   //   AliFMDCorrector
 139   if (what & AliForwardCorrectionManager::kAcceptance &&
 140       !fcm.GetAcceptance()) {
 141     AliFatal("No acceptance corrections");
 142     return false;
 143   }
 144   return true;
 145 }
 146 //____________________________________________________________________
 147 Bool_t
 148 AliForwardMultiplicityBase::ReadCorrections(const TAxis*& pe,
 149                                             const TAxis*& pv,
 150                                             Bool_t        mc)
 151 {
 152   //
 153   // Read corrections
 154   //
 155   //
 156
 157   UInt_t what = AliForwardCorrectionManager::kAll;
 158   if (!fEnableLowFlux)
 159     what ^= AliForwardCorrectionManager::kDoubleHit;
 160   if (!GetCorrections().IsUseVertexBias())
 161     what ^= AliForwardCorrectionManager::kVertexBias;
 162   if (!GetCorrections().IsUseAcceptance())
 163     what ^= AliForwardCorrectionManager::kAcceptance;
 164   if (!GetCorrections().IsUseMergingEfficiency())
 165     what ^= AliForwardCorrectionManager::kMergingEfficiency;
 166   DGUARD(fDebug,1,"Read corrections 0x%x", what);
 167
 168   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 169   if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
 170                 GetEventInspector().GetEnergy(),
 171                 GetEventInspector().GetField(),
 172                 mc,
 173                 what)) {
 174     AliWarning("Failed to read in some corrections, making task zombie");
 175     return false;
 176   }
 177   if (!CheckCorrections(what)) return false;
 178
 179   // Sett our persistency pointer
 180   // fCorrManager = &fcm;
 181
 182   // Get the eta axis from the secondary maps - if read in
 183   if (!pe) {
 184     pe = fcm.GetEtaAxis();
 185     if (!pe) AliFatal("No eta axis defined");
 186   }
 187   // Get the vertex axis from the secondary maps - if read in
 188   if (!pv) {
 189     pv = fcm.GetVertexAxis();
 190     if (!pv) AliFatal("No vertex axis defined");
 191   }
 192
 193   return true;
 194 }
 195 //____________________________________________________________________
 196 AliESDEvent*
 197 AliForwardMultiplicityBase::GetESDEvent()
 198 {
 199   //
 200   // Get the ESD event. IF this is the first event, initialise
 201   //
 202   DGUARD(fDebug,1,"Get the ESD event");
 203   if (IsZombie()) return 0;
 204   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 205   if (!esd) {
 206     AliWarning("No ESD event found for input event");
 207     return 0;
 208   }
 209
 210   // --- Load the data -----------------------------------------------
 211   LoadBranches();
 212
 213   // On the first event, initialize the parameters
 214   if (fFirstEvent && esd->GetESDRun()) {
 215     GetEventInspector().ReadRunDetails(esd);
 216
 217     AliInfo(Form("Initializing with parameters from the ESD:\n"
 218                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 219                  "         AliESDEvent::GetBeamType()     ->%s\n"
 220                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 221                  "         AliESDEvent::GetMagneticField()->%f\n"
 222                  "         AliESDEvent::GetRunNumber()    ->%d",
 223                  esd->GetBeamEnergy(),
 224                  esd->GetBeamType(),
 225                  esd->GetCurrentL3(),
 226                  esd->GetMagneticField(),
 227                  esd->GetRunNumber()));
 228
 229     fFirstEvent = false;
 230
 231     GetEventPlaneFinder().SetRunNumber(esd->GetRunNumber());
 232     if (!SetupForData()) {
 233       AliError("Failed to initialize sub-algorithms, making this a zombie");
 234       esd = 0; // Make sure we do nothing on this event
 235       Info("GetESDEvent", "ESD event pointer %p", esd);
 236       SetZombie(true);
 237       // return 0;
 238     }
 239   }
 240   return esd;
 241 }
 242 //____________________________________________________________________
 243 void
 244 AliForwardMultiplicityBase::MarkEventForStore() const
 245 {
 246   // Make sure the AOD tree is filled
 247   DGUARD(fDebug,3,"Mark AOD event for storage");
 248   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 249   AliAODHandler*      ah =
 250     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 251   if (!ah)
 252     AliFatal("No AOD output handler set in analysis manager");
 253
 254   ah->SetFillAOD(kTRUE);
 255 }
 256
 257 //____________________________________________________________________
 258 Bool_t
 259 AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input,
 260                                              TList*       output,
 261                                              Double_t&    nTr,
 262                                              Double_t&    nTrVtx,
 263                                              Double_t&    nAcc)
 264 {
 265   // Get our histograms from the container
 266   TH1I* hEventsTr    = 0;
 267   TH1I* hEventsTrVtx = 0;
 268   TH1I* hEventsAcc   = 0;
 269   TH1I* hTriggers    = 0;
 270   if (!GetEventInspector().FetchHistograms(input,
 271                                            hEventsTr,
 272                                            hEventsTrVtx,
 273                                            hEventsAcc,
 274                                            hTriggers)) {
 275     AliError(Form("Didn't get histograms from event selector "
 276                   "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
 277                   hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
 278     input->ls();
 279     return false;
 280   }
 281   nTr             = hEventsTr->Integral();
 282   nTrVtx          = hEventsTrVtx->Integral();
 283   nAcc            = hEventsAcc->Integral();
 284   Double_t vtxEff = nTrVtx / nTr;
 285   TH2D*   hData   = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
 286   if (!hData) {
 287     AliError(Form("Couldn't get our summed histogram from output "
 288                   "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
 289     input->ls();
 290     return false;
 291   }
 292
 293   Int_t nY      = hData->GetNbinsY();
 294   TH1D* dNdeta  = hData->ProjectionX("dNdeta",  1,     nY, "e");
 295   TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1,     nY, "e");
 296   TH1D* norm    = hData->ProjectionX("norm",    0,     0,  "");
 297   TH1D* phi     = hData->ProjectionX("phi",     nY+1,  nY+1,  "");
 298   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 299   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 300   dNdeta->SetMarkerColor(kRed+1);
 301   dNdeta->SetMarkerStyle(20);
 302   dNdeta->SetDirectory(0);
 303
 304   dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
 305   dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 306   dNdeta_->SetMarkerColor(kMagenta+1);
 307   dNdeta_->SetMarkerStyle(21);
 308   dNdeta_->SetDirectory(0);
 309
 310   norm->SetTitle("Normalization to #eta coverage");
 311   norm->SetYTitle("#eta coverage");
 312   norm->SetMarkerColor(kBlue+1);
 313   norm->SetMarkerStyle(21);
 314   norm->SetFillColor(kBlue+1);
 315   norm->SetFillStyle(3005);
 316   norm->SetDirectory(0);
 317
 318   phi->SetTitle("Normalization to #phi acceptance");
 319   phi->SetYTitle("#phi acceptance");
 320   phi->SetMarkerColor(kGreen+1);
 321   phi->SetMarkerStyle(20);
 322   phi->SetFillColor(kGreen+1);
 323   phi->SetFillStyle(3004);
 324   // phi->Scale(1. / nAcc);
 325   phi->SetDirectory(0);
 326
 327   // dNdeta->Divide(norm);
 328   dNdeta->Divide(phi);
 329   dNdeta->SetStats(0);
 330   dNdeta->Scale(vtxEff, "width");
 331
 332   dNdeta_->Divide(norm);
 333   dNdeta_->SetStats(0);
 334   dNdeta_->Scale(vtxEff, "width");
 335
 336   output->Add(dNdeta);
 337   output->Add(dNdeta_);
 338   output->Add(norm);
 339   output->Add(phi);
 340
 341   return true;
 342 }
 343
 344
 345 //____________________________________________________________________
 346 void
 347 AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
 348                                            const char*  inName,
 349                                            TList*       output,
 350                                            const char*  outName,
 351                                            Int_t        style) const
 352 {
 353   // Make dN/deta for each ring found in the input list.
 354   //
 355   // A stack of all the dN/deta is also made for easy drawing.
 356   //
 357   // Note, that the distributions are normalised to the number of
 358   // observed events only - they should be corrected for
 359   DGUARD(fDebug,3,"Make first-shot ring dN/deta");
 360
 361   if (!input) return;
 362   TList* list = static_cast<TList*>(input->FindObject(inName));
 363   if (!list) {
 364     AliWarning(Form("No list %s found in %s", inName, input->GetName()));
 365     return;
 366   }
 367
 368   TList* out = new TList;
 369   out->SetName(outName);
 370   out->SetOwner();
 371   output->Add(out);
 372
 373   THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
 374   const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 375   const char** ptr         = names;
 376
 377   while (*ptr) {
 378     TList* thisList = new TList;
 379     thisList->SetOwner();
 380     thisList->SetName(*ptr);
 381     out->Add(thisList);
 382
 383     TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
 384     if (!h) {
 385       AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
 386       ptr++;
 387       continue;
 388     }
 389     TH2D* sumPhi = static_cast<TH2D*>(h->Clone("sum_phi"));
 390     sumPhi->SetDirectory(0);
 391     thisList->Add(sumPhi);
 392
 393     TH2D* sumEta = static_cast<TH2D*>(h->Clone("sum_eta"));
 394     sumEta->SetDirectory(0);
 395     thisList->Add(sumEta);
 396
 397     Int_t nY   = sumEta->GetNbinsY();
 398     TH1D* etaCov =static_cast<TH1D*>(h->ProjectionX("etaCov", 0,    0,    ""));
 399     TH1D* phiAcc =static_cast<TH1D*>(h->ProjectionX("phiAcc", nY+1, nY+1, ""));
 400
 401     etaCov->SetTitle("Normalization to #eta coverage");
 402     etaCov->SetYTitle("#eta coverage");
 403     etaCov->SetMarkerColor(kBlue+1);
 404     etaCov->SetFillColor(kBlue+1);
 405     etaCov->SetFillStyle(3005);
 406     etaCov->SetDirectory(0);
 407
 408     phiAcc->SetTitle("Normalization to #phi acceptance");
 409     phiAcc->SetYTitle("#phi acceptance");
 410     phiAcc->SetMarkerColor(kGreen+1);
 411     phiAcc->SetFillColor(kGreen+1);
 412     phiAcc->SetFillStyle(3004);
 413     // phiAcc->Scale(1. / nAcc);
 414     phiAcc->SetDirectory(0);
 415
 416     // Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1);
 417     for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) {
 418       for (Int_t j = 1; j <= nY; j++) {
 419         Double_t c = sumEta->GetBinContent(i, j);
 420         Double_t e = sumEta->GetBinError(i, j);
 421         Double_t a = etaCov->GetBinContent(i);
 422         Double_t p = phiAcc->GetBinContent(i);
 423         // Double_t t = p; // * a
 424         sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a);
 425         sumEta->SetBinError(  i, j, a <= 0 ? 0 : e / a);
 426         sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p);
 427         sumPhi->SetBinError(  i, j, p <= 0 ? 0 : e / p);
 428       }
 429     }
 430     // etaCov->Scale(s);
 431     // phiAcc->Scale(s);
 432
 433     TH1D* resPhi  =static_cast<TH1D*>(sumPhi->ProjectionX("dndeta_phi",
 434                                                           1,nY,"e"));
 435     resPhi->SetMarkerStyle(style);
 436     resPhi->SetDirectory(0);
 437     resPhi->Scale(1, "width");
 438
 439     TH1D* resEta  =static_cast<TH1D*>(sumEta->ProjectionX("dndeta_eta",
 440                                                           1,nY,"e"));
 441     resEta->SetMarkerStyle(style);
 442     resEta->SetDirectory(0);
 443     resEta->Scale(1, "width");
 444
 445     thisList->Add(resEta);
 446     thisList->Add(etaCov);
 447     thisList->Add(resPhi);
 448     thisList->Add(phiAcc);
 449     dndetaRings->Add(resPhi);
 450     ptr++;
 451   }
 452   out->Add(dndetaRings);
 453 }
 454 //____________________________________________________________________
 455 void
 456 AliForwardMultiplicityBase::Print(Option_t* option) const
 457 {
 458   //
 459   // Print information
 460   //
 461   // Parameters:
 462   //    option Not used
 463   //
 464
 465   std::cout << ClassName() << ": " << GetName() << "\n"
 466             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 467             << "\n"
 468             << "  Off-line trigger mask:  0x"
 469             << std::hex     << std::setfill('0')
 470             << std::setw (8) << fOfflineTriggerMask
 471             << std::dec     << std::setfill (' ') << std::endl;
 472   gROOT->IncreaseDirLevel();
 473   if (fCorrManager) fCorrManager->Print();
 474   else
 475     std::cout << "  Correction manager not set yet" << std::endl;
 476   GetEventInspector()   .Print(option);
 477   GetSharingFilter()    .Print(option);
 478   GetDensityCalculator().Print(option);
 479   GetCorrections()      .Print(option);
 480   GetHistCollector()    .Print(option);
 481   GetEventPlaneFinder() .Print(option);
 482   gROOT->DecreaseDirLevel();
 483 }
 484
 485
 486 //
 487 // EOF
 488 //