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