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   // fTree = new TTree("T", "T");
 188   // fTree->Branch("forward", &fAODFMD);
 189
 190   PostData(1, fList);
 191   // PostData(2, fTree);
 192 }
 193 //____________________________________________________________________
 194 void
 195 AliForwardMultiplicity::UserExec(Option_t*)
 196 {
 197   // Get the input data
 198   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 199   if (!esd) {
 200     AliWarning("No ESD event found for input event");
 201     return;
 202   }
 203
 204 #if 0
 205   static Int_t nEvents = 0;
 206   nEvents++;
 207   if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
 208 #endif
 209
 210   // On the first event, initialize the parameters
 211   if (fFirstEvent) {
 212     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 213     pars->SetParametersFromESD(esd);
 214     pars->PrintStatus();
 215     fFirstEvent = false;
 216
 217     InitializeSubs();
 218   }
 219   // Clear stuff
 220   fHistos.Clear();
 221   fESDFMD.Clear();
 222   fAODFMD.Clear();
 223
 224   // Read trigger information from the ESD and store in AOD object
 225   UInt_t triggers = 0;
 226   if (!AliForwardUtil::ReadTriggers(esd, triggers, fHTriggers)) {
 227 #ifdef VERBOSE
 228     AliWarning("Failed to read triggers from ESD");
 229 #endif
 230     return;
 231   }
 232   fAODFMD.SetTriggerBits(triggers);
 233
 234   // Mark this event for storage
 235   MarkEventForStore();
 236
 237   // Check if this is a high-flux event
 238   const AliMultiplicity* testmult = esd->GetMultiplicity();
 239   if (!testmult) {
 240 #ifdef VERBOSE
 241     AliWarning("No central multiplicity object found");
 242 #endif
 243     return;
 244   }
 245   Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 246
 247   // Get the FMD ESD data
 248   AliESDFMD* esdFMD = esd->GetFMDData();
 249   if (!esdFMD) {
 250 #ifdef VERBOSE
 251     AliWarning("No FMD data found in ESD");
 252 #endif
 253     return;
 254   }
 255
 256   // Get the vertex information
 257   Double_t vz   = 0;
 258   Bool_t   vzOk = AliForwardUtil::ReadVertex(esd, vz);
 259
 260   fHEventsTr->Fill(vz);
 261   if (!vzOk) {
 262 #ifdef VERBOSE
 263     AliWarning("Failed to read vertex from ESD");
 264 #endif
 265     return;
 266   }
 267   fHEventsTrVtx->Fill(vz);
 268
 269   // Get the vertex bin
 270   Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
 271   fAODFMD.SetIpZ(vz);
 272   if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
 273 #if 0
 274     AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 275                     vz, fHEventsTr->GetXaxis()->GetXmin(),
 276                     fHEventsTr->GetXaxis()->GetXmax()));
 277 #endif
 278     return;
 279   }
 280
 281   // Apply the sharing filter (or hit merging or clustering if you like)
 282   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 283 #ifdef VERBOSE
 284     AliWarning("Sharing filter failed!");
 285 #endif
 286     return;
 287   }
 288
 289   // Calculate the inclusive charged particle density
 290   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 291     AliWarning("Density calculator failed!");
 292     return;
 293   }
 294
 295   // Do the secondary and other corrections.
 296   if (!fCorrections.Correct(fHistos, ivz)) {
 297     AliWarning("Corrections failed");
 298     return;
 299   }
 300
 301   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 302     AliWarning("Histogram collector failed");
 303     return;
 304   }
 305
 306   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 307     fHData->Add(&(fAODFMD.GetHistogram()));
 308 }
 309
 310 //____________________________________________________________________
 311 void
 312 AliForwardMultiplicity::Terminate(Option_t*)
 313 {
 314   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 315   if (!list) {
 316     AliError("No output list defined");
 317     return;
 318   }
 319
 320   // Get our histograms from the container
 321   TH1I* hEventsTr    = static_cast<TH1I*>(list->FindObject("nEventsTr"));
 322   TH1I* hEventsTrVtx = static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 323   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 324
 325
 326   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 327   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 328   TH1D* norm   = hData->ProjectionX("dNdeta", 0, 1,  "");
 329   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 330   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 331   dNdeta->Divide(norm);
 332   dNdeta->SetStats(0);
 333   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 334                 "width");
 335   list->Add(dNdeta);
 336
 337   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 338   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 339   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 340 }
 341
 342 //____________________________________________________________________
 343 void
 344 AliForwardMultiplicity::MarkEventForStore() const
 345 {
 346   // Make sure the AOD tree is filled
 347   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 348   AliAODHandler*      ah =
 349     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 350   if (!ah)
 351     AliFatal("No AOD output handler set in analysis manager");
 352
 353   ah->SetFillAOD(kTRUE);
 354 }
 355
 356 //____________________________________________________________________
 357 void
 358 AliForwardMultiplicity::Print(Option_t*) const
 359 {
 360 }
 361
 362 //
 363 // EOF
 364 //