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("nEvents", "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   // fHEventsTr->Sumw2();
 124
 125   fHEventsTrVtx = new TH1I("nEventsTrVtx",
 126                            "Number of events w/trigger and vertex",
 127                            pars->GetNvtxBins(),
 128                            -pars->GetVtxCutZ(),
 129                            pars->GetVtxCutZ());
 130   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 131   fHEventsTrVtx->SetYTitle("# of events");
 132   fHEventsTrVtx->SetFillColor(kBlue+1);
 133   fHEventsTrVtx->SetFillStyle(3001);
 134   fHEventsTrVtx->SetDirectory(0);
 135   // fHEventsTrVtx->Sumw2();
 136
 137
 138   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
 139   fHTriggers->SetFillColor(kRed+1);
 140   fHTriggers->SetFillStyle(3001);
 141   fHTriggers->SetStats(0);
 142   fHTriggers->SetDirectory(0);
 143   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
 144   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
 145   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
 146   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
 147   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
 148   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
 149   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
 150   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
 151
 152   TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
 153   fHistos.Init(e);
 154   fAODFMD.Init(e);
 155
 156   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
 157   fHData->SetStats(0);
 158   fHData->SetDirectory(0);
 159   fSharingFilter.Init();
 160   fHistCollector.Init(*(fHEventsTr->GetXaxis()));
 161 }
 162
 163 //____________________________________________________________________
 164 void
 165 AliForwardMultiplicity::UserCreateOutputObjects()
 166 {
 167   fList = new TList;
 168
 169   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 170   AliAODHandler*      ah =
 171     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 172   if (!ah)
 173     AliFatal("No AOD output handler set in analysis manager");
 174
 175
 176   TObject* obj = &fAODFMD;
 177   ah->AddBranch("AliAODForwardMult", &obj);
 178
 179   // fTree = new TTree("T", "T");
 180   // fTree->Branch("forward", &fAODFMD);
 181
 182   PostData(1, fList);
 183   // PostData(2, fTree);
 184 }
 185 //____________________________________________________________________
 186 void
 187 AliForwardMultiplicity::UserExec(Option_t*)
 188 {
 189   // Get the input data
 190   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 191   if (!esd) {
 192     AliWarning("No ESD event found for input event");
 193     return;
 194   }
 195
 196 #if 0
 197   static Int_t nEvents = 0;
 198   nEvents++;
 199   if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
 200 #endif
 201
 202   // On the first event, initialize the parameters
 203   if (fFirstEvent) {
 204     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 205     pars->SetParametersFromESD(esd);
 206     pars->PrintStatus();
 207     fFirstEvent = false;
 208
 209     InitializeSubs();
 210   }
 211   // Clear stuff
 212   fHistos.Clear();
 213   fESDFMD.Clear();
 214   fAODFMD.Clear();
 215
 216   // Read trigger information from the ESD and store in AOD object
 217   if (!ReadTriggers(esd)) {
 218 #ifdef VERBOSE
 219     AliWarning("Failed to read triggers from ESD");
 220 #endif
 221     return;
 222   }
 223
 224   // Mark this event for storage
 225   MarkEventForStore();
 226
 227   // Check if this is a high-flux event
 228   const AliMultiplicity* testmult = esd->GetMultiplicity();
 229   if (!testmult) {
 230 #ifdef VERBOSE
 231     AliWarning("No central multiplicity object found");
 232 #endif
 233     return;
 234   }
 235   Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 236
 237   // Get the FMD ESD data
 238   AliESDFMD* esdFMD = esd->GetFMDData();
 239   if (!esdFMD) {
 240 #ifdef VERBOSE
 241     AliWarning("No FMD data found in ESD");
 242 #endif
 243     return;
 244   }
 245
 246   // Get the vertex information
 247   Double_t vz   = 0;
 248   Bool_t   vzOk = ReadVertex(esd, vz);
 249
 250   fHEventsTr->Fill(vz);
 251   if (!vzOk) {
 252 #ifdef VERBOSE
 253     AliWarning("Failed to read vertex from ESD");
 254 #endif
 255     return;
 256   }
 257   fHEventsTrVtx->Fill(vz);
 258
 259   // Get the vertex bin
 260   Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
 261   fAODFMD.SetIpZ(vz);
 262   if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
 263 #if 0
 264     AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 265                     vz, fHEventsTr->GetXaxis()->GetXmin(),
 266                     fHEventsTr->GetXaxis()->GetXmax()));
 267 #endif
 268     return;
 269   }
 270
 271   // Apply the sharing filter (or hit merging or clustering if you like)
 272   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
 273 #ifdef VERBOSE
 274     AliWarning("Sharing filter failed!");
 275 #endif
 276     return;
 277   }
 278
 279   // Calculate the inclusive charged particle density
 280   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 281     AliWarning("Density calculator failed!");
 282     return;
 283   }
 284
 285   // Do the secondary and other corrections.
 286   if (!fCorrections.Correct(fHistos, ivz)) {
 287     AliWarning("Corrections failed");
 288     return;
 289   }
 290
 291   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
 292     AliWarning("Histogram collector failed");
 293     return;
 294   }
 295
 296   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
 297     fHData->Add(&(fAODFMD.GetHistogram()));
 298 }
 299
 300 //____________________________________________________________________
 301 void
 302 AliForwardMultiplicity::Terminate(Option_t*)
 303 {
 304   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 305   if (!list) {
 306     AliError("No output list defined");
 307     return;
 308   }
 309   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
 310   TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
 311   TH1D* norm   = fHData->ProjectionX("dNdeta", 0, 1,  "");
 312   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
 313   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 314   dNdeta->Divide(norm);
 315   dNdeta->SetStats(0);
 316   dNdeta->Scale(Double_t(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
 317                 "width");
 318
 319   list->Add(fHEventsTr);
 320   list->Add(fHEventsTrVtx);
 321   list->Add(fHTriggers);
 322   list->Add(fHData);
 323   list->Add(dNdeta);
 324
 325   TList* last = new TList;
 326   last->SetName("LastEvent");
 327   list->Add(last);
 328   last->Add(&fAODFMD.GetHistogram());
 329   last->Add(fHistos.fFMD1i);
 330   last->Add(fHistos.fFMD2i);
 331   last->Add(fHistos.fFMD2o);
 332   last->Add(fHistos.fFMD3i);
 333   last->Add(fHistos.fFMD3o);
 334
 335
 336   fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
 337   fSharingFilter.Output(list);
 338
 339   fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
 340   fDensityCalculator.Output(list);
 341
 342   fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
 343   fCorrections.Output(list);
 344 }
 345
 346 //____________________________________________________________________
 347 void
 348 AliForwardMultiplicity::MarkEventForStore() const
 349 {
 350   // Make sure the AOD tree is filled
 351   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 352   AliAODHandler*      ah =
 353     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 354   if (!ah)
 355     AliFatal("No AOD output handler set in analysis manager");
 356
 357   ah->SetFillAOD(kTRUE);
 358 }
 359 //____________________________________________________________________
 360 Bool_t
 361 AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
 362 {
 363   // Get the analysis manager - should always be there
 364   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 365   if (!am) {
 366     AliWarning("No analysis manager defined!");
 367     return kFALSE;
 368   }
 369
 370   // Get the input handler - should always be there
 371   AliInputEventHandler* ih =
 372     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 373   if (!ih) {
 374     AliWarning("No input handler");
 375     return kFALSE;
 376   }
 377
 378   // Get the physics selection - add that by using the macro
 379   // AddTaskPhysicsSelection.C
 380   AliPhysicsSelection* ps =
 381     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
 382   if (!ps) {
 383     AliWarning("No physics selection");
 384     return kFALSE;
 385   }
 386
 387   // Check if this is a collision candidate (INEL)
 388   Bool_t inel = ps->IsCollisionCandidate(esd);
 389   if (inel) {
 390     fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
 391     fHTriggers->Fill(.5);
 392   }
 393
 394
 395   // IF this is inel, see if we have a tracklet
 396   if (inel) {
 397     const AliMultiplicity* spdmult = esd->GetMultiplicity();
 398     if (!spdmult) {
 399       AliWarning("No SPD multiplicity");
 400     }
 401     else {
 402       Int_t n = spdmult->GetNumberOfTracklets();
 403       for (Int_t j = 0; j < n; j++) {
 404         if(TMath::Abs(spdmult->GetEta(j)) < 1) {
 405           fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
 406           fHTriggers->Fill(1.5);
 407           break;
 408         }
 409       }
 410     }
 411   }
 412
 413   // Analyse some trigger stuff
 414   AliTriggerAnalysis ta;
 415   if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
 416     fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
 417     fHTriggers->Fill(2.5);
 418   }
 419
 420   // Get trigger stuff
 421   TString triggers = esd->GetFiredTriggerClasses();
 422   if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
 423     fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
 424     fHTriggers->Fill(3.5);
 425   }
 426
 427   if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
 428     fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
 429     fHTriggers->Fill(4.5);
 430   }
 431
 432   if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
 433     fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
 434     fHTriggers->Fill(5.5);
 435   }
 436
 437
 438   if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
 439     fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
 440     fHTriggers->Fill(6.5);
 441   }
 442
 443   if (triggers.Contains("CINT1-E-NOPF-ALL")) {
 444     fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
 445     fHTriggers->Fill(7.5);
 446   }
 447
 448   return kTRUE;
 449 }
 450 //____________________________________________________________________
 451 Bool_t
 452 AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
 453 {
 454   // Get the vertex
 455   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
 456   if (!vertex) {
 457 #ifdef VERBOSE
 458     AliWarning("No SPD vertex found in ESD");
 459 #endif
 460     return kFALSE;
 461   }
 462
 463   // Check that enough tracklets contributed
 464   if(vertex->GetNContributors() <= 0) {
 465 #ifdef VERBOSE
 466     AliWarning(Form("Number of contributors to vertex is %d<=0",
 467                     vertex->GetNContributors()));
 468 #endif
 469     return kFALSE;
 470   }
 471
 472   // Check that the uncertainty isn't too large
 473   if (vertex->GetZRes() > 0.1) {
 474 #ifdef VERBOSE
 475     AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1",
 476                     vertex->GetZRes()));
 477 #endif
 478     return kFALSE;
 479   }
 480
 481   // Get the z coordiante
 482   vz = vertex->GetZ();
 483
 484   return kTRUE;
 485 }
 486
 487 //____________________________________________________________________
 488 void
 489 AliForwardMultiplicity::Print(Option_t*) const
 490 {
 491 }
 492
 493 //
 494 // EOF
 495 //