PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx

   1 #include "AliForwardMultiplicityTask.h"
   2 #include "AliTriggerAnalysis.h"
   3 #include "AliPhysicsSelection.h"
   4 #include "AliLog.h"
   5 // #include "AliFMDAnaParameters.h"
   6 #include "AliESDEvent.h"
   7 #include "AliAODHandler.h"
   8 #include "AliMultiplicity.h"
   9 #include "AliInputEventHandler.h"
  10 #include "AliForwardCorrectionManager.h"
  11 #include "AliAnalysisManager.h"
  12 #include <TH1.h>
  13 #include <TDirectory.h>
  14 #include <TTree.h>
  15 #include <TROOT.h>
  16 #include <iostream>
  17 #include <iomanip>
  18
  19 //====================================================================
  20 AliForwardMultiplicityTask::AliForwardMultiplicityTask()
  21   : AliAnalysisTaskSE(),
  22     fEnableLowFlux(true),
  23     fHData(0),
  24     fFirstEvent(true),
  25     fESDFMD(),
  26     fHistos(),
  27     fAODFMD(),
  28     fEventInspector(),
  29     fEnergyFitter(),
  30     fSharingFilter(),
  31     fDensityCalculator(),
  32     fCorrections(),
  33     fHistCollector(),
  34     fList(0)
  35 {
  36 }
  37
  38 //____________________________________________________________________
  39 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
  40   : AliAnalysisTaskSE(name),
  41     fEnableLowFlux(true),
  42     fHData(0),
  43     fFirstEvent(true),
  44     fESDFMD(),
  45     fHistos(),
  46     fAODFMD(kTRUE),
  47     fEventInspector("event"),
  48     fEnergyFitter("energy"),
  49     fSharingFilter("sharing"),
  50     fDensityCalculator("density"),
  51     fCorrections("corrections"),
  52     fHistCollector("collector"),
  53     fList(0)
  54 {
  55   DefineOutput(1, TList::Class());
  56 }
  57
  58 //____________________________________________________________________
  59 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
  60   : AliAnalysisTaskSE(o),
  61     fEnableLowFlux(o.fEnableLowFlux),
  62     fHData(o.fHData),
  63     fFirstEvent(true),
  64     fESDFMD(o.fESDFMD),
  65     fHistos(o.fHistos),
  66     fAODFMD(o.fAODFMD),
  67     fEventInspector(o.fEventInspector),
  68     fEnergyFitter(o.fEnergyFitter),
  69     fSharingFilter(o.fSharingFilter),
  70     fDensityCalculator(o.fDensityCalculator),
  71     fCorrections(o.fCorrections),
  72     fHistCollector(o.fHistCollector),
  73     fList(o.fList)
  74 {
  75   DefineOutput(1, TList::Class());
  76 }
  77
  78 //____________________________________________________________________
  79 AliForwardMultiplicityTask&
  80 AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
  81 {
  82   AliAnalysisTaskSE::operator=(o);
  83
  84   fEnableLowFlux     = o.fEnableLowFlux;
  85   fHData             = o.fHData;
  86   fFirstEvent        = o.fFirstEvent;
  87   fEventInspector    = o.fEventInspector;
  88   fEnergyFitter      = o.fEnergyFitter;
  89   fSharingFilter     = o.fSharingFilter;
  90   fDensityCalculator = o.fDensityCalculator;
  91   fCorrections       = o.fCorrections;
  92   fHistCollector     = o.fHistCollector;
  93   fHistos            = o.fHistos;
  94   fAODFMD            = o.fAODFMD;
  95   fList              = o.fList;
  96
  97   return *this;
  98 }
  99
 100 //____________________________________________________________________
 101 void
 102 AliForwardMultiplicityTask::SetDebug(Int_t dbg)
 103 {
 104   fEventInspector.SetDebug(dbg);
 105   fEnergyFitter.SetDebug(dbg);
 106   fSharingFilter.SetDebug(dbg);
 107   fDensityCalculator.SetDebug(dbg);
 108   fCorrections.SetDebug(dbg);
 109   fHistCollector.SetDebug(dbg);
 110 }
 111
 112 //____________________________________________________________________
 113 void
 114 AliForwardMultiplicityTask::Init()
 115 {
 116   fFirstEvent = true;
 117 }
 118
 119 //____________________________________________________________________
 120 void
 121 AliForwardMultiplicityTask::InitializeSubs()
 122 {
 123
 124   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 125   // pars->Init(kTRUE);
 126
 127   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 128   fcm.Init(fEventInspector.GetCollisionSystem(),
 129            fEventInspector.GetEnergy(),
 130            fEventInspector.GetField());
 131   // Check that we have the energy loss fits, needed by
 132   //   AliFMDSharingFilter
 133   //   AliFMDDensityCalculator
 134   if (!fcm.GetELossFit()) {
 135     AliFatal(Form("No energy loss fits"));
 136     return;
 137   }
 138   // Check that we have the double hit correction - (optionally) used by
 139   //  AliFMDDensityCalculator
 140   if (!fcm.GetDoubleHit()) {
 141     AliWarning("No double hit corrections");
 142   }
 143   // Check that we have the secondary maps, needed by
 144   //   AliFMDCorrections
 145   //   AliFMDHistCollector
 146   if (!fcm.GetSecondaryMap()) {
 147     AliFatal("No secondary corrections");
 148     return;
 149   }
 150   // Check that we have the vertex bias correction, needed by
 151   //   AliFMDCorrections
 152   if (!fcm.GetVertexBias()) {
 153     AliFatal("No event vertex bias corrections");
 154     return;
 155   }
 156   // Check that we have the merging efficiencies, optionally used by
 157   //   AliFMDCorrections
 158   if (!fcm.GetMergingEfficiency()) {
 159     AliWarning("No merging efficiencies");
 160   }
 161
 162
 163   const TAxis* pe = fcm.GetEtaAxis();
 164   const TAxis* pv = fcm.GetVertexAxis();
 165   if (!pe) AliFatal("No eta axis defined");
 166   if (!pv) AliFatal("No vertex axis defined");
 167
 168   // TAxis e(pars->GetNetaBins(),  pars->GetEtaMin(),  pars->GetEtaMax());
 169   // TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
 170   // pv->Dump();
 171
 172   fHistos.Init(*pe);
 173   fAODFMD.Init(*pe);
 174
 175   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 176   fHData->SetStats(0);
 177   fHData->SetDirectory(0);
 178   fList->Add(fHData);
 179
 180   fEnergyFitter.Init(*pe);
 181   fEventInspector.Init(*pv);
 182   fHistCollector.Init(*pv);
 183
 184   this->Print();
 185 }
 186
 187 //____________________________________________________________________
 188 void
 189 AliForwardMultiplicityTask::UserCreateOutputObjects()
 190 {
 191   fList = new TList;
 192
 193   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 194   AliAODHandler*      ah =
 195     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 196   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 197
 198
 199   TObject* obj = &fAODFMD;
 200   ah->AddBranch("AliAODForwardMult", &obj);
 201
 202   fEventInspector.DefineOutput(fList);
 203   fEnergyFitter.DefineOutput(fList);
 204   fSharingFilter.DefineOutput(fList);
 205   fDensityCalculator.DefineOutput(fList);
 206   fCorrections.DefineOutput(fList);
 207
 208   PostData(1, fList);
 209 }
 210 //____________________________________________________________________
 211 void
 212 AliForwardMultiplicityTask::UserExec(Option_t*)
 213 {
 214   // static Int_t cnt = 0;
 215   // cnt++;
 216   // Get the input data
 217   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 218   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
 219   if (!esd) {
 220     AliWarning("No ESD event found for input event");
 221     return;
 222   }
 223
 224   // On the first event, initialize the parameters
 225   if (fFirstEvent && esd->GetESDRun()) {
 226     fEventInspector.ReadRunDetails(esd);
 227
 228     AliInfo(Form("Initializing with parameters from the ESD:\n"
 229                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 230                  "         AliESDEvent::GetBeamType()     ->%s\n"
 231                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 232                  "         AliESDEvent::GetMagneticField()->%f\n"
 233                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 234                  esd->GetBeamEnergy(),
 235                  esd->GetBeamType(),
 236                  esd->GetCurrentL3(),
 237                  esd->GetMagneticField(),
 238                  esd->GetRunNumber()));
 239
 240
 241
 242     // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 243     // pars->SetParametersFromESD(esd);
 244     // pars->PrintStatus();
 245     fFirstEvent = false;
 246
 247     InitializeSubs();
 248   }
 249   // Clear stuff
 250   fHistos.Clear();
 251   fESDFMD.Clear();
 252   fAODFMD.Clear();
 253
 254   Bool_t   lowFlux  = kFALSE;
 255   UInt_t   triggers = 0;
 256   UShort_t ivz      = 0;
 257   Double_t vz       = 0;
 258   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
 259   if (found & AliFMDEventInspector::kNoEvent)    return;
 260   if (found & AliFMDEventInspector::kNoTriggers) return;
 261
 262   // Set trigger bits, and mark this event for storage
 263   fAODFMD.SetTriggerBits(triggers);
 264   MarkEventForStore();
 265
 266   if (found & AliFMDEventInspector::kNoSPD)     return;
 267   if (found & AliFMDEventInspector::kNoFMD)     return;
 268   if (found & AliFMDEventInspector::kNoVertex)  return;
 269   fAODFMD.SetIpZ(vz);
 270
 271   if (found & AliFMDEventInspector::kBadVertex) return;
 272
 273   // We we do not want to use low flux specific code, we disable it here.
 274   if (!fEnableLowFlux) lowFlux = false;
 275
 276   // Get FMD data
 277   AliESDFMD* esdFMD = esd->GetFMDData();
 278   // Apply the sharing filter (or hit merging or clustering if you like)
 279   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
 280     AliWarning("Sharing filter failed!");
 281     return;
 282   }
 283
 284   // Do the energy stuff
 285   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
 286     AliWarning("Energy fitter failed");
 287     return;
 288   }
 289
 290   // Calculate the inclusive charged particle density
 291   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 292     AliWarning("Density calculator failed!");
 293     return;
 294   }
 295
 296   // Do the secondary and other corrections.
 297   if (!fCorrections.Correct(fHistos, ivz)) {
 298     AliWarning("Corrections failed");
 299     return;
 300   }
 301
 302   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 303     AliWarning("Histogram collector failed");
 304     return;
 305   }
 306
 307   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 308     fHData->Add(&(fAODFMD.GetHistogram()));
 309
 310   PostData(1, fList);
 311 }
 312
 313 //____________________________________________________________________
 314 void
 315 AliForwardMultiplicityTask::Terminate(Option_t*)
 316 {
 317   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 318   if (!list) {
 319     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 320     if (GetOutputData(1)) GetOutputData(1)->Print();
 321     return;
 322   }
 323
 324   // Get our histograms from the container
 325   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 326   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 327   TH1I* hTriggers    = 0;
 328   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 329                                        hEventsTrVtx, hTriggers)) {
 330     AliError(Form("Didn't get histograms from event selector "
 331                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 332                   hEventsTr, hEventsTrVtx));
 333     list->ls();
 334     return;
 335   }
 336
 337   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 338   if (!hData) {
 339     AliError(Form("Couldn't get our summed histogram from output "
 340                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 341     list->ls();
 342     return;
 343   }
 344
 345   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 346   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 347   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
 348   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 349   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 350   dNdeta->Divide(norm);
 351   dNdeta->SetStats(0);
 352   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 353                 "width");
 354   list->Add(dNdeta);
 355   list->Add(norm);
 356
 357   fEnergyFitter.Fit(list);
 358   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 359   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 360   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 361 }
 362
 363 //____________________________________________________________________
 364 void
 365 AliForwardMultiplicityTask::MarkEventForStore() const
 366 {
 367   // Make sure the AOD tree is filled
 368   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 369   AliAODHandler*      ah =
 370     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 371   if (!ah)
 372     AliFatal("No AOD output handler set in analysis manager");
 373
 374   ah->SetFillAOD(kTRUE);
 375 }
 376
 377 //____________________________________________________________________
 378 void
 379 AliForwardMultiplicityTask::Print(Option_t* option) const
 380 {
 381   std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n"
 382             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 383             << std::endl;
 384   gROOT->IncreaseDirLevel();
 385   fEventInspector   .Print(option);
 386   fEnergyFitter     .Print(option);
 387   fSharingFilter    .Print(option);
 388   fDensityCalculator.Print(option);
 389   fCorrections      .Print(option);
 390   fHistCollector    .Print(option);
 391   gROOT->DecreaseDirLevel();
 392 }
 393
 394 //
 395 // EOF
 396 //