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