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