PWG2/FORWARD/analysis2/AliForwardMultiplicity.cxx

   1 #include "AliForwardMultiplicity.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 <TH1.h>
  11 #include <TDirectory.h>
  12 #include <TTree.h>
  13
  14 //====================================================================
  15 AliForwardMultiplicity::AliForwardMultiplicity()
  16   : AliAnalysisTaskSE(),
  17     fHEventsTr(0),
  18     fHEventsTrVtx(0),
  19     fHTriggers(0),
  20     fHData(0),
  21     fFirstEvent(true),
  22     fLowFluxCut(1000),
  23     fESDFMD(),
  24     fHistos(),
  25     fAODFMD(),
  26     fSharingFilter(),
  27     fDensityCalculator(),
  28     fCorrections(),
  29     fHistCollector(),
  30     fList(0),
  31     fTree(0)
  32 {
  33 }
  34
  35 //____________________________________________________________________
  36 AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
  37   : AliAnalysisTaskSE(name),
  38     fHEventsTr(0),
  39     fHEventsTrVtx(0),
  40     fHTriggers(0),
  41     fHData(0),
  42     fFirstEvent(true),
  43     fLowFluxCut(1000),
  44     fESDFMD(),
  45     fHistos(),
  46     fAODFMD(kTRUE),
  47     fSharingFilter("sharing"),
  48     fDensityCalculator("density"),
  49     fCorrections("corrections"),
  50     fHistCollector("collector"),
  51     fList(0),
  52     fTree(0)
  53 {
  54   DefineOutput(1, TList::Class());
  55   // DefineOutput(2, TTree::Class());
  56 }
  57
  58 //____________________________________________________________________
  59 AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
  60   : AliAnalysisTaskSE(o),
  61     fHEventsTr(o.fHEventsTr),
  62     fHEventsTrVtx(o.fHEventsTrVtx),
  63     fHTriggers(o.fHTriggers),
  64     fHData(o.fHData),
  65     fFirstEvent(true),
  66     fLowFluxCut(1000),
  67     fESDFMD(o.fESDFMD),
  68     fHistos(o.fHistos),
  69     fAODFMD(o.fAODFMD),
  70     fSharingFilter(o.fSharingFilter),
  71     fDensityCalculator(o.fDensityCalculator),
  72     fCorrections(o.fCorrections),
  73     fHistCollector(o.fHistCollector),
  74     fList(o.fList),
  75     fTree(o.fTree)
  76 {
  77 }
  78
  79 //____________________________________________________________________
  80 AliForwardMultiplicity&
  81 AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
  82 {
  83   fHEventsTr         = o.fHEventsTr;
  84   fHEventsTrVtx      = o.fHEventsTrVtx;
  85   fHTriggers         = o.fHTriggers;
  86   fHData             = o.fHData;
  87   fFirstEvent        = o.fFirstEvent;
  88   fSharingFilter     = o.fSharingFilter;
  89   fDensityCalculator = o.fDensityCalculator;
  90   fCorrections       = o.fCorrections;
  91   fHistCollector     = o.fHistCollector;
  92   fHistos            = o.fHistos;
  93   fAODFMD            = o.fAODFMD;
  94   fList              = o.fList;
  95   fTree              = o.fTree;
  96
  97   return *this;
  98 }
  99
 100 //____________________________________________________________________
 101 void
 102 AliForwardMultiplicity::Init()
 103 {
 104   fFirstEvent = true;
 105 }
 106
 107 //____________________________________________________________________
 108 void
 109 AliForwardMultiplicity::InitializeSubs()
 110 {
 111   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 112   pars->Init(kTRUE);
 113
 114   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
 115                       pars->GetNvtxBins(),
 116                       -pars->GetVtxCutZ(),
 117                       pars->GetVtxCutZ());
 118   fHEventsTr->SetXTitle("v_{z} [cm]");
 119   fHEventsTr->SetYTitle("# of events");
 120   fHEventsTr->SetFillColor(kRed+1);
 121   fHEventsTr->SetFillStyle(3001);
 122   fHEventsTr->SetDirectory(0);
 123   fList->Add(fHEventsTr);
 124   // fHEventsTr->Sumw2();
 125
 126   fHEventsTrVtx = new TH1I("nEventsTrVtx",
 127                            "Number of events w/trigger and vertex",
 128                            pars->GetNvtxBins(),
 129                            -pars->GetVtxCutZ(),
 130                            pars->GetVtxCutZ());
 131   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 132   fHEventsTrVtx->SetYTitle("# of events");
 133   fHEventsTrVtx->SetFillColor(kBlue+1);
 134   fHEventsTrVtx->SetFillStyle(3001);
 135   fHEventsTrVtx->SetDirectory(0);
 136   fList->Add(fHEventsTrVtx);
 137   // fHEventsTrVtx->Sumw2();
 138
 139
 140   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
 141   fHTriggers->SetFillColor(kRed+1);
 142   fHTriggers->SetFillStyle(3001);
 143   fHTriggers->SetStats(0);
 144   fHTriggers->SetDirectory(0);
 145   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
 146   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
 147   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
 148   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
 149   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
 150   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
 151   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
 152   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
 153   fList->Add(fHTriggers);
 154
 155   TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
 156   fHistos.Init(e);
 157   fAODFMD.Init(e);
 158
 159   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 160   fHData->SetStats(0);
 161   fHData->SetDirectory(0);
 162   fList->Add(fHData);
 163
 164   fHistCollector.Init(*(fHEventsTr->GetXaxis()));
 165 }
 166
 167 //____________________________________________________________________
 168 void
 169 AliForwardMultiplicity::UserCreateOutputObjects()
 170 {
 171   fList = new TList;
 172
 173   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 174   AliAODHandler*      ah =
 175     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 176   if (!ah)
 177     AliFatal("No AOD output handler set in analysis manager");
 178
 179
 180   TObject* obj = &fAODFMD;
 181   ah->AddBranch("AliAODForwardMult", &obj);
 182
 183   fSharingFilter.DefineOutput(fList);
 184   fDensityCalculator.DefineOutput(fList);
 185   fCorrections.DefineOutput(fList);
 186 }
 187 //____________________________________________________________________
 188 void
 189 AliForwardMultiplicity::UserExec(Option_t*)
 190 {
 191   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 192
 193   // Get the input data
 194   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 195   if (!esd) {
 196     AliWarning("No ESD event found for input event");
 197     return;
 198   }
 199
 200 #if 0
 201   static Int_t nEvents = 0;
 202   nEvents++;
 203   if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
 204 #endif
 205
 206   // On the first event, initialize the parameters
 207   if (fFirstEvent) {
 208     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 209     pars->SetParametersFromESD(esd);
 210     pars->PrintStatus();
 211     fFirstEvent = false;
 212
 213     InitializeSubs();
 214   }
 215   // Clear stuff
 216   fHistos.Clear();
 217   fESDFMD.Clear();
 218   fAODFMD.Clear();
 219
 220   // Read trigger information from the ESD and store in AOD object
 221   UInt_t triggers = 0;
 222   if (!AliForwardUtil::ReadTriggers(esd, triggers, fHTriggers)) {
 223     if (am->GetDebugLevel() > 2)
 224       AliWarning("Failed to read triggers from ESD");
 225     return;
 226   }
 227   fAODFMD.SetTriggerBits(triggers);
 228
 229   // Mark this event for storage
 230   MarkEventForStore();
 231
 232   // Check if this is a high-flux event
 233   const AliMultiplicity* testmult = esd->GetMultiplicity();
 234   if (!testmult) {
 235     if (am->GetDebugLevel() > 3)
 236       AliWarning("No central multiplicity object found");
 237     return;
 238   }
 239   Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 240
 241   // Get the FMD ESD data
 242   AliESDFMD* esdFMD = esd->GetFMDData();
 243   if (!esdFMD) {
 244     if (am->GetDebugLevel() > 3)
 245       AliWarning("No FMD data found in ESD");
 246     return;
 247   }
 248
 249   // Get the vertex information
 250   Double_t vz   = 0;
 251   Bool_t   vzOk = AliForwardUtil::ReadVertex(esd, vz);
 252
 253   fHEventsTr->Fill(vz);
 254   if (!vzOk) {
 255     if (am->GetDebugLevel() > 3)
 256       AliWarning("Failed to read vertex from ESD");
 257     return;
 258   }
 259   fHEventsTrVtx->Fill(vz);
 260
 261   // Get the vertex bin
 262   Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
 263   fAODFMD.SetIpZ(vz);
 264   if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
 265     if (am->GetDebugLevel() > 3)
 266       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 267                       vz, fHEventsTr->GetXaxis()->GetXmin(),
 268                       fHEventsTr->GetXaxis()->GetXmax()));
 269     return;
 270   }
 271   if (am->GetDebugLevel() > 1)
 272     AliInfo(Form("Event has %d SPD tracklets, cut is %d, this is a %s event",
 273                  testmult->GetNumberOfTracklets(),
 274                  fLowFluxCut, (lowFlux ? "low" : "high")));
 275   if (am->GetDebugLevel() > 1)
 276     AliInfo(Form("Events vertex @ %f (bin %d), in range", vz, ivz));
 277
 278
 279   // Apply the sharing filter (or hit merging or clustering if you like)
 280   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 281     AliWarning("Sharing filter failed!");
 282     return;
 283   }
 284
 285   // Calculate the inclusive charged particle density
 286   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 287     AliWarning("Density calculator failed!");
 288     return;
 289   }
 290
 291   // Do the secondary and other corrections.
 292   if (!fCorrections.Correct(fHistos, ivz)) {
 293     AliWarning("Corrections failed");
 294     return;
 295   }
 296
 297   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 298     AliWarning("Histogram collector failed");
 299     return;
 300   }
 301
 302   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 303     fHData->Add(&(fAODFMD.GetHistogram()));
 304
 305   PostData(1, fList);
 306 }
 307
 308 //____________________________________________________________________
 309 void
 310 AliForwardMultiplicity::Terminate(Option_t*)
 311 {
 312   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 313   if (!list) {
 314     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 315     if (GetOutputData(1)) GetOutputData(1)->Print();
 316     return;
 317   }
 318
 319   // Get our histograms from the container
 320   TH1I* hEventsTr    = static_cast<TH1I*>(list->FindObject("nEventsTr"));
 321   TH1I* hEventsTrVtx = static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 322   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 323   if (!hData || !hEventsTr || !hEventsTrVtx) {
 324     AliError(Form("one or more histograms could not be found in output "
 325                   "list %s (hEventsTr=%p,hEventsTrVtx=%p,d2Ndetadphi=%p)",
 326                   list->GetName(), hEventsTr, hEventsTrVtx, hData));
 327     list->ls();
 328     return;
 329   }
 330
 331   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 332   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 333   TH1D* norm   = hData->ProjectionX("dNdeta", 0, 1,  "");
 334   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 335   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 336   dNdeta->Divide(norm);
 337   dNdeta->SetStats(0);
 338   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 339                 "width");
 340   list->Add(dNdeta);
 341
 342   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 343   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 344   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 345 }
 346
 347 //____________________________________________________________________
 348 void
 349 AliForwardMultiplicity::MarkEventForStore() const
 350 {
 351   // Make sure the AOD tree is filled
 352   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 353   AliAODHandler*      ah =
 354     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 355   if (!ah)
 356     AliFatal("No AOD output handler set in analysis manager");
 357
 358   ah->SetFillAOD(kTRUE);
 359 }
 360
 361 //____________________________________________________________________
 362 void
 363 AliForwardMultiplicity::Print(Option_t*) const
 364 {
 365 }
 366
 367 //
 368 // EOF
 369 //