PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx

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