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(0,0,"Named Construction 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   // On the first event, initialize the parameters
 211   if (fFirstEvent && esd->GetESDRun()) {
 212     GetEventInspector().ReadRunDetails(esd);
 213
 214     AliInfo(Form("Initializing with parameters from the ESD:\n"
 215                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 216                  "         AliESDEvent::GetBeamType()     ->%s\n"
 217                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 218                  "         AliESDEvent::GetMagneticField()->%f\n"
 219                  "         AliESDEvent::GetRunNumber()    ->%d",
 220                  esd->GetBeamEnergy(),
 221                  esd->GetBeamType(),
 222                  esd->GetCurrentL3(),
 223                  esd->GetMagneticField(),
 224                  esd->GetRunNumber()));
 225
 226     fFirstEvent = false;
 227
 228     GetEventPlaneFinder().SetRunNumber(esd->GetRunNumber());
 229     if (!InitializeSubs()) {
 230       AliError("Failed to initialize sub-algorithms, making this a zombie");
 231       esd = 0; // Make sure we do nothing on this event
 232       Info("GetESDEvent", "ESD event pointer %p", esd);
 233       SetZombie(true);
 234       // return 0;
 235     }
 236   }
 237   return esd;
 238 }
 239 //____________________________________________________________________
 240 void
 241 AliForwardMultiplicityBase::MarkEventForStore() const
 242 {
 243   // Make sure the AOD tree is filled
 244   DGUARD(fDebug,3,"Mark AOD event for storage");
 245   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 246   AliAODHandler*      ah =
 247     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 248   if (!ah)
 249     AliFatal("No AOD output handler set in analysis manager");
 250
 251   ah->SetFillAOD(kTRUE);
 252 }
 253
 254 //____________________________________________________________________
 255 void
 256 AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
 257                                            const char*  inName,
 258                                            TList*       output,
 259                                            const char*  outName,
 260                                            Int_t        style) const
 261 {
 262   // Make dN/deta for each ring found in the input list.
 263   //
 264   // A stack of all the dN/deta is also made for easy drawing.
 265   //
 266   // Note, that the distributions are normalised to the number of
 267   // observed events only - they should be corrected for
 268   DGUARD(fDebug,3,"Make first-shot ring dN/deta");
 269
 270   if (!input) return;
 271   TList* list = static_cast<TList*>(input->FindObject(inName));
 272   if (!list) {
 273     AliWarning(Form("No list %s found in %s", inName, input->GetName()));
 274     return;
 275   }
 276
 277   TList* out = new TList;
 278   out->SetName(outName);
 279   out->SetOwner();
 280   output->Add(out);
 281
 282   THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
 283   const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 284   const char** ptr         = names;
 285
 286   while (*ptr) {
 287     TList* thisList = new TList;
 288     thisList->SetOwner();
 289     thisList->SetName(*ptr);
 290     out->Add(thisList);
 291
 292     TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
 293     if (!h) {
 294       AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
 295       ptr++;
 296       continue;
 297     }
 298     TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
 299     copy->SetDirectory(0);
 300     thisList->Add(copy);
 301
 302     TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
 303     for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 304       for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
 305         Double_t c = copy->GetBinContent(i, j);
 306         Double_t e = copy->GetBinError(i, j);
 307         Double_t a = norm->GetBinContent(i);
 308         copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
 309         copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
 310       }
 311     }
 312
 313     TH1D* res  =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
 314                                                      h->GetNbinsY(),"e"));
 315     TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
 316     res->SetTitle(*ptr);
 317     res->Scale(1., "width");
 318     copy->Scale(1., "width");
 319
 320     if(norm->GetMaximum() > 0) {
 321       proj->Scale(1. / norm->GetMaximum(), "width");
 322       norm->Scale(1. / norm->GetMaximum());
 323     }
 324
 325     res->SetMarkerStyle(style);
 326     norm->SetDirectory(0);
 327     res->SetDirectory(0);
 328     proj->SetDirectory(0);
 329     thisList->Add(norm);
 330     thisList->Add(res);
 331     thisList->Add(proj);
 332     dndetaRings->Add(res);
 333     ptr++;
 334   }
 335   out->Add(dndetaRings);
 336 }
 337
 338 //____________________________________________________________________
 339 void
 340 AliForwardMultiplicityBase::Print(Option_t* option) const
 341 {
 342   //
 343   // Print information
 344   //
 345   // Parameters:
 346   //    option Not used
 347   //
 348
 349   std::cout << ClassName() << ": " << GetName() << "\n"
 350             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 351             << "\n"
 352             << "  Off-line trigger mask:  0x"
 353             << std::hex     << std::setfill('0')
 354             << std::setw (8) << fOfflineTriggerMask
 355             << std::dec     << std::setfill (' ') << std::endl;
 356   gROOT->IncreaseDirLevel();
 357   if (fCorrManager) fCorrManager->Print();
 358   else
 359     std::cout << "  Correction manager not set yet" << std::endl;
 360   GetEventInspector()   .Print(option);
 361   GetSharingFilter()    .Print(option);
 362   GetDensityCalculator().Print(option);
 363   GetCorrections()      .Print(option);
 364   GetHistCollector()    .Print(option);
 365   GetEventPlaneFinder() .Print(option);
 366   gROOT->DecreaseDirLevel();
 367 }
 368
 369
 370 //
 371 // EOF
 372 //