PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx

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