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