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