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