PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx

   1 //====================================================================
   2 #include "AliBasedNdetaTask.h"
   3 #include <TMath.h>
   4 #include <TH2D.h>
   5 #include <TH1D.h>
   6 #include <THStack.h>
   7 #include <TList.h>
   8 #include <TParameter.h>
   9 #include <AliAnalysisManager.h>
  10 #include <AliAODEvent.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODInputHandler.h>
  13 #include "AliForwardUtil.h"
  14 #include "AliAODForwardMult.h"
  15 #include <TFile.h>
  16 #include <TStyle.h>
  17
  18 //____________________________________________________________________
  19 AliBasedNdetaTask::AliBasedNdetaTask()
  20   : AliAnalysisTaskSE(),
  21     fSums(0),           // Container of sums
  22     fOutput(0),         // Container of output
  23     fVtxMin(0),         // Minimum v_z
  24     fVtxMax(0),         // Maximum v_z
  25     fTriggerMask(0),    // Trigger mask
  26     fRebin(0),          // Rebinning factor
  27     fCutEdges(false),
  28     fSymmetrice(true),
  29     fCorrEmpty(true),
  30     fUseROOTProj(false),
  31     fTriggerEff(1),
  32     fShapeCorr(0),
  33     fListOfCentralities(0),
  34     fSNNString(0),
  35     fSysString(0),
  36     fCent(0),
  37     fCentAxis(0),
  38     fNormalizationScheme(kFull),
  39     fSchemeString(0),
  40     fTriggerString(0),
  41     fFinalMCCorrFile(""),
  42     fIsESD(false)
  43 {
  44   //
  45   // Constructor
  46   //
  47   DGUARD(0,0,"Default construction of AliBasedNdetaTask");
  48 }
  49
  50 //____________________________________________________________________
  51 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  52   : AliAnalysisTaskSE(name),
  53     fSums(0),           // Container of sums
  54     fOutput(0),         // Container of output
  55     fVtxMin(-10),       // Minimum v_z
  56     fVtxMax(10),        // Maximum v_z
  57     fTriggerMask(AliAODForwardMult::kInel),
  58     fRebin(5),          // Rebinning factor
  59     fCutEdges(false),
  60     fSymmetrice(true),
  61     fCorrEmpty(true),
  62     fUseROOTProj(false),
  63     fTriggerEff(1),
  64     fShapeCorr(0),
  65     fListOfCentralities(0),
  66     fSNNString(0),
  67     fSysString(0),
  68     fCent(0),
  69     fCentAxis(0),
  70     fNormalizationScheme(kFull),
  71     fSchemeString(0),
  72     fTriggerString(0),
  73     fFinalMCCorrFile(""),
  74     fIsESD(false)
  75 {
  76   //
  77   // Constructor
  78   //
  79   DGUARD(0,0,"Named construction of AliBasedNdetaTask: %s", name);
  80   fListOfCentralities = new TObjArray(1);
  81
  82   // Set the normalisation scheme
  83   SetNormalizationScheme(kFull);
  84
  85   // Set the trigger mask
  86   SetTriggerMask(AliAODForwardMult::kInel);
  87
  88   // Output slot #1 writes into a TH1 container
  89   DefineOutput(1, TList::Class());
  90   DefineOutput(2, TList::Class());
  91 }
  92
  93 //____________________________________________________________________
  94 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
  95   : AliAnalysisTaskSE(o),
  96     fSums(o.fSums),             // TList* - Container of sums
  97     fOutput(o.fOutput),         // Container of output
  98     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
  99     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
 100     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
 101     fRebin(o.fRebin),           // Int_t - Rebinning factor
 102     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
 103     fSymmetrice(o.fSymmetrice),
 104     fCorrEmpty(o.fCorrEmpty),
 105     fUseROOTProj(o.fUseROOTProj),
 106     fTriggerEff(o.fTriggerEff),
 107     fShapeCorr(o.fShapeCorr),
 108     fListOfCentralities(o.fListOfCentralities),
 109     fSNNString(o.fSNNString),
 110     fSysString(o.fSysString),
 111     fCent(o.fCent),
 112     fCentAxis(o.fCentAxis),
 113     fNormalizationScheme(o.fNormalizationScheme),
 114     fSchemeString(o.fSchemeString),
 115     fTriggerString(o.fTriggerString),
 116     fFinalMCCorrFile(o.fFinalMCCorrFile),
 117     fIsESD(o.fIsESD)
 118 {
 119   DGUARD(0,0,"Copy construction of AliBasedNdetaTask");
 120 }
 121
 122 //____________________________________________________________________
 123 AliBasedNdetaTask::~AliBasedNdetaTask()
 124 {
 125   //
 126   // Destructor
 127   //
 128   DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
 129   if (fSums) {
 130     fSums->Delete();
 131     delete fSums;
 132     fSums = 0;
 133   }
 134   if (fOutput) {
 135     fOutput->Delete();
 136     delete fOutput;
 137     fOutput = 0;
 138   }
 139 }
 140
 141 //________________________________________________________________________
 142 void
 143 AliBasedNdetaTask::SetDebugLevel(Int_t lvl)
 144 {
 145   AliAnalysisTaskSE::SetDebugLevel(lvl);
 146   for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) {
 147     CentralityBin* bin =
 148       static_cast<CentralityBin*>(fListOfCentralities->At(i));
 149     bin->SetDebugLevel(lvl);
 150   }
 151 }
 152
 153 //________________________________________________________________________
 154 void
 155 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
 156 {
 157   DGUARD(fDebug,3,"Set centrality axis, %d bins", n);
 158   if (!fCentAxis) {
 159     fCentAxis = new TAxis();
 160     fCentAxis->SetName("centAxis");
 161     fCentAxis->SetTitle("Centrality [%]");
 162   }
 163   TArrayD dbins(n+1);
 164   for (UShort_t i = 0; i <= n; i++)
 165     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
 166   fCentAxis->Set(n, dbins.GetArray());
 167 }
 168
 169 //________________________________________________________________________
 170 void
 171 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
 172 {
 173   //
 174   // Add a centrality bin
 175   //
 176   // Parameters:
 177   //    low  Low cut
 178   //    high High cut
 179   //
 180   DGUARD(fDebug,3,"Add a centrality bin [%f,%f] @ %d", low, high, at);
 181   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
 182   bin->SetDebugLevel(fDebug);
 183   fListOfCentralities->AddAtAndExpand(bin, at);
 184 }
 185
 186 //________________________________________________________________________
 187 AliBasedNdetaTask::CentralityBin*
 188 AliBasedNdetaTask::MakeCentralityBin(const char* name,
 189                                      Short_t low, Short_t high) const
 190 {
 191   //
 192   // Make a centrality bin
 193   //
 194   // Parameters:
 195   //    name  Name used for histograms
 196   //    low   Low cut in percent
 197   //    high  High cut in percent
 198   //
 199   // Return:
 200   //    A newly created centrality bin
 201   //
 202   DGUARD(fDebug,3,"Make a centrality bin [%f,%f]: %s", low, high, name);
 203   return new CentralityBin(name, low, high);
 204 }
 205 //________________________________________________________________________
 206 void
 207 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
 208 {
 209   //
 210   // Set normalisation scheme
 211   //
 212   DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
 213   UShort_t    scheme = 0;
 214   TString     twhat(what);
 215   twhat.ToUpper();
 216   TObjString* opt;
 217   TIter       next(twhat.Tokenize(" ,|"));
 218   while ((opt = static_cast<TObjString*>(next()))) {
 219     TString s(opt->GetString());
 220     if      (s.IsNull()) continue;
 221     Bool_t add = true;
 222     switch (s[0]) {
 223     case '-': add = false; // Fall through
 224     case '+': s.Remove(0,1);  // Remove character
 225     }
 226     UShort_t bit = 0;
 227     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
 228     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
 229     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
 230     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
 231     else if (s.CompareTo("FULL")      == 0) bit = kFull;
 232     else if (s.CompareTo("NONE")      == 0) bit = kNone;
 233     else if (s.CompareTo("ZEROBIN")   == 0) bit = kZeroBin;
 234     else
 235       Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
 236     if (add) scheme |= bit;
 237     else     scheme ^= bit;
 238   }
 239   SetNormalizationScheme(scheme);
 240 }
 241 //________________________________________________________________________
 242 void
 243 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
 244 {
 245   DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
 246   fNormalizationScheme = scheme;
 247   TString tit = "";
 248   if (scheme == kFull) tit = "FULL";
 249   else {
 250     if (scheme & kEventLevel)        tit.Append("EVENT ");
 251     if (scheme & kShape)             tit.Append("SHAPE ");
 252     if (scheme & kBackground)        tit.Append("BACKGROUND ");
 253     if (scheme & kTriggerEfficiency) tit.Append("TRIGGER ");
 254     if (scheme & kZeroBin)           tit.Append("ZEROBIN ");
 255   }
 256   tit = tit.Strip(TString::kBoth);
 257   if (!fSchemeString) fSchemeString = new TNamed("scheme", "");
 258   fSchemeString->SetTitle(tit);
 259   fSchemeString->SetUniqueID(fNormalizationScheme);
 260 }
 261 //________________________________________________________________________
 262 void
 263 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 264 {
 265   //
 266   // Set the trigger maskl
 267   //
 268   // Parameters:
 269   //    mask Trigger mask
 270   //
 271   DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
 272   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
 273 }
 274 //________________________________________________________________________
 275 void
 276 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
 277 {
 278   DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
 279   fTriggerMask = mask;
 280   TString tit(AliAODForwardMult::GetTriggerString(mask));
 281   tit = tit.Strip(TString::kBoth);
 282   if (!fTriggerString) fTriggerString = new TNamed("trigger", "");
 283   fTriggerString->SetTitle(tit);
 284   fTriggerString->SetUniqueID(fTriggerMask);
 285 }
 286
 287 //________________________________________________________________________
 288 void
 289 AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
 290 {
 291   //
 292   // Set the shape correction (a.k.a., track correction) for selected
 293   // trigger(s)
 294   //
 295   // Parameters:
 296   //    h Correction
 297   //
 298   DGUARD(fDebug,3,"Set the shape correction: %p", c);
 299   if (!c) return;
 300   fShapeCorr = static_cast<TH1*>(c->Clone());
 301   fShapeCorr->SetDirectory(0);
 302 }
 303
 304 //________________________________________________________________________
 305 void
 306 AliBasedNdetaTask::UserCreateOutputObjects()
 307 {
 308   //
 309   // Create output objects.
 310   //
 311   // This is called once per slave process
 312   //
 313   DGUARD(fDebug,1,"Create user ouput object");
 314   fSums = new TList;
 315   fSums->SetName(Form("%s_sums", GetName()));
 316   fSums->SetOwner();
 317
 318   // Automatically add 'all' centrality bin if nothing has been defined.
 319   AddCentralityBin(0, 0, 0);
 320   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
 321     const TArrayD* bins = fCentAxis->GetXbins();
 322     Int_t          nbin = fCentAxis->GetNbins();
 323     for (Int_t i = 0; i < nbin; i++)
 324       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
 325   }
 326   if (fCentAxis) fSums->Add(fCentAxis);
 327
 328
 329   // Centrality histogram
 330   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
 331   fCent->SetDirectory(0);
 332   fCent->SetXTitle(0);
 333   fSums->Add(fCent);
 334
 335   // Loop over centrality bins
 336   TIter next(fListOfCentralities);
 337   CentralityBin* bin = 0;
 338   while ((bin = static_cast<CentralityBin*>(next())))
 339     bin->CreateOutputObjects(fSums);
 340
 341   // Check that we have an AOD input handler
 342   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 343   AliAODInputHandler* ah =
 344     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
 345   if (!ah) {
 346     fIsESD = true;
 347     const char* depName = GetTitle();
 348     AliAnalysisTask* dep = am->GetTask(depName); // "Central"
 349     if (!dep) {
 350       AliErrorF("The multiplicity task %s wasn't added to the train", depName);
 351       fIsESD = false;
 352     }
 353   }
 354   if (!fIsESD && !ah) AliFatal("No AOD input handler set in analysis manager");
 355
 356   // Post data for ALL output slots >0 here, to get at least an empty histogram
 357   PostData(1, fSums);
 358 }
 359 //____________________________________________________________________
 360 AliAODEvent*
 361 AliBasedNdetaTask::GetAODEvent(bool isESD)
 362 {
 363   // Get the AOD event
 364   //
 365   // Parameters:
 366   //   isESD     Whether we have ESD input or not
 367   //
 368   // Return:
 369   //   Pointer to AOD event or null
 370   DGUARD(fDebug,1,"Get the AOD event (%s analysis)", isESD ? "ESD" : "AOD");
 371   AliAODEvent* aod = 0;
 372
 373   if (!isESD) aod = dynamic_cast<AliAODEvent*>(InputEvent());
 374   else        aod = AODEvent();
 375   if (!aod) AliError("Cannot get the AOD event");
 376
 377   return aod;
 378 }
 379
 380 //____________________________________________________________________
 381 void
 382 AliBasedNdetaTask::UserExec(Option_t *)
 383 {
 384   //
 385   // Process a single event
 386   //
 387   // Parameters:
 388   //    option Not used
 389   //
 390   // Main loop
 391   DGUARD(fDebug,1,"Analyse the AOD event");
 392   AliAODEvent* aod = GetAODEvent(fIsESD);
 393   if (!aod) return;
 394
 395
 396   TObject* obj = aod->FindListObject("Forward");
 397   if (!obj) {
 398     AliWarning("No forward object found");
 399     return;
 400   }
 401   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 402
 403   // Fill centrality histogram
 404   Float_t cent = forward->GetCentrality();
 405   fCent->Fill(cent);
 406
 407   // Get the histogram(s)
 408   TH2D* data   = GetHistogram(aod, false);
 409   TH2D* dataMC = GetHistogram(aod, true);
 410
 411   Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
 412                    !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
 413
 414
 415   // Loop over centrality bins
 416   CentralityBin* allBin =
 417     static_cast<CentralityBin*>(fListOfCentralities->At(0));
 418   allBin->ProcessEvent(forward, fTriggerMask, isZero,
 419                        fVtxMin, fVtxMax, data, dataMC);
 420
 421   // Find this centrality bin
 422   if (fCentAxis && fCentAxis->GetNbins() > 0) {
 423     Int_t          icent   = fCentAxis->FindBin(cent);
 424     CentralityBin* thisBin = 0;
 425     if (icent >= 1 && icent <= fCentAxis->GetNbins())
 426       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
 427     if (thisBin)
 428       thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin, fVtxMax,
 429                             data, dataMC);
 430   }
 431
 432   // Here, we get the update
 433   if (!fSNNString) {
 434     UShort_t sNN = forward->GetSNN();
 435     fSNNString = new TNamed("sNN", "");
 436     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
 437     fSNNString->SetUniqueID(sNN);
 438     fSums->Add(fSNNString);
 439
 440     UShort_t sys = forward->GetSystem();
 441     fSysString = new TNamed("sys", "");
 442     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
 443     fSysString->SetUniqueID(sys);
 444     fSums->Add(fSysString);
 445
 446     fSums->Add(fSchemeString);
 447     fSums->Add(fTriggerString);
 448
 449     // Print();
 450   }
 451
 452   PostData(1, fSums);
 453 }
 454
 455 //________________________________________________________________________
 456 void
 457 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 458                                           const char* title, const char* ytitle)
 459 {
 460   //
 461   // Set histogram graphical options, etc.
 462   //
 463   // Parameters:
 464   //    h       Histogram to modify
 465   //    colour  Marker color
 466   //    marker  Marker style
 467   //    title   Title of histogram
 468   //    ytitle  Title on y-axis.
 469   //
 470   h->SetTitle(title);
 471   h->SetMarkerColor(colour);
 472   h->SetMarkerStyle(marker);
 473   h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
 474   h->SetFillStyle(0);
 475   h->SetYTitle(ytitle);
 476   h->GetXaxis()->SetTitleFont(132);
 477   h->GetXaxis()->SetLabelFont(132);
 478   h->GetXaxis()->SetNdivisions(10);
 479   h->GetYaxis()->SetTitleFont(132);
 480   h->GetYaxis()->SetLabelFont(132);
 481   h->GetYaxis()->SetNdivisions(10);
 482   h->GetYaxis()->SetDecimals();
 483   h->SetStats(0);
 484 }
 485
 486 //________________________________________________________________________
 487 void
 488 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm)
 489 {
 490   // Normalize to the acceptance -
 491   // dndeta->Divide(accNorm);
 492   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 493     Double_t a = norm->GetBinContent(i);
 494     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
 495       if (a <= 0) {
 496         copy->SetBinContent(i,j,0);
 497         copy->SetBinError(i,j,0);
 498         continue;
 499       }
 500       Double_t c = copy->GetBinContent(i, j);
 501       Double_t e = copy->GetBinError(i, j);
 502       copy->SetBinContent(i, j, c / a);
 503       copy->SetBinError(i, j, e / a);
 504     }
 505   }
 506 }
 507
 508 //________________________________________________________________________
 509 TH1D*
 510 AliBasedNdetaTask::ProjectX(const TH2D* h,
 511                             const char* name,
 512                             Int_t firstbin,
 513                             Int_t lastbin,
 514                             bool  useRoot,
 515                             bool  corr,
 516                             bool  error)
 517 {
 518   //
 519   // Project onto the X axis
 520   //
 521   // Parameters:
 522   //    h         2D histogram
 523   //    name      New name
 524   //    firstbin  First bin to use
 525   //    lastbin   Last bin to use
 526   //    error     Whether to calculate errors
 527   //
 528   // Return:
 529   //    Newly created histogram or null
 530   //
 531   if (!h) return 0;
 532   if (useRoot)
 533     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 534
 535   TAxis* xaxis = h->GetXaxis();
 536   TAxis* yaxis = h->GetYaxis();
 537   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 538                           xaxis->GetXmin(), xaxis->GetXmax());
 539   static_cast<const TAttLine*>(h)->Copy(*ret);
 540   static_cast<const TAttFill*>(h)->Copy(*ret);
 541   static_cast<const TAttMarker*>(h)->Copy(*ret);
 542   ret->GetXaxis()->ImportAttributes(xaxis);
 543
 544   Int_t first = firstbin;
 545   Int_t last  = lastbin;
 546   if (first < 0)                         first = 0;
 547   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
 548   if      (last  < 0)                    last  = yaxis->GetNbins();
 549   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
 550   if (last-first < 0) {
 551     AliWarningGeneral("AliBasedNdetaTask",
 552                       Form("Nothing to project [%d,%d]", first, last));
 553     return 0;
 554
 555   }
 556
 557   // Loop over X bins
 558   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 559   Int_t ybins = (last-first+1);
 560   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 561     Double_t content = 0;
 562     Double_t error2  = 0;
 563     Int_t    nbins   = 0;
 564
 565
 566     for (Int_t ybin = first; ybin <= last; ybin++) {
 567       Double_t c1 = h->GetCellContent(xbin, ybin);
 568       Double_t e1 = h->GetCellError(xbin, ybin);
 569
 570       // Ignore empty bins
 571       if (c1 < 1e-12) continue;
 572       if (e1 < 1e-12) {
 573         if (error) continue;
 574         e1 = 1;
 575       }
 576
 577       content    += c1;
 578       error2     += e1*e1;
 579       nbins++;
 580     } // for (ybin)
 581     if(content > 0 && nbins > 0) {
 582       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 583 #if 0
 584       AliWarningGeneral(ret->GetName(),
 585                         Form("factor @ %d is %d/%d -> %f",
 586                              xbin, ybins, nbins, factor));
 587 #endif
 588       if (error) {
 589         // Calculate weighted average
 590         ret->SetBinContent(xbin, content * factor);
 591         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 592       }
 593       else
 594         ret->SetBinContent(xbin, factor * content);
 595     }
 596   } // for (xbin)
 597
 598   return ret;
 599 }
 600
 601 //________________________________________________________________________
 602 void
 603 AliBasedNdetaTask::Terminate(Option_t *)
 604 {
 605   //
 606   // Called at end of event processing..
 607   //
 608   // This is called once in the master
 609   //
 610   // Parameters:
 611   //    option Not used
 612   //
 613   // Draw result to screen, or perform fitting, normalizations Called
 614   // once at the end of the query
 615   DGUARD(fDebug,1,"Process final merged results");
 616
 617   fSums = dynamic_cast<TList*> (GetOutputData(1));
 618   if(!fSums) {
 619     AliError("Could not retrieve TList fSums");
 620     return;
 621   }
 622
 623   fOutput = new TList;
 624   fOutput->SetName(Form("%s_result", GetName()));
 625   fOutput->SetOwner();
 626
 627   fSNNString     = static_cast<TNamed*>(fSums->FindObject("sNN"));
 628   fSysString     = static_cast<TNamed*>(fSums->FindObject("sys"));
 629   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
 630   fSchemeString  = static_cast<TNamed*>(fSums->FindObject("scheme"));
 631   fTriggerString = static_cast<TNamed*>(fSums->FindObject("trigger"));
 632
 633   if(fSysString && fSNNString &&
 634      fSysString->GetUniqueID() == AliForwardUtil::kPP)
 635     LoadNormalizationData(fSysString->GetUniqueID(),
 636                           fSNNString->GetUniqueID());
 637   // Print before we loop
 638   Print();
 639
 640   // Loop over centrality bins
 641   TIter next(fListOfCentralities);
 642   CentralityBin* bin = 0;
 643   gStyle->SetPalette(1);
 644   THStack* dndetaStack        = new THStack("dndeta", "dN/d#eta");
 645   THStack* dndetaStackRebin   = new THStack(Form("dndeta_rebin%02d", fRebin),
 646                                             "dN_{ch}/d#eta");
 647   THStack* dndetaMCStack      = new THStack("dndetaMC", "dN_{ch}/d#eta");
 648   THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
 649                                             "dN_{ch}/d#eta");
 650
 651   TList* mclist = 0;
 652   TList* truthlist = 0;
 653
 654   if (fFinalMCCorrFile.Contains(".root")) {
 655     TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
 656     if(ftest) {
 657       mclist    = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
 658       truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
 659     }
 660     else
 661       AliWarning("MC analysis file invalid - no final MC correction possible");
 662   }
 663   Int_t style = GetMarker();
 664   Int_t color = GetColor();
 665
 666   AliInfo(Form("Marker style=%d, color=%d", style, color));
 667   while ((bin = static_cast<CentralityBin*>(next()))) {
 668
 669     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
 670              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
 671              fTriggerMask, style, color, mclist, truthlist);
 672     if (fCentAxis && bin->IsAllBin()) continue;
 673     TH1* dndeta      = bin->GetResult(0, false, "");
 674     TH1* dndetaSym   = bin->GetResult(0, true,  "");
 675     TH1* dndetaMC    = bin->GetResult(0, false, "MC");
 676     TH1* dndetaMCSym = bin->GetResult(0, true,  "MC");
 677     if (dndeta)      dndetaStack->Add(dndeta);
 678     if (dndetaSym)   dndetaStack->Add(dndetaSym);
 679     if (dndetaMC)    dndetaMCStack->Add(dndetaMC);
 680     if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
 681     if (fRebin > 1) {
 682       dndeta      = bin->GetResult(fRebin, false, "");
 683       dndetaSym   = bin->GetResult(fRebin, true,  "");
 684       dndetaMC    = bin->GetResult(fRebin, false, "MC");
 685       dndetaMCSym = bin->GetResult(fRebin, true,  "MC");
 686       if (dndeta)      dndetaStackRebin->Add(dndeta);
 687       if (dndetaSym)   dndetaStackRebin->Add(dndetaSym);
 688       if (dndetaMC)    dndetaMCStackRebin->Add(dndetaMC);
 689       if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
 690     }
 691   }
 692   // Output the stack
 693   fOutput->Add(dndetaStack);
 694
 695   // If available output rebinned stack
 696   if (!dndetaStackRebin->GetHists() ||
 697       dndetaStackRebin->GetHists()->GetEntries() <= 0) {
 698     AliWarning("No rebinned histograms found");
 699     delete dndetaStackRebin;
 700     dndetaStackRebin = 0;
 701   }
 702   if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
 703
 704   // If available, output track-ref stack
 705   if (!dndetaMCStack->GetHists() ||
 706       dndetaMCStack->GetHists()->GetEntries() <= 0) {
 707     AliWarning("No MC histograms found");
 708     delete dndetaMCStack;
 709     dndetaMCStack = 0;
 710   }
 711   if (dndetaMCStack) fOutput->Add(dndetaMCStack);
 712
 713   // If available, output rebinned track-ref stack
 714   if (!dndetaMCStackRebin->GetHists() ||
 715       dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
 716     AliWarning("No rebinned MC histograms found");
 717     delete dndetaMCStackRebin;
 718     dndetaMCStackRebin = 0;
 719   }
 720   if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
 721
 722   // Output collision energy string
 723   if (fSNNString) fOutput->Add(fSNNString->Clone());
 724
 725   // Output collision system string
 726   if (fSysString) fOutput->Add(fSysString->Clone());
 727
 728   // Output centrality axis
 729   if (fCentAxis) fOutput->Add(fCentAxis);
 730
 731   // Output trigger string
 732   if (fTriggerString) fOutput->Add(fTriggerString->Clone());
 733
 734   // Normalization string
 735   if (fSchemeString) fOutput->Add(fSchemeString->Clone());
 736
 737   // Output vertex axis
 738   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 739   vtxAxis->SetName("vtxAxis");
 740   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 741   fOutput->Add(vtxAxis);
 742
 743   // Output trigger efficiency and shape correction
 744   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
 745   if (fShapeCorr) fOutput->Add(fShapeCorr);
 746
 747   TNamed* options = new TNamed("options","");
 748   TString str;
 749   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
 750   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
 751   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
 752   options->SetTitle(str);
 753   fOutput->Add(options);
 754
 755   PostData(2, fOutput);
 756 }
 757 //________________________________________________________________________
 758 void
 759 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 760 {
 761   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
 762   DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
 763          sys, energy);
 764   TString type("pp");
 765   TString snn("900");
 766   if(energy == 7000) snn.Form("7000");
 767   if(energy == 2750) snn.Form("2750");
 768
 769   if(fShapeCorr &&  (fTriggerEff != 1)) {
 770     AliInfo("Objects already set for normalization - no action taken");
 771     return;
 772   }
 773
 774   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
 775                                 "Normalization/normalizationHists_%s_%s.root",
 776                                 type.Data(),snn.Data()));
 777   if(!fin) {
 778     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
 779     return;
 780   }
 781
 782   // Shape correction
 783   if ((fNormalizationScheme & kShape) && !fShapeCorr) {
 784     TString trigName("All");
 785     if (fTriggerMask == AliAODForwardMult::kInel ||
 786         fTriggerMask == AliAODForwardMult::kNClusterGt0)
 787       trigName = "Inel";
 788     else if (fTriggerMask == AliAODForwardMult::kNSD)
 789       trigName = "NSD";
 790     else if (fTriggerMask == AliAODForwardMult::kInelGt0)
 791       trigName = "InelGt0";
 792     else {
 793       AliWarning(Form("Normalization for trigger %s not known, using all",
 794                       AliAODForwardMult::GetTriggerString(fTriggerMask)));
 795     }
 796
 797     TString shapeCorName(Form("h%sNormalization", trigName.Data()));
 798     TH2F*   shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
 799     if (shapeCor) SetShapeCorrection(shapeCor);
 800     else {
 801       AliWarning(Form("No shape correction found for %s", trigName.Data()));
 802     }
 803   }
 804
 805   // Trigger efficiency
 806   TString effName(Form("%sTriggerEff",
 807                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 808                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 809                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
 810                        "inelgt0" : "all"));
 811   TParameter<float>* eff = 0;
 812   if (fNormalizationScheme & kTriggerEfficiency)
 813     eff = static_cast<TParameter<float>*>(fin->Get(effName));
 814   Double_t trigEff = eff ? eff->GetVal() : 1;
 815   if (fTriggerEff != 1) SetTriggerEff(trigEff);
 816   if (fTriggerEff < 0)  fTriggerEff = 1;
 817
 818   // TEMPORARY FIX
 819   // Rescale the shape correction by the trigger efficiency
 820   if (fShapeCorr) {
 821     AliWarning(Form("Rescaling shape correction by trigger efficiency: "
 822                     "1/E_X=1/%f", trigEff));
 823     fShapeCorr->Scale(1. / trigEff);
 824   }
 825
 826   // Print - out
 827   if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
 828 }
 829
 830
 831 //________________________________________________________________________
 832 void
 833 AliBasedNdetaTask::Print(Option_t*) const
 834 {
 835   //
 836   // Print information
 837   //
 838   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
 839             << std::boolalpha
 840             << " Trigger:              " << (fTriggerString ?
 841                                              fTriggerString->GetTitle() :
 842                                              "none") << "\n"
 843             << " Vertex range:         [" << fVtxMin << ":" << fVtxMax << "]\n"
 844             << " Rebin factor:         " << fRebin << "\n"
 845             << " Cut edges:            " << fCutEdges << "\n"
 846             << " Symmertrice:          " << fSymmetrice << "\n"
 847             << " Use TH2::ProjectionX: " << fUseROOTProj << "\n"
 848             << " Correct for empty:    " << fCorrEmpty << "\n"
 849             << " Normalization scheme: " << (fSchemeString ?
 850                                              fSchemeString->GetTitle() :
 851                                              "none") <<"\n"
 852             << " Trigger efficiency:   " << fTriggerEff << "\n"
 853             << " Shape correction:     " << (fShapeCorr ?
 854                                            fShapeCorr->GetName() :
 855                                            "none") << "\n"
 856             << " sqrt(s_NN):           " << (fSNNString ?
 857                                              fSNNString->GetTitle() :
 858                                            "unknown") << "\n"
 859             << " Collision system:     " << (fSysString ?
 860                                              fSysString->GetTitle() :
 861                                            "unknown") << "\n"
 862             << " Centrality bins:      " << (fCentAxis ? "" : "none");
 863   if (fCentAxis) {
 864     Int_t           nBins = fCentAxis->GetNbins();
 865     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
 866     for (Int_t i = 0; i <= nBins; i++)
 867       std::cout << (i==0 ? "" : "-") << bins[i];
 868   }
 869   std::cout << std::noboolalpha << std::endl;
 870
 871 }
 872
 873 //________________________________________________________________________
 874 TH1D*
 875 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
 876 {
 877   //
 878   // Make a copy of the input histogram and rebin that histogram
 879   //
 880   // Parameters:
 881   //    h  Histogram to rebin
 882   //
 883   // Return:
 884   //    New (rebinned) histogram
 885   //
 886   if (rebin <= 1) return 0;
 887
 888   Int_t nBins = h->GetNbinsX();
 889   if(nBins % rebin != 0) {
 890     AliWarningGeneral("AliBasedNdetaTask",
 891                       Form("Rebin factor %d is not a devisor of current number "
 892                            "of bins %d in the histogram %s",
 893                            rebin, nBins, h->GetName()));
 894     return 0;
 895   }
 896
 897   // Make a copy
 898   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
 899                                                h->GetName(), rebin)));
 900   tmp->Rebin(rebin);
 901   tmp->SetDirectory(0);
 902
 903   // The new number of bins
 904   Int_t nBinsNew = nBins / rebin;
 905   for(Int_t i = 1;i<= nBinsNew; i++) {
 906     Double_t content = 0;
 907     Double_t sumw    = 0;
 908     Double_t wsum    = 0;
 909     Int_t    nbins   = 0;
 910     for(Int_t j = 1; j<=rebin;j++) {
 911       Int_t    bin = (i-1)*rebin + j;
 912       Double_t c   =  h->GetBinContent(bin);
 913       if (c <= 0) continue;
 914
 915       if (cutEdges) {
 916         if (h->GetBinContent(bin+1)<=0 ||
 917             h->GetBinContent(bin-1)<=0) {
 918 #if 0
 919           AliWarningGeneral("AliBasedNdetaTask",
 920                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
 921                                  bin, c, h->GetName(),
 922                                  bin+1, h->GetBinContent(bin+1),
 923                                  bin-1, h->GetBinContent(bin-1)));
 924 #endif
 925           continue;
 926         }
 927       }
 928       Double_t e =  h->GetBinError(bin);
 929       Double_t w =  1 / (e*e); // 1/c/c
 930       content    += c;
 931       sumw       += w;
 932       wsum       += w * c;
 933       nbins++;
 934     }
 935
 936     if(content > 0 && nbins > 0) {
 937       tmp->SetBinContent(i, wsum / sumw);
 938       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 939     }
 940   }
 941
 942   return tmp;
 943 }
 944
 945 //__________________________________________________________________
 946 TH1*
 947 AliBasedNdetaTask::Symmetrice(const TH1* h)
 948 {
 949   //
 950   // Make an extension of @a h to make it symmetric about 0
 951   //
 952   // Parameters:
 953   //    h Histogram to symmertrice
 954   //
 955   // Return:
 956   //    Symmetric extension of @a h
 957   //
 958   Int_t nBins = h->GetNbinsX();
 959   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 960   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 961   s->Reset();
 962   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 963   s->SetMarkerColor(h->GetMarkerColor());
 964   s->SetMarkerSize(h->GetMarkerSize());
 965   s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
 966   s->SetFillColor(h->GetFillColor());
 967   s->SetFillStyle(h->GetFillStyle());
 968   s->SetDirectory(0);
 969
 970   // Find the first and last bin with data
 971   Int_t first = nBins+1;
 972   Int_t last  = 0;
 973   for (Int_t i = 1; i <= nBins; i++) {
 974     if (h->GetBinContent(i) <= 0) continue;
 975     first = TMath::Min(first, i);
 976     last  = TMath::Max(last,  i);
 977   }
 978
 979   Double_t xfirst = h->GetBinCenter(first-1);
 980   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 981   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 982   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 983     s->SetBinContent(j, h->GetBinContent(i));
 984     s->SetBinError(j, h->GetBinError(i));
 985   }
 986   // Fill in overlap bin
 987   s->SetBinContent(l2+1, h->GetBinContent(first));
 988   s->SetBinError(l2+1, h->GetBinError(first));
 989   return s;
 990 }
 991
 992 //__________________________________________________________________
 993 Int_t
 994 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
 995 {
 996   Int_t  base   = bits & (0xFE);
 997   Bool_t hollow = bits & kHollow;
 998   switch (base) {
 999   case kCircle:       return (hollow ? 24 : 20);
1000   case kSquare:       return (hollow ? 25 : 21);
1001   case kUpTriangle:   return (hollow ? 26 : 22);
1002   case kDownTriangle: return (hollow ? 32 : 23);
1003   case kDiamond:      return (hollow ? 27 : 33);
1004   case kCross:        return (hollow ? 28 : 34);
1005   case kStar:         return (hollow ? 30 : 29);
1006   }
1007   return 1;
1008 }
1009 //__________________________________________________________________
1010 UShort_t
1011 AliBasedNdetaTask::GetMarkerBits(Int_t style)
1012 {
1013   UShort_t bits = 0;
1014   switch (style) {
1015   case 24: case 25: case 26: case 27: case 28: case 30: case 32:
1016     bits |= kHollow; break;
1017   }
1018   switch (style) {
1019   case 20: case 24: bits |= kCircle;       break;
1020   case 21: case 25: bits |= kSquare;       break;
1021   case 22: case 26: bits |= kUpTriangle;   break;
1022   case 23: case 32: bits |= kDownTriangle; break;
1023   case 27: case 33: bits |= kDiamond;      break;
1024   case 28: case 34: bits |= kCross;        break;
1025   case 29: case 30: bits |= kStar;         break;
1026   }
1027   return bits;
1028 }
1029 //__________________________________________________________________
1030 Int_t
1031 AliBasedNdetaTask::FlipHollowStyle(Int_t style)
1032 {
1033   UShort_t bits = GetMarkerBits(style);
1034   Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
1035   return ret;
1036 }
1037
1038 //====================================================================
1039 void
1040 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1041 {
1042   DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1043   TString n(GetHistName(0));
1044   TString n0(GetHistName(1));
1045   const char* postfix = GetTitle();
1046
1047   fSum = static_cast<TH2D*>(data->Clone(n));
1048   if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1049   fSum->SetDirectory(0);
1050   fSum->SetMarkerColor(col);
1051   fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1052   fSum->Reset();
1053   list->Add(fSum);
1054
1055   fSum0 = static_cast<TH2D*>(data->Clone(n0));
1056   if (postfix)
1057     fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1058   else
1059     fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1060   fSum0->SetDirectory(0);
1061   fSum0->SetMarkerColor(col);
1062   fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1063   fSum0->Reset();
1064   list->Add(fSum0);
1065
1066   fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1067   fEvents->SetDirectory(0);
1068   fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1069   fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1070   list->Add(fEvents);
1071 }
1072
1073 //____________________________________________________________________
1074 TString
1075 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1076 {
1077   TString n(GetName());
1078   if      (what == 1) n.Append("0");
1079   else if (what == 2) n.Append("Events");
1080   const char* postfix = GetTitle();
1081   if (postfix && postfix[0] != '\0')  n.Append(postfix);
1082   return n;
1083 }
1084
1085 //____________________________________________________________________
1086 void
1087 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1088 {
1089   DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1090   if (isZero) fSum0->Add(data);
1091   else        fSum->Add(data);
1092   fEvents->Fill(isZero ? 1 : 0);
1093 }
1094
1095 //____________________________________________________________________
1096 TH2D*
1097 AliBasedNdetaTask::Sum::GetSum(const TList* input,
1098                                TList*       output,
1099                                Double_t&    ntotal,
1100                                Double_t     epsilon0,
1101                                Double_t     epsilon,
1102                                Int_t        marker,
1103                                Bool_t       rootProj,
1104                                Bool_t       corrEmpty) const
1105 {
1106   DGUARD(fDebug,2,"Get sums from input list %s", input->GetName());
1107   TH2D* sum      = static_cast<TH2D*>(input->FindObject(GetHistName(0)));
1108   TH2D* sum0     = static_cast<TH2D*>(input->FindObject(GetHistName(1)));
1109   TH1I* events   = static_cast<TH1I*>(input->FindObject(GetHistName(2)));
1110   if (!sum || !sum0 || !events) {
1111     AliWarning(Form("Failed to find one or more histograms: "
1112                     "%s (%p) %s (%p) %s (%p)",
1113                     GetHistName(0).Data(), sum,
1114                     GetHistName(1).Data(), sum0,
1115                     GetHistName(2).Data(), events));
1116     return 0;
1117   }
1118
1119   TH2D* ret      = static_cast<TH2D*>(sum->Clone(sum->GetName()));
1120   ret->SetDirectory(0);
1121   ret->Reset();
1122   Int_t n        = Int_t(events->GetBinContent(1));
1123   Int_t n0       = Int_t(events->GetBinContent(2));
1124
1125   AliInfo(Form("Adding histograms %s and %s with weights %f and %f resp.",
1126                sum0->GetName(), sum->GetName(), 1./epsilon, 1./epsilon0));
1127   // Generate merged histogram
1128   ret->Add(sum0, sum, 1. / epsilon0, 1. / epsilon);
1129   ntotal = n / epsilon + n0 / epsilon0;
1130
1131   TList* out = new TList;
1132   out->SetOwner();
1133   const char* postfix = GetTitle();
1134   if (!postfix) postfix = "";
1135   out->SetName(Form("partial%s", postfix));
1136   output->Add(out);
1137
1138   // Now make copies, normalize them, and store in output list
1139   TH2D* sumCopy  = static_cast<TH2D*>(sum->Clone("sum"));
1140   TH2D* sum0Copy = static_cast<TH2D*>(sum0->Clone("sum0"));
1141   TH2D* retCopy  = static_cast<TH2D*>(ret->Clone("sumAll"));
1142   sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1143   sumCopy->SetDirectory(0);
1144   sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1145   sum0Copy->SetDirectory(0);
1146   retCopy->SetMarkerStyle(marker);
1147   retCopy->SetDirectory(0);
1148
1149   TH1D* norm    = ProjectX(sum,  "norm",    0, 0, rootProj, corrEmpty, false);
1150   TH1D* norm0   = ProjectX(sum0, "norm0",   0, 0, rootProj, corrEmpty, false);
1151   TH1D* normAll = ProjectX(ret,  "normAll", 0, 0, rootProj, corrEmpty, false);
1152   norm->SetDirectory(0);
1153   norm0->SetDirectory(0);
1154   normAll->SetDirectory(0);
1155
1156   ScaleToCoverage(sumCopy, norm);
1157   ScaleToCoverage(sum0Copy, norm0);
1158   ScaleToCoverage(retCopy, normAll);
1159
1160   Int_t nY = sum->GetNbinsY();
1161   TH1D* sumCopyPx  = ProjectX(sumCopy,  "average",    1, nY,rootProj,corrEmpty);
1162   TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0",   1, nY,rootProj,corrEmpty);
1163   TH1D* retCopyPx  = ProjectX(retCopy,  "averageAll", 1, nY,rootProj,corrEmpty);
1164   sumCopyPx->SetDirectory(0);
1165   sum0CopyPx->SetDirectory(0);
1166   retCopyPx->SetDirectory(0);
1167
1168   // Scale our 1D histograms
1169   sumCopyPx->Scale(1., "width");
1170   sum0CopyPx->Scale(1., "width");
1171   retCopyPx->Scale(1., "width");
1172
1173   AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1174                sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1175
1176   // Scale the normalization - they should be 1 at the maximum
1177   norm->Scale(n > 0   ? 1. / n  : 1);
1178   norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1179   normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1180
1181   out->Add(sumCopy);
1182   out->Add(sum0Copy);
1183   out->Add(retCopy);
1184   out->Add(sumCopyPx);
1185   out->Add(sum0CopyPx);
1186   out->Add(retCopyPx);
1187   out->Add(norm);
1188   out->Add(norm0);
1189   out->Add(normAll);
1190
1191   AliInfo(Form("Returning  (1/%f * %s + 1/%f * %s), "
1192                "1/%f * %d + 1/%f * %d = %d",
1193                epsilon0, sum0->GetName(), epsilon, sum->GetName(),
1194                epsilon0, n0, epsilon, n, int(ntotal)));
1195 #if 0
1196   for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1197     Double_t nc  = sum->GetBinContent(i, 0);
1198     Double_t nc0 = sum0->GetBinContent(i, 0);
1199     ret->SetBinContent(i, 0, nc + nc0); // Just count events
1200   }
1201 #endif
1202
1203   return ret;
1204 }
1205
1206 //====================================================================
1207 AliBasedNdetaTask::CentralityBin::CentralityBin()
1208   : TNamed("", ""),
1209     fSums(0),
1210     fOutput(0),
1211     fSum(0),
1212     fSumMC(0),
1213     fTriggers(0),
1214     fLow(0),
1215     fHigh(0),
1216     fDoFinalMCCorrection(false),
1217     fDebug(0)
1218 {
1219   //
1220   // Constructor
1221   //
1222   DGUARD(0,0,"Centrality bin default construction");
1223 }
1224 //____________________________________________________________________
1225 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
1226                                                 Short_t low, Short_t high)
1227   : TNamed(name, ""),
1228     fSums(0),
1229     fOutput(0),
1230     fSum(0),
1231     fSumMC(0),
1232     fTriggers(0),
1233     fLow(low),
1234     fHigh(high),
1235     fDoFinalMCCorrection(false),
1236     fDebug(0)
1237 {
1238   //
1239   // Constructor
1240   //
1241   // Parameters:
1242   //    name Name used for histograms (e.g., Forward)
1243   //    low  Lower centrality cut in percent
1244   //    high Upper centrality cut in percent
1245   //
1246   DGUARD(0,0,"Named Centrality bin construction: %s [%d,%d]", name);
1247   if (low <= 0 && high <= 0) {
1248     fLow  = 0;
1249     fHigh = 0;
1250     SetTitle("All centralities");
1251   }
1252   else {
1253     fLow  = low;
1254     fHigh = high;
1255     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1256   }
1257 }
1258 //____________________________________________________________________
1259 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1260   : TNamed(o),
1261     fSums(o.fSums),
1262     fOutput(o.fOutput),
1263     fSum(o.fSum),
1264     fSumMC(o.fSumMC),
1265     fTriggers(o.fTriggers),
1266     fLow(o.fLow),
1267     fHigh(o.fHigh),
1268     fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1269     fDebug(o.fDebug)
1270 {
1271   //
1272   // Copy constructor
1273   //
1274   // Parameters:
1275   //    other Object to copy from
1276   //
1277   DGUARD(0,0,"Copy Centrality bin construction");
1278 }
1279 //____________________________________________________________________
1280 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1281 {
1282   //
1283   // Destructor
1284   //
1285   DGUARD(fDebug,3,"Centrality bin desctruction");
1286   if (fSums) fSums->Delete();
1287   if (fOutput) fOutput->Delete();
1288 }
1289
1290 //____________________________________________________________________
1291 AliBasedNdetaTask::CentralityBin&
1292 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1293 {
1294   //
1295   // Assignment operator
1296   //
1297   // Parameters:
1298   //    other Object to assign from
1299   //
1300   // Return:
1301   //    Reference to this
1302   //
1303   DGUARD(fDebug,3,"Centrality bin assignment");
1304   if (&o == this) return *this;
1305   SetName(o.GetName());
1306   SetTitle(o.GetTitle());
1307   fSums      = o.fSums;
1308   fOutput    = o.fOutput;
1309   fSum       = o.fSum;
1310   fSumMC     = o.fSumMC;
1311   fTriggers  = o.fTriggers;
1312   fLow       = o.fLow;
1313   fHigh      = o.fHigh;
1314   fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1315
1316   return *this;
1317 }
1318 //____________________________________________________________________
1319 Int_t
1320 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const
1321 {
1322   if (IsAllBin()) return fallback;
1323   Float_t  fc       = (fLow+double(fHigh-fLow)/2) / 100;
1324   Int_t    nCol     = gStyle->GetNumberOfColors();
1325   Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
1326   Int_t    col      = gStyle->GetColorPalette(icol);
1327   return col;
1328 }
1329 //____________________________________________________________________
1330 const char*
1331 AliBasedNdetaTask::CentralityBin::GetListName() const
1332 {
1333   //
1334   // Get the list name
1335   //
1336   // Return:
1337   //    List Name
1338   //
1339   if (IsAllBin()) return "all";
1340   return Form("cent%03d_%03d", fLow, fHigh);
1341 }
1342 //____________________________________________________________________
1343 void
1344 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
1345 {
1346   //
1347   // Create output objects
1348   //
1349   // Parameters:
1350   //    dir   Parent list
1351   //
1352   DGUARD(fDebug,1,"Create centrality bin output objects");
1353   fSums = new TList;
1354   fSums->SetName(GetListName());
1355   fSums->SetOwner();
1356   dir->Add(fSums);
1357
1358   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
1359   fTriggers->SetDirectory(0);
1360   fSums->Add(fTriggers);
1361 }
1362 //____________________________________________________________________
1363 void
1364 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1365 {
1366   fDebug = lvl;
1367   if (fSum) fSum->fDebug = lvl;
1368   if (fSumMC) fSumMC->fDebug = lvl;
1369 }
1370
1371 //____________________________________________________________________
1372 void
1373 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1374 {
1375   //
1376   // Create sum histogram
1377   //
1378   // Parameters:
1379   //    data  Data histogram to clone
1380   //    mc    (optional) MC histogram to clone
1381   //
1382   DGUARD(fDebug,1,"Create centrality bin sums from %s", data->GetName());
1383   if (data) {
1384     fSum = new Sum(GetName(),"");
1385     fSum->Init(fSums, data, GetColor());
1386     fSum->fDebug = fDebug;
1387   }
1388
1389   // If no MC data is given, then do not create MC sum histogram
1390   if (!mc) return;
1391
1392   fSumMC = new Sum(GetName(), "MC");
1393   fSumMC->Init(fSums, mc, GetColor());
1394   fSumMC->fDebug = fDebug;
1395 }
1396
1397 //____________________________________________________________________
1398 Bool_t
1399 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1400                                              Int_t triggerMask,
1401                                              Double_t vzMin, Double_t vzMax)
1402 {
1403   //
1404   // Check the trigger, vertex, and centrality
1405   //
1406   // Parameters:
1407   //    aod Event input
1408   //
1409   // Return:
1410   //    true if the event is to be used
1411   //
1412   if (!forward) return false;
1413
1414   DGUARD(fDebug,2,"Check the event");
1415   // We do not check for centrality here - it's already done
1416   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
1417 }
1418
1419
1420 //____________________________________________________________________
1421 void
1422 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1423                                                Int_t triggerMask, Bool_t isZero,
1424                                                Double_t vzMin, Double_t vzMax,
1425                                                const TH2D* data, const TH2D* mc)
1426 {
1427   //
1428   // Process an event
1429   //
1430   // Parameters:
1431   //    forward     Forward data (for trigger, vertex, & centrality)
1432   //    triggerMask Trigger mask
1433   //    vzMin       Minimum IP z coordinate
1434   //    vzMax       Maximum IP z coordinate
1435   //    data        Data histogram
1436   //    mc          MC histogram
1437   //
1438   DGUARD(fDebug,1,"Process one event for %s a given centrality bin",
1439          data->GetName());
1440   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1441   if (!data) return;
1442   if (!fSum) CreateSums(data, mc);
1443
1444   fSum->Add(data, isZero);
1445   if (mc) fSumMC->Add(mc, isZero);
1446 }
1447
1448 //________________________________________________________________________
1449 Double_t
1450 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1451                                                 UShort_t    scheme,
1452                                                 Double_t    trigEff,
1453                                                 Double_t&   ntotal) const
1454 {
1455   //
1456   // Calculate normalization
1457   //
1458   // Parameters:
1459   //    t       Trigger histogram
1460   //    scheme  Normaliztion scheme
1461   //    trigEff From MC
1462   //    ntotal  On return, contains the number of events.
1463   //
1464   DGUARD(fDebug,1,"Normalize centrality bin with %s", t.GetName());
1465   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1466   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1467   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1468   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1469   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1470   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1471   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1472   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1473   Double_t nAccepted   = ntotal; // t.GetBinContent(AliAODForwardMult::kAccepted);
1474   ntotal               = 0;
1475
1476   if (nTriggered <= 0.1) {
1477     AliError("Number of triggered events <= 0");
1478     return -1;
1479   }
1480   if (nWithVertex <= 0.1) {
1481     AliError("Number of events with vertex <= 0");
1482     return -1;
1483   }
1484   ntotal          = nAccepted;
1485   Double_t vtxEff = nWithVertex / nTriggered;
1486   Double_t scaler = 1;
1487   Double_t beta   = nA + nC - 2*nE;
1488
1489   if (scheme & kEventLevel && !(scheme & kZeroBin)) {
1490     ntotal = nAccepted / vtxEff;
1491     scaler = vtxEff;
1492     AliInfo(Form("Calculating event normalisation as\n"
1493                  " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1494                  Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1495                  ntotal, scaler));
1496
1497     if (scheme & kBackground) {
1498       //          1            E_V             E_V
1499       //   s = --------- = ------------- = ------------
1500       //        1 - beta   1 - beta E_V    1 - beta N_V
1501       //       ---  ----       --------        ---- ---
1502       //       E_V  N_V          N_V           N_V  N_T
1503       //
1504       //          E_V
1505       //     = ------------
1506       //        1 - beta
1507       //            ----
1508       //             N_T
1509       //
1510       ntotal -= nAccepted * beta / nWithVertex;
1511       // This one is direct and correct.
1512       // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1513       // A simpler expresion
1514       scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1515       AliInfo(Form("Calculating event normalisation as\n"
1516                    " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1517                    " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1518                    Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1519                    nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1520                    Int_t(nWithVertex), ntotal, scaler));
1521     }
1522   }
1523   if (scheme & kZeroBin) {
1524     // Calculate as
1525     //
1526     //  N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1527     //    = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1528     //
1529     //  s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1530     //    = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1531     //
1532     if (!(scheme & kBackground)) beta = 0;
1533     ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1534                                          - beta / nWithVertex));
1535     scaler = nWithVertex / (nWithVertex +
1536                             1/trigEff * (nTriggered-nWithVertex-beta));
1537     AliInfo(Form("Calculating event normalisation as\n"
1538                  "  beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1539                  "  N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1540                  "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1541                  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1542                  Int_t(nAccepted), trigEff, Int_t(nTriggered),
1543                  Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1544                  ntotal, scaler));
1545   }
1546   if (scheme & kTriggerEfficiency && !(scheme & kZeroBin)) {
1547     ntotal /= trigEff;
1548     scaler *= trigEff;
1549     AliInfo(Form("Correcting for trigger efficiency:\n"
1550                  " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1551                  trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1552   }
1553
1554   AliInfo(Form("\n"
1555                " Total of        %9d events for %s\n"
1556                "  of these       %9d have an offline trigger\n"
1557                "  of these N_T = %9d has the selected trigger\n"
1558                "  of these N_V = %9d has a vertex\n"
1559                "  of these N_A = %9d were in the selected range\n"
1560                "  Triggers by hardware type:\n"
1561                "    N_b   =  %9d\n"
1562                "    N_ac  =  %9d (%d+%d)\n"
1563                "    N_e   =  %9d\n"
1564                "  Vertex efficiency:          %f\n"
1565                "  Trigger efficiency:         %f\n"
1566                "  Total number of events: N = %f\n"
1567                "  Scaler (N_A/N):             %f",
1568                Int_t(nAll), GetTitle(), Int_t(nOffline),
1569                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1570                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1571                vtxEff, trigEff, ntotal, scaler));
1572   return scaler;
1573 }
1574
1575 //________________________________________________________________________
1576 const char*
1577 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1578                                                 Bool_t sym,
1579                                                 const char* postfix) const
1580 {
1581   static TString n;
1582   n = Form("dndeta%s%s",GetName(), postfix);
1583   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1584   if (sym)       n.Append("_mirror");
1585   return n.Data();
1586 }
1587 //________________________________________________________________________
1588 TH1*
1589 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1590                                             Bool_t sym,
1591                                             const char* postfix) const
1592 {
1593   if (!fOutput) {
1594     AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
1595                     fLow, fHigh));
1596     return 0;
1597   }
1598   TString  n = GetResultName(rebin, sym, postfix);
1599   TObject* o = fOutput->FindObject(n.Data());
1600   if (!o) {
1601     // AliWarning(Form("Object %s not found in output list", n.Data()));
1602     return 0;
1603   }
1604   return static_cast<TH1*>(o);
1605 }
1606
1607 //________________________________________________________________________
1608 void
1609 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1610                                              const char* postfix,
1611                                              bool        rootProj,
1612                                              bool        corrEmpty,
1613                                              const TH1*  shapeCorr,
1614                                              Double_t    scaler,
1615                                              bool        symmetrice,
1616                                              Int_t       rebin,
1617                                              bool        cutEdges,
1618                                              Int_t       marker,
1619                                              Int_t       color,
1620                                              TList*      mclist,
1621                                              TList*      truthlist)
1622 {
1623   //
1624   // Generate the dN/deta result from input
1625   //
1626   // Parameters:
1627   //     sum        Sum of 2D hists
1628   //     postfix    Post fix on names
1629   //     rootProj   Whether to use ROOT TH2::ProjectionX
1630   //     corrEmpty  Correct for empty bins
1631   //     shapeCorr  Shape correction to use
1632   //     scaler     Event-level normalization scaler
1633   //     symmetrice Whether to make symmetric extensions
1634   //     rebin      Whether to rebin
1635   //     cutEdges   Whether to cut edges when rebinning
1636   //
1637   DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1638   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1639                                                      GetName(), postfix)));
1640   TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0,
1641                            rootProj, corrEmpty, false);
1642   accNorm->SetDirectory(0);
1643
1644   // ---- Scale by shape correction ----------------------------------
1645   if (shapeCorr) copy->Divide(shapeCorr);
1646   else AliInfo("No shape correction specified, or disabled");
1647
1648   // --- Normalize to the coverage -----------------------------------
1649   ScaleToCoverage(copy, accNorm);
1650
1651   // --- Event-level normalization -----------------------------------
1652   copy->Scale(scaler);
1653
1654   // --- Project on X axis -------------------------------------------
1655   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1656                           1, sum->GetNbinsY(), rootProj, corrEmpty);
1657   dndeta->SetDirectory(0);
1658   // Event-level normalization
1659   dndeta->Scale(1., "width");
1660   copy->Scale(1., "width");
1661
1662   TH1D* dndetaMCCorrection = 0;
1663   TList* centlist          = 0;
1664   TH1D* dndetaMCtruth      = 0;
1665   TList* truthcentlist     = 0;
1666
1667   // Possible final correction to <MC analysis> / <MC truth>
1668   if(mclist)
1669     centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1670   if(centlist)
1671     dndetaMCCorrection = static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",GetName(), postfix)));
1672   if(truthlist)
1673     truthcentlist = static_cast<TList*> (truthlist->FindObject(GetListName()));
1674   if(truthcentlist)
1675     dndetaMCtruth =  static_cast<TH1D*> (truthcentlist->FindObject("dndetaTruth"));
1676   //std::cout<<dndetaMCCorrection<<"  "<<dndetaMCtruth<<std::endl;
1677   if(dndetaMCCorrection && dndetaMCtruth) {
1678     AliInfo("Correcting with final MC correction");
1679     dndetaMCCorrection->Divide(dndetaMCtruth);
1680     dndeta->Divide(dndetaMCCorrection);
1681
1682     //std::cout<<"histo "<<Form("dndeta%s%s",GetName(), postfix)<<"  "<<GetListName()<<"  "<<dndetaMCCorrection<<std::endl;
1683     //std::cout<<"truth "<<GetListName()<<"  "<<dndetaMCtruth<<std::endl;
1684
1685   }
1686   else AliInfo("No final MC correction applied");
1687
1688   // --- Set some histogram attributes -------------------------------
1689   TString post;
1690   Int_t rColor = GetColor(color);
1691   if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
1692   SetHistogramAttributes(dndeta,  rColor, marker,
1693                          Form("ALICE %s%s", GetName(), post.Data()));
1694   SetHistogramAttributes(accNorm, rColor, marker,
1695                          Form("ALICE %s normalisation%s",
1696                               GetName(), post.Data()));
1697
1698   // --- Make symmetric extensions and rebinnings --------------------
1699   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
1700   fOutput->Add(dndeta);
1701   fOutput->Add(accNorm);
1702   fOutput->Add(copy);
1703   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1704   if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1705 }
1706
1707 //________________________________________________________________________
1708 void
1709 AliBasedNdetaTask::CentralityBin::End(TList*      sums,
1710                                       TList*      results,
1711                                       UShort_t    scheme,
1712                                       const TH1*  shapeCorr,
1713                                       Double_t    trigEff,
1714                                       Bool_t      symmetrice,
1715                                       Int_t       rebin,
1716                                       Bool_t      rootProj,
1717                                       Bool_t      corrEmpty,
1718                                       Bool_t      cutEdges,
1719                                       Int_t       triggerMask,
1720                                       Int_t       marker,
1721                                       Int_t       color,
1722                                       TList*      mclist,
1723                                       TList*      truthlist)
1724 {
1725   //
1726   // End of processing
1727   //
1728   // Parameters:
1729   //    sums        List of sums
1730   //    results     Output list of results
1731   //    shapeCorr   Shape correction or nil
1732   //    trigEff     Trigger efficiency
1733   //    symmetrice  Whether to symmetrice the results
1734   //    rebin       Whether to rebin the results
1735   //    corrEmpty   Whether to correct for empty bins
1736   //    cutEdges    Whether to cut edges when rebinning
1737   //    triggerMask Trigger mask
1738   //
1739   DGUARD(fDebug,1,"End centrality bin procesing");
1740
1741   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1742   if(!fSums) {
1743     AliError("Could not retrieve TList fSums");
1744     return;
1745   }
1746
1747   fOutput = new TList;
1748   fOutput->SetName(GetListName());
1749   fOutput->SetOwner();
1750   results->Add(fOutput);
1751
1752   if (!fSum) {
1753     AliInfo("This task did not produce any output");
1754     return;
1755   }
1756
1757   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1758
1759   if (!fTriggers) {
1760     AliError("Couldn't find histogram 'triggers' in list");
1761     return;
1762   }
1763   if (!fSum) {
1764     AliError(Form("No sum object for %s", GetName()));
1765     return;
1766   }
1767
1768   // --- Get normalization scaler ------------------------------------
1769   Double_t epsilonT  = trigEff;
1770   Double_t epsilonT0 = trigEff;
1771   // TEMPORARY FIX
1772   if (triggerMask == AliAODForwardMult::kNSD) {
1773     // This is a local change
1774     epsilonT = 0.96;
1775     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1776   }
1777   else if (triggerMask == AliAODForwardMult::kInel) {
1778     // This is a local change
1779     epsilonT = 0.934;
1780     AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
1781   }
1782   if (scheme & kZeroBin) {
1783     if (triggerMask==AliAODForwardMult::kInel)
1784       epsilonT0 = 0.785021; // 0.100240;
1785     else if (triggerMask==AliAODForwardMult::kInelGt0)
1786       epsilonT0 = 0;
1787     else if (triggerMask==AliAODForwardMult::kNSD)
1788       epsilonT0 = .706587;
1789     epsilonT = 1;
1790     AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
1791                     epsilonT0));
1792   }
1793
1794   // Get our histograms
1795   Double_t nSum   = 0;
1796   TH2D*    sum    = fSum->GetSum(fSums, fOutput, nSum, epsilonT0, 1,
1797                                  marker, rootProj, corrEmpty);
1798   Double_t nSumMC = 0;
1799   TH2D*    sumMC  = 0;
1800   if (fSumMC) sumMC = fSumMC->GetSum(fSums, fOutput, nSumMC,
1801                                      epsilonT0, 1, marker,
1802                                      rootProj, corrEmpty);
1803   if (!sum) {
1804     AliError("Failed to get sum from summer - bailing out");
1805     return;
1806   }
1807
1808   Double_t ntotal = nSum;
1809   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1810   if (scaler < 0) {
1811     AliError("Failed to calculate normalization - bailing out");
1812     return;
1813   }
1814   fOutput->Add(fTriggers->Clone());
1815
1816   // --- Make result and store ---------------------------------------
1817   MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1818              scaler, symmetrice, rebin, cutEdges, marker, color,
1819              mclist, truthlist);
1820
1821   // --- Process result from TrackRefs -------------------------------
1822   if (sumMC)
1823     MakeResult(sumMC, "MC", rootProj, corrEmpty,
1824                (scheme & kShape) ? shapeCorr : 0,
1825                scaler, symmetrice, rebin, cutEdges,
1826                GetMarkerStyle(GetMarkerBits(marker)+4), color,
1827                mclist, truthlist);
1828
1829   // Temporary stuff
1830   // if (!IsAllBin()) return;
1831
1832 }
1833
1834 //
1835 // EOF
1836 //