PWG2/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 "AliESDEvent.h"
  28 #include <TROOT.h>
  29 #include <TAxis.h>
  30 #include <THStack.h>
  31 #include <iostream>
  32 #include <iomanip>
  33
  34 //====================================================================
  35 AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name)
  36   : AliAnalysisTaskSE(name),
  37     fEnableLowFlux(false),
  38     fFirstEvent(true),
  39     fCorrManager(0)
  40 {
  41   // Set our persistent pointer
  42   fCorrManager = &AliForwardCorrectionManager::Instance();
  43   fBranchNames =
  44     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  45     "AliESDFMD.,SPDVertex.,PrimaryVertex.";
  46 }
  47
  48 //____________________________________________________________________
  49 AliForwardMultiplicityBase&
  50 AliForwardMultiplicityBase::operator=(const AliForwardMultiplicityBase& o)
  51 {
  52   if (&o == this) return *this;
  53   fEnableLowFlux = o.fEnableLowFlux;
  54   fFirstEvent    = o.fFirstEvent;
  55   fCorrManager   = o.fCorrManager;
  56   return *this;
  57 }
  58
  59 //____________________________________________________________________
  60 Bool_t
  61 AliForwardMultiplicityBase::CheckCorrections(UInt_t what) const
  62 {
  63   //
  64   // Check if all needed corrections are there and accounted for.  If not,
  65   // do a Fatal exit
  66   //
  67   // Parameters:
  68   //    what Which corrections is needed
  69   //
  70   // Return:
  71   //    true if all present, false otherwise
  72   //
  73
  74   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  75   // Check that we have the energy loss fits, needed by
  76   //   AliFMDSharingFilter
  77   //   AliFMDDensityCalculator
  78   if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) {
  79     AliFatal(Form("No energy loss fits"));
  80     return false;
  81   }
  82   // Check that we have the double hit correction - (optionally) used by
  83   //  AliFMDDensityCalculator
  84   if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit()) {
  85     AliFatal("No double hit corrections");
  86     return false;
  87   }
  88   // Check that we have the secondary maps, needed by
  89   //   AliFMDCorrector
  90   //   AliFMDHistCollector
  91   if (what & AliForwardCorrectionManager::kSecondaryMap &&
  92       !fcm.GetSecondaryMap()) {
  93     AliFatal("No secondary corrections");
  94     return false;
  95   }
  96   // Check that we have the vertex bias correction, needed by
  97   //   AliFMDCorrector
  98   if (what & AliForwardCorrectionManager::kVertexBias &&
  99       !fcm.GetVertexBias()) {
 100     AliFatal("No event vertex bias corrections");
 101     return false;
 102   }
 103   // Check that we have the merging efficiencies, optionally used by
 104   //   AliFMDCorrector
 105   if (what & AliForwardCorrectionManager::kMergingEfficiency &&
 106       !fcm.GetMergingEfficiency()) {
 107     AliFatal("No merging efficiencies");
 108     return false;
 109   }
 110   // Check that we have the acceptance correction, needed by
 111   //   AliFMDCorrector
 112   if (what & AliForwardCorrectionManager::kAcceptance &&
 113       !fcm.GetAcceptance()) {
 114     AliFatal("No acceptance corrections");
 115     return false;
 116   }
 117   return true;
 118 }
 119 //____________________________________________________________________
 120 Bool_t
 121 AliForwardMultiplicityBase::ReadCorrections(const TAxis*& pe,
 122                                             const TAxis*& pv,
 123                                             Bool_t        mc)
 124 {
 125   //
 126   // Read corrections
 127   //
 128   //
 129   UInt_t what = AliForwardCorrectionManager::kAll;
 130   if (!fEnableLowFlux)
 131     what ^= AliForwardCorrectionManager::kDoubleHit;
 132   if (!GetCorrections().IsUseVertexBias())
 133     what ^= AliForwardCorrectionManager::kVertexBias;
 134   if (!GetCorrections().IsUseAcceptance())
 135     what ^= AliForwardCorrectionManager::kAcceptance;
 136   if (!GetCorrections().IsUseMergingEfficiency())
 137     what ^= AliForwardCorrectionManager::kMergingEfficiency;
 138
 139   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 140   if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
 141                 GetEventInspector().GetEnergy(),
 142                 GetEventInspector().GetField(),
 143                 mc,
 144                 what)) return false;
 145   if (!CheckCorrections(what)) return false;
 146
 147   // Sett our persistency pointer
 148   // fCorrManager = &fcm;
 149
 150   // Get the eta axis from the secondary maps - if read in
 151   if (!pe) {
 152     pe = fcm.GetEtaAxis();
 153     if (!pe) AliFatal("No eta axis defined");
 154   }
 155   // Get the vertex axis from the secondary maps - if read in
 156   if (!pv) {
 157     pv = fcm.GetVertexAxis();
 158     if (!pv) AliFatal("No vertex axis defined");
 159   }
 160
 161   return true;
 162 }
 163 //____________________________________________________________________
 164 AliESDEvent*
 165 AliForwardMultiplicityBase::GetESDEvent()
 166 {
 167   //
 168   // Get the ESD event. IF this is the first event, initialise
 169   //
 170   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 171   if (!esd) {
 172     AliWarning("No ESD event found for input event");
 173     return 0;
 174   }
 175
 176   // On the first event, initialize the parameters
 177   if (fFirstEvent && esd->GetESDRun()) {
 178     GetEventInspector().ReadRunDetails(esd);
 179
 180     AliInfo(Form("Initializing with parameters from the ESD:\n"
 181                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 182                  "         AliESDEvent::GetBeamType()     ->%s\n"
 183                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 184                  "         AliESDEvent::GetMagneticField()->%f\n"
 185                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 186                  esd->GetBeamEnergy(),
 187                  esd->GetBeamType(),
 188                  esd->GetCurrentL3(),
 189                  esd->GetMagneticField(),
 190                  esd->GetRunNumber()));
 191
 192     fFirstEvent = false;
 193
 194     InitializeSubs();
 195   }
 196   return esd;
 197 }
 198 //____________________________________________________________________
 199 void
 200 AliForwardMultiplicityBase::MarkEventForStore() const
 201 {
 202   // Make sure the AOD tree is filled
 203   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 204   AliAODHandler*      ah =
 205     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 206   if (!ah)
 207     AliFatal("No AOD output handler set in analysis manager");
 208
 209   ah->SetFillAOD(kTRUE);
 210 }
 211
 212 //____________________________________________________________________
 213 void
 214 AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
 215                                            const char*  inName,
 216                                            TList*       output,
 217                                            const char*  outName,
 218                                            Int_t        style) const
 219 {
 220   // Make dN/deta for each ring found in the input list.
 221   //
 222   // A stack of all the dN/deta is also made for easy drawing.
 223   //
 224   // Note, that the distributions are normalised to the number of
 225   // observed events only - they should be corrected for
 226   if (!input) return;
 227   TList* list = static_cast<TList*>(input->FindObject(inName));
 228   if (!list) {
 229     AliWarning(Form("No list %s found in %s", inName, input->GetName()));
 230     return;
 231   }
 232
 233   TList* out = new TList;
 234   out->SetName(outName);
 235   out->SetOwner();
 236   output->Add(out);
 237
 238   THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
 239   const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 240   const char** ptr         = names;
 241
 242   while (*ptr) {
 243     TList* thisList = new TList;
 244     thisList->SetOwner();
 245     thisList->SetName(*ptr);
 246     out->Add(thisList);
 247
 248     TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
 249     if (!h) {
 250       AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
 251       ptr++;
 252       continue;
 253     }
 254     TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
 255     copy->SetDirectory(0);
 256     thisList->Add(copy);
 257
 258     TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
 259     for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 260       for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
 261         Double_t c = copy->GetBinContent(i, j);
 262         Double_t e = copy->GetBinError(i, j);
 263         Double_t a = norm->GetBinContent(i);
 264         copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
 265         copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
 266       }
 267     }
 268
 269     TH1D* res  =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
 270                                                      h->GetNbinsY(),"e"));
 271     TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
 272     res->SetTitle(*ptr);
 273     res->Scale(1., "width");
 274     copy->Scale(1., "width");
 275     proj->Scale(1. / norm->GetMaximum(), "width");
 276     norm->Scale(1. / norm->GetMaximum());
 277
 278     res->SetMarkerStyle(style);
 279     norm->SetDirectory(0);
 280     res->SetDirectory(0);
 281     proj->SetDirectory(0);
 282     thisList->Add(norm);
 283     thisList->Add(res);
 284     thisList->Add(proj);
 285     dndetaRings->Add(res);
 286     ptr++;
 287   }
 288   out->Add(dndetaRings);
 289 }
 290
 291 //____________________________________________________________________
 292 void
 293 AliForwardMultiplicityBase::Print(Option_t* option) const
 294 {
 295   //
 296   // Print information
 297   //
 298   // Parameters:
 299   //    option Not used
 300   //
 301
 302   std::cout << ClassName() << ": " << GetName() << "\n"
 303             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 304             << "\n"
 305             << "  Off-line trigger mask:  0x"
 306             << std::hex     << std::setfill('0')
 307             << std::setw (8) << fOfflineTriggerMask
 308             << std::dec     << std::setfill (' ') << std::endl;
 309   gROOT->IncreaseDirLevel();
 310   if (fCorrManager) fCorrManager->Print();
 311   else
 312     std::cout << "  Correction manager not set yet" << std::endl;
 313   GetEventInspector()   .Print(option);
 314   GetSharingFilter()    .Print(option);
 315   GetDensityCalculator().Print(option);
 316   GetCorrections()      .Print(option);
 317   GetHistCollector()    .Print(option);
 318   gROOT->DecreaseDirLevel();
 319 }
 320
 321
 322 //
 323 // EOF
 324 //