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     fHData(0),
  18     fFirstEvent(true),
  19     fESDFMD(),
  20     fHistos(),
  21     fAODFMD(),
  22     fEventInspector(),
  23     fEnergyFitter(),
  24     fSharingFilter(),
  25     fDensityCalculator(),
  26     fCorrections(),
  27     fHistCollector(),
  28     fList(0)
  29 {
  30 }
  31
  32 //____________________________________________________________________
  33 AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
  34   : AliAnalysisTaskSE(name),
  35     fHData(0),
  36     fFirstEvent(true),
  37     fESDFMD(),
  38     fHistos(),
  39     fAODFMD(kTRUE),
  40     fEventInspector("event"),
  41     fEnergyFitter("energy"),
  42     fSharingFilter("sharing"),
  43     fDensityCalculator("density"),
  44     fCorrections("corrections"),
  45     fHistCollector("collector"),
  46     fList(0)
  47 {
  48   DefineOutput(1, TList::Class());
  49 }
  50
  51 //____________________________________________________________________
  52 AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
  53   : AliAnalysisTaskSE(o),
  54     fHData(o.fHData),
  55     fFirstEvent(true),
  56     fESDFMD(o.fESDFMD),
  57     fHistos(o.fHistos),
  58     fAODFMD(o.fAODFMD),
  59     fEventInspector(o.fEventInspector),
  60     fEnergyFitter(o.fEnergyFitter),
  61     fSharingFilter(o.fSharingFilter),
  62     fDensityCalculator(o.fDensityCalculator),
  63     fCorrections(o.fCorrections),
  64     fHistCollector(o.fHistCollector),
  65     fList(o.fList)
  66 {
  67 }
  68
  69 //____________________________________________________________________
  70 AliForwardMultiplicity&
  71 AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
  72 {
  73   AliAnalysisTaskSE::operator=(o);
  74
  75   fHData             = o.fHData;
  76   fFirstEvent        = o.fFirstEvent;
  77   fEventInspector    = o.fEventInspector;
  78   fEnergyFitter      = o.fEnergyFitter;
  79   fSharingFilter     = o.fSharingFilter;
  80   fDensityCalculator = o.fDensityCalculator;
  81   fCorrections       = o.fCorrections;
  82   fHistCollector     = o.fHistCollector;
  83   fHistos            = o.fHistos;
  84   fAODFMD            = o.fAODFMD;
  85   fList              = o.fList;
  86
  87   return *this;
  88 }
  89
  90 //____________________________________________________________________
  91 void
  92 AliForwardMultiplicity::SetDebug(Int_t dbg)
  93 {
  94   fEventInspector.SetDebug(dbg);
  95   fEnergyFitter.SetDebug(dbg);
  96   fSharingFilter.SetDebug(dbg);
  97   fDensityCalculator.SetDebug(dbg);
  98   fCorrections.SetDebug(dbg);
  99   fHistCollector.SetDebug(dbg);
 100 }
 101
 102 //____________________________________________________________________
 103 void
 104 AliForwardMultiplicity::Init()
 105 {
 106   fFirstEvent = true;
 107 }
 108
 109 //____________________________________________________________________
 110 void
 111 AliForwardMultiplicity::InitializeSubs()
 112 {
 113   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 114   pars->Init(kTRUE);
 115
 116
 117   TAxis e(pars->GetNetaBins(),  pars->GetEtaMin(),  pars->GetEtaMax());
 118   TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
 119
 120   fHistos.Init(e);
 121   fAODFMD.Init(e);
 122
 123   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 124   fHData->SetStats(0);
 125   fHData->SetDirectory(0);
 126   fList->Add(fHData);
 127
 128   fEnergyFitter.Init(e);
 129   fEventInspector.Init(v);
 130   fHistCollector.Init(v);
 131 }
 132
 133 //____________________________________________________________________
 134 void
 135 AliForwardMultiplicity::UserCreateOutputObjects()
 136 {
 137   fList = new TList;
 138
 139   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 140   AliAODHandler*      ah = dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 141   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 142
 143
 144   TObject* obj = &fAODFMD;
 145   ah->AddBranch("AliAODForwardMult", &obj);
 146
 147   fEventInspector.DefineOutput(fList);
 148   fEnergyFitter.DefineOutput(fList);
 149   fSharingFilter.DefineOutput(fList);
 150   fDensityCalculator.DefineOutput(fList);
 151   fCorrections.DefineOutput(fList);
 152 }
 153 //____________________________________________________________________
 154 void
 155 AliForwardMultiplicity::UserExec(Option_t*)
 156 {
 157   // Get the input data
 158   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 159   if (!esd) {
 160     AliWarning("No ESD event found for input event");
 161     return;
 162   }
 163
 164   // On the first event, initialize the parameters
 165   if (fFirstEvent && esd->GetESDRun()) {
 166     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 167     AliInfo(Form("Initializing with parameters from the ESD:\n"
 168                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 169                  "         AliESDEvent::GetBeamType()     ->%s\n"
 170                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 171                  "         AliESDEvent::GetMagneticField()->%f\n"
 172                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 173                  esd->GetBeamEnergy(),
 174                  esd->GetBeamType(),
 175                  esd->GetCurrentL3(),
 176                  esd->GetMagneticField(),
 177                  esd->GetRunNumber()));
 178     pars->SetParametersFromESD(esd);
 179     pars->PrintStatus();
 180     fFirstEvent = false;
 181
 182     InitializeSubs();
 183   }
 184   // Clear stuff
 185   fHistos.Clear();
 186   fESDFMD.Clear();
 187   fAODFMD.Clear();
 188
 189   Bool_t   lowFlux  = kFALSE;
 190   UInt_t   triggers = 0;
 191   Int_t    ivz      = -1;
 192   Double_t vz       = 0;
 193   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
 194   if (found & AliFMDEventInspector::kNoEvent)    return;
 195   if (found & AliFMDEventInspector::kNoTriggers) return;
 196
 197   // Set trigger bits, and mark this event for storage
 198   fAODFMD.SetTriggerBits(triggers);
 199   MarkEventForStore();
 200
 201   if (found & AliFMDEventInspector::kNoSPD)     return;
 202   if (found & AliFMDEventInspector::kNoFMD)     return;
 203   if (found & AliFMDEventInspector::kNoVertex)  return;
 204   fAODFMD.SetIpZ(vz);
 205
 206   if (found & AliFMDEventInspector::kBadVertex) return;
 207
 208   // Get FMD data
 209   AliESDFMD* esdFMD = esd->GetFMDData();
 210   // Apply the sharing filter (or hit merging or clustering if you like)
 211   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
 212     AliWarning("Sharing filter failed!");
 213     return;
 214   }
 215
 216   // Do the energy stuff
 217   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
 218     AliWarning("Energy fitter failed");
 219     return;
 220   }
 221
 222   // Calculate the inclusive charged particle density
 223   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 224     AliWarning("Density calculator failed!");
 225     return;
 226   }
 227
 228   // Do the secondary and other corrections.
 229   if (!fCorrections.Correct(fHistos, ivz)) {
 230     AliWarning("Corrections failed");
 231     return;
 232   }
 233
 234   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 235     AliWarning("Histogram collector failed");
 236     return;
 237   }
 238
 239   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 240     fHData->Add(&(fAODFMD.GetHistogram()));
 241
 242   PostData(1, fList);
 243 }
 244
 245 //____________________________________________________________________
 246 void
 247 AliForwardMultiplicity::Terminate(Option_t*)
 248 {
 249   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 250   if (!list) {
 251     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 252     if (GetOutputData(1)) GetOutputData(1)->Print();
 253     return;
 254   }
 255
 256   // Get our histograms from the container
 257   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 258   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 259   TH1I* hTriggers    = 0;
 260   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 261                                        hEventsTrVtx, hTriggers)) {
 262     AliError(Form("Didn't get histograms from event selector "
 263                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 264                   hEventsTr, hEventsTrVtx));
 265     list->ls();
 266     return;
 267   }
 268
 269   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
 270   if (!hData) {
 271     AliError(Form("Couldn't get our summed histogram from output "
 272                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
 273     list->ls();
 274     return;
 275   }
 276
 277   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 278   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
 279   TH1D* norm   = hData->ProjectionX("dNdeta", 0, 1,  "");
 280   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 281   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 282   dNdeta->Divide(norm);
 283   dNdeta->SetStats(0);
 284   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
 285                 "width");
 286   list->Add(dNdeta);
 287
 288   fEnergyFitter.Fit(list);
 289   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
 290   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
 291   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
 292 }
 293
 294 //____________________________________________________________________
 295 void
 296 AliForwardMultiplicity::MarkEventForStore() const
 297 {
 298   // Make sure the AOD tree is filled
 299   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 300   AliAODHandler*      ah =
 301     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 302   if (!ah)
 303     AliFatal("No AOD output handler set in analysis manager");
 304
 305   ah->SetFillAOD(kTRUE);
 306 }
 307
 308 //____________________________________________________________________
 309 void
 310 AliForwardMultiplicity::Print(Option_t*) const
 311 {
 312 }
 313
 314 //
 315 // EOF
 316 //