PWG2/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 <TParameter.h>
   9 #include <AliAnalysisManager.h>
  10 #include <AliAODEvent.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODInputHandler.h>
  13 #include "AliForwardUtil.h"
  14 #include "AliAODForwardMult.h"
  15 #include <TFile.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     fShapeCorr(0),
  32     fListOfCentralities(0),
  33     fSNNString(0),
  34     fSysString(0),
  35     fCent(0),
  36     fCentAxis(0),
  37     fNormalizationScheme(kFull),
  38     fSchemeString(0),
  39     fTriggerString(0)
  40 {
  41   //
  42   // Constructor
  43   //
  44 }
  45
  46 //____________________________________________________________________
  47 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  48   : AliAnalysisTaskSE(name),
  49     fSums(0),           // Container of sums
  50     fOutput(0),         // Container of output
  51     fVtxMin(-10),       // Minimum v_z
  52     fVtxMax(10),        // Maximum v_z
  53     fTriggerMask(AliAODForwardMult::kInel),
  54     fRebin(5),          // Rebinning factor
  55     fCutEdges(false),
  56     fSymmetrice(true),
  57     fCorrEmpty(true),
  58     fUseROOTProj(false),
  59     fTriggerEff(1),
  60     fShapeCorr(0),
  61     fListOfCentralities(0),
  62     fSNNString(0),
  63     fSysString(0),
  64     fCent(0),
  65     fCentAxis(0),
  66     fNormalizationScheme(kFull),
  67     fSchemeString(0),
  68     fTriggerString(0)
  69 {
  70   //
  71   // Constructor
  72   //
  73   fListOfCentralities = new TObjArray(1);
  74
  75   // Set the normalisation scheme
  76   SetNormalizationScheme(kFull);
  77
  78   // Set the trigger mask
  79   SetTriggerMask(AliAODForwardMult::kInel);
  80
  81   // Output slot #1 writes into a TH1 container
  82   DefineOutput(1, TList::Class());
  83   DefineOutput(2, TList::Class());
  84 }
  85
  86 //____________________________________________________________________
  87 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
  88   : AliAnalysisTaskSE(o),
  89     fSums(o.fSums),             // TList* - Container of sums
  90     fOutput(o.fOutput),         // Container of output
  91     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
  92     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
  93     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
  94     fRebin(o.fRebin),           // Int_t - Rebinning factor
  95     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
  96     fSymmetrice(o.fSymmetrice),
  97     fCorrEmpty(o.fCorrEmpty),
  98     fUseROOTProj(o.fUseROOTProj),
  99     fTriggerEff(o.fTriggerEff),
 100     fShapeCorr(o.fShapeCorr),
 101     fListOfCentralities(o.fListOfCentralities),
 102     fSNNString(o.fSNNString),
 103     fSysString(o.fSysString),
 104     fCent(o.fCent),
 105     fCentAxis(o.fCentAxis),
 106     fNormalizationScheme(o.fNormalizationScheme),
 107     fSchemeString(o.fSchemeString),
 108     fTriggerString(o.fTriggerString)
 109 {}
 110
 111 //____________________________________________________________________
 112 AliBasedNdetaTask::~AliBasedNdetaTask()
 113 {
 114   //
 115   // Destructor
 116   //
 117   if (fSums) {
 118     fSums->Delete();
 119     delete fSums;
 120     fSums = 0;
 121   }
 122   if (fOutput) {
 123     fOutput->Delete();
 124     delete fOutput;
 125     fOutput = 0;
 126   }
 127 }
 128
 129 //________________________________________________________________________
 130 void
 131 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
 132 {
 133   if (!fCentAxis) {
 134     fCentAxis = new TAxis();
 135     fCentAxis->SetName("centAxis");
 136     fCentAxis->SetTitle("Centrality [%]");
 137   }
 138   TArrayD dbins(n+1);
 139   for (UShort_t i = 0; i <= n; i++)
 140     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
 141   fCentAxis->Set(n, dbins.GetArray());
 142 }
 143
 144 //________________________________________________________________________
 145 void
 146 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
 147 {
 148   //
 149   // Add a centrality bin
 150   //
 151   // Parameters:
 152   //    low  Low cut
 153   //    high High cut
 154   //
 155   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
 156   AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
 157   fListOfCentralities->AddAtAndExpand(bin, at);
 158 }
 159
 160 //________________________________________________________________________
 161 AliBasedNdetaTask::CentralityBin*
 162 AliBasedNdetaTask::MakeCentralityBin(const char* name,
 163                                      Short_t low, Short_t high) const
 164 {
 165   //
 166   // Make a centrality bin
 167   //
 168   // Parameters:
 169   //    name  Name used for histograms
 170   //    low   Low cut in percent
 171   //    high  High cut in percent
 172   //
 173   // Return:
 174   //    A newly created centrality bin
 175   //
 176   return new CentralityBin(name, low, high);
 177 }
 178 //________________________________________________________________________
 179 void
 180 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
 181 {
 182   //
 183   // Set normalisation scheme
 184   //
 185   UShort_t    scheme = 0;
 186   TString     twhat(what);
 187   twhat.ToUpper();
 188   TObjString* opt;
 189   TIter       next(twhat.Tokenize(" ,|"));
 190   while ((opt = static_cast<TObjString*>(next()))) {
 191     TString s(opt->GetString());
 192     if      (s.IsNull()) continue;
 193     Bool_t add = true;
 194     switch (s[0]) {
 195     case '-': add = false; // Fall through
 196     case '+': s.Remove(0,1);  // Remove character
 197     }
 198     UShort_t bit = 0;
 199     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
 200     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
 201     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
 202     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
 203     else if (s.CompareTo("FULL")      == 0) bit = kFull;
 204     else if (s.CompareTo("NONE")      == 0) bit = kNone;
 205     else
 206       Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
 207     if (add) scheme |= bit;
 208     else     scheme ^= bit;
 209   }
 210   SetNormalizationScheme(scheme);
 211 }
 212 //________________________________________________________________________
 213 void
 214 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
 215 {
 216   fNormalizationScheme = scheme;
 217   TString tit = "";
 218   if (scheme == kFull) tit = "FULL";
 219   else {
 220     if (scheme & kEventLevel)        tit.Append("EVENT ");
 221     if (scheme & kShape)             tit.Append("SHAPE ");
 222     if (scheme & kBackground)        tit.Append("BACKGROUND ");
 223     if (scheme & kTriggerEfficiency) tit.Append("TRIGGER");
 224   }
 225   tit = tit.Strip(TString::kBoth);
 226   if (!fSchemeString) fSchemeString = new TNamed("scheme", "");
 227   fSchemeString->SetTitle(tit);
 228   fSchemeString->SetUniqueID(fNormalizationScheme);
 229 }
 230 //________________________________________________________________________
 231 void
 232 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 233 {
 234   //
 235   // Set the trigger maskl
 236   //
 237   // Parameters:
 238   //    mask Trigger mask
 239   //
 240   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
 241 }
 242 //________________________________________________________________________
 243 void
 244 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
 245 {
 246   fTriggerMask = mask;
 247   TString tit(AliAODForwardMult::GetTriggerString(mask));
 248   tit = tit.Strip(TString::kBoth);
 249   if (!fTriggerString) fTriggerString = new TNamed("trigger", "");
 250   fTriggerString->SetTitle(tit);
 251   fTriggerString->SetUniqueID(fTriggerMask);
 252 }
 253
 254 //________________________________________________________________________
 255 void
 256 AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
 257 {
 258   //
 259   // Set the shape correction (a.k.a., track correction) for selected
 260   // trigger(s)
 261   //
 262   // Parameters:
 263   //    h Correction
 264   //
 265   if (!c) return;
 266   fShapeCorr = static_cast<TH1*>(c->Clone());
 267   fShapeCorr->SetDirectory(0);
 268 }
 269
 270 //________________________________________________________________________
 271 void
 272 AliBasedNdetaTask::UserCreateOutputObjects()
 273 {
 274   //
 275   // Create output objects.
 276   //
 277   // This is called once per slave process
 278   //
 279   fSums = new TList;
 280   fSums->SetName(Form("%s_sums", GetName()));
 281   fSums->SetOwner();
 282
 283   // Automatically add 'all' centrality bin if nothing has been defined.
 284   AddCentralityBin(0, 0, 0);
 285   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
 286     const TArrayD* bins = fCentAxis->GetXbins();
 287     Int_t          nbin = fCentAxis->GetNbins();
 288     for (Int_t i = 0; i < nbin; i++)
 289       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
 290   }
 291   fSums->Add(fCentAxis);
 292
 293
 294   // Centrality histogram
 295   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
 296   fCent->SetDirectory(0);
 297   fCent->SetXTitle(0);
 298   fSums->Add(fCent);
 299
 300   // Loop over centrality bins
 301   TIter next(fListOfCentralities);
 302   CentralityBin* bin = 0;
 303   while ((bin = static_cast<CentralityBin*>(next())))
 304     bin->CreateOutputObjects(fSums);
 305
 306   // Check that we have an AOD input handler
 307   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 308   AliAODInputHandler* ah =
 309     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
 310   if (!ah) AliFatal("No AOD input handler set in analysis manager");
 311
 312   // Post data for ALL output slots >0 here, to get at least an empty histogram
 313   PostData(1, fSums);
 314 }
 315 //____________________________________________________________________
 316 void
 317 AliBasedNdetaTask::UserExec(Option_t *)
 318 {
 319   //
 320   // Process a single event
 321   //
 322   // Parameters:
 323   //    option Not used
 324   //
 325   // Main loop
 326   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 327   if (!aod) {
 328     AliError("Cannot get the AOD event");
 329     return;
 330   }
 331
 332   TObject* obj = aod->FindListObject("Forward");
 333   if (!obj) {
 334     AliWarning("No forward object found");
 335     return;
 336   }
 337   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 338
 339   // Fill centrality histogram
 340   Float_t cent = forward->GetCentrality();
 341   fCent->Fill(cent);
 342
 343   // Get the histogram(s)
 344   TH2D* data   = GetHistogram(aod, false);
 345   TH2D* dataMC = GetHistogram(aod, true);
 346
 347   // Loop over centrality bins
 348   CentralityBin* allBin =
 349     static_cast<CentralityBin*>(fListOfCentralities->At(0));
 350   allBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
 351
 352   // Find this centrality bin
 353   if (fCentAxis && fCentAxis->GetNbins() > 0) {
 354     Int_t          icent   = fCentAxis->FindBin(cent);
 355     CentralityBin* thisBin = 0;
 356     if (icent >= 1 && icent <= fCentAxis->GetNbins())
 357       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
 358     if (thisBin)
 359       thisBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax,
 360                             data, dataMC);
 361   }
 362
 363   // Here, we get the update
 364   if (!fSNNString) {
 365     UShort_t sNN = forward->GetSNN();
 366     fSNNString = new TNamed("sNN", "");
 367     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
 368     fSNNString->SetUniqueID(sNN);
 369     fSums->Add(fSNNString);
 370
 371     UShort_t sys = forward->GetSystem();
 372     fSysString = new TNamed("sys", "");
 373     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
 374     fSysString->SetUniqueID(sys);
 375     fSums->Add(fSysString);
 376
 377     fSums->Add(fSchemeString);
 378     fSums->Add(fTriggerString);
 379
 380     // Print();
 381   }
 382
 383   PostData(1, fSums);
 384 }
 385
 386 //________________________________________________________________________
 387 void
 388 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 389                                           const char* title, const char* ytitle)
 390 {
 391   //
 392   // Set histogram graphical options, etc.
 393   //
 394   // Parameters:
 395   //    h       Histogram to modify
 396   //    colour  Marker color
 397   //    marker  Marker style
 398   //    title   Title of histogram
 399   //    ytitle  Title on y-axis.
 400   //
 401   h->SetTitle(title);
 402   h->SetMarkerColor(colour);
 403   h->SetMarkerStyle(marker);
 404   h->SetMarkerSize(1);
 405   h->SetFillStyle(0);
 406   h->SetYTitle(ytitle);
 407   h->GetXaxis()->SetTitleFont(132);
 408   h->GetXaxis()->SetLabelFont(132);
 409   h->GetXaxis()->SetNdivisions(10);
 410   h->GetYaxis()->SetTitleFont(132);
 411   h->GetYaxis()->SetLabelFont(132);
 412   h->GetYaxis()->SetNdivisions(10);
 413   h->GetYaxis()->SetDecimals();
 414   h->SetStats(0);
 415 }
 416
 417 //________________________________________________________________________
 418 TH1D*
 419 AliBasedNdetaTask::ProjectX(const TH2D* h,
 420                             const char* name,
 421                             Int_t firstbin,
 422                             Int_t lastbin,
 423                             bool  useRoot,
 424                             bool  corr,
 425                             bool  error)
 426 {
 427   //
 428   // Project onto the X axis
 429   //
 430   // Parameters:
 431   //    h         2D histogram
 432   //    name      New name
 433   //    firstbin  First bin to use
 434   //    lastbin   Last bin to use
 435   //    error     Whether to calculate errors
 436   //
 437   // Return:
 438   //    Newly created histogram or null
 439   //
 440   if (!h) return 0;
 441   if (useRoot)
 442     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 443
 444   TAxis* xaxis = h->GetXaxis();
 445   TAxis* yaxis = h->GetYaxis();
 446   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 447                           xaxis->GetXmin(), xaxis->GetXmax());
 448   static_cast<const TAttLine*>(h)->Copy(*ret);
 449   static_cast<const TAttFill*>(h)->Copy(*ret);
 450   static_cast<const TAttMarker*>(h)->Copy(*ret);
 451   ret->GetXaxis()->ImportAttributes(xaxis);
 452
 453   Int_t first = firstbin;
 454   Int_t last  = lastbin;
 455   if (first < 0)                         first = 0;
 456   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
 457   if      (last  < 0)                    last  = yaxis->GetNbins();
 458   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
 459   if (last-first < 0) {
 460     AliWarningGeneral("AliBasedNdetaTask",
 461                       Form("Nothing to project [%d,%d]", first, last));
 462     return 0;
 463
 464   }
 465
 466   // Loop over X bins
 467   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 468   Int_t ybins = (last-first+1);
 469   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 470     Double_t content = 0;
 471     Double_t error2  = 0;
 472     Int_t    nbins   = 0;
 473
 474
 475     for (Int_t ybin = first; ybin <= last; ybin++) {
 476       Double_t c1 = h->GetCellContent(xbin, ybin);
 477       Double_t e1 = h->GetCellError(xbin, ybin);
 478
 479       // Ignore empty bins
 480       if (c1 < 1e-12) continue;
 481       if (e1 < 1e-12) {
 482         if (error) continue;
 483         e1 = 1;
 484       }
 485
 486       content    += c1;
 487       error2     += e1*e1;
 488       nbins++;
 489     } // for (ybin)
 490     if(content > 0 && nbins > 0) {
 491       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 492 #if 0
 493       AliWarningGeneral(ret->GetName(),
 494                         Form("factor @ %d is %d/%d -> %f",
 495                              xbin, ybins, nbins, factor));
 496 #endif
 497       if (error) {
 498         // Calculate weighted average
 499         ret->SetBinContent(xbin, content * factor);
 500         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 501       }
 502       else
 503         ret->SetBinContent(xbin, factor * content);
 504     }
 505   } // for (xbin)
 506
 507   return ret;
 508 }
 509
 510 //________________________________________________________________________
 511 void
 512 AliBasedNdetaTask::Terminate(Option_t *)
 513 {
 514   //
 515   // Called at end of event processing..
 516   //
 517   // This is called once in the master
 518   //
 519   // Parameters:
 520   //    option Not used
 521   //
 522   // Draw result to screen, or perform fitting, normalizations Called
 523   // once at the end of the query
 524
 525   fSums = dynamic_cast<TList*> (GetOutputData(1));
 526   if(!fSums) {
 527     AliError("Could not retrieve TList fSums");
 528     return;
 529   }
 530
 531   fOutput = new TList;
 532   fOutput->SetName(Form("%s_result", GetName()));
 533   fOutput->SetOwner();
 534
 535   fSNNString     = static_cast<TNamed*>(fSums->FindObject("sNN"));
 536   fSysString     = static_cast<TNamed*>(fSums->FindObject("sys"));
 537   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
 538   fSchemeString  = static_cast<TNamed*>(fSums->FindObject("scheme"));
 539   fTriggerString = static_cast<TNamed*>(fSums->FindObject("trigger"));
 540
 541   if(fSysString && fSNNString &&
 542      fSysString->GetUniqueID() == AliForwardUtil::kPP)
 543     LoadNormalizationData(fSysString->GetUniqueID(),
 544                           fSNNString->GetUniqueID());
 545   // Print before we loop
 546   Print();
 547
 548   // Loop over centrality bins
 549   TIter next(fListOfCentralities);
 550   CentralityBin* bin = 0;
 551   while ((bin = static_cast<CentralityBin*>(next())))
 552     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
 553              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
 554              fTriggerMask, GetColor(), GetMarker());
 555
 556   // Output collision energy string
 557   if (fSNNString) fOutput->Add(fSNNString->Clone());
 558
 559   // Output collision system string
 560   if (fSysString) fOutput->Add(fSysString->Clone());
 561
 562   // Output centrality axis
 563   if (fCentAxis) fOutput->Add(fCentAxis);
 564
 565   // Output trigger string
 566   if (fTriggerString) fOutput->Add(fTriggerString->Clone());
 567
 568   // Normalization string
 569   if (fSchemeString) fOutput->Add(fSchemeString->Clone());
 570
 571   // Output vertex axis
 572   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 573   vtxAxis->SetName("vtxAxis");
 574   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 575   fOutput->Add(vtxAxis);
 576
 577   // Output trigger efficiency and shape correction
 578   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
 579   if (fShapeCorr) fOutput->Add(fShapeCorr);
 580
 581   TNamed* options = new TNamed("options","");
 582   TString str;
 583   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
 584   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
 585   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
 586   options->SetTitle(str);
 587   fOutput->Add(options);
 588
 589   PostData(2, fOutput);
 590 }
 591 //________________________________________________________________________
 592 void
 593 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 594 {
 595   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
 596   TString type("pp");
 597   TString snn("900");
 598   if(energy == 7000) snn.Form("7000");
 599   if(energy == 2750) snn.Form("2750");
 600
 601   if(fShapeCorr && (fTriggerEff != 1)) {
 602     AliInfo("Objects already set for normalization - no action taken");
 603     return;
 604   }
 605
 606   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/"
 607                                 "Normalization/normalizationHists_%s_%s.root",
 608                                 type.Data(),snn.Data()));
 609   if(!fin) {
 610     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
 611     return;
 612   }
 613
 614   // Shape correction
 615   TString shapeCorName(Form("h%sNormalization",
 616                             fTriggerMask == AliAODForwardMult::kInel ? "Inel" :
 617                             fTriggerMask == AliAODForwardMult::kNSD ? "NSD" :
 618                             fTriggerMask == AliAODForwardMult::kInelGt0 ?
 619                             "InelGt0" : "All"));
 620   TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
 621   if (shapeCor) SetShapeCorrection(shapeCor);
 622
 623   // Trigger efficiency
 624   TString effName(Form("%sTriggerEff",
 625                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 626                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 627                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
 628                        "inelgt0" : "all"));
 629   TParameter<float>* eff = static_cast<TParameter<float>*>(fin->Get(effName));
 630   if (eff) SetTriggerEff(eff->GetVal());
 631
 632   // TEMPORARY FIX
 633   // Rescale the shape correction by the trigger efficiency
 634   AliWarning(Form("Rescaling shape correction by trigger efficiency: "
 635                   "1/E_X=1/%f", fTriggerEff));
 636   fShapeCorr->Scale(1. / fTriggerEff);
 637
 638   // Print - out
 639   if (shapeCor && eff) AliInfo("Loaded objects for normalization.");
 640 }
 641
 642
 643 //________________________________________________________________________
 644 void
 645 AliBasedNdetaTask::Print(Option_t*) const
 646 {
 647   //
 648   // Print information
 649   //
 650   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
 651             << std::boolalpha
 652             << " Trigger:              " << (fTriggerString ?
 653                                              fTriggerString->GetTitle() :
 654                                              "none") << "\n"
 655             << " Vertex range:         [" << fVtxMin << ":" << fVtxMax << "]\n"
 656             << " Rebin factor:         " << fRebin << "\n"
 657             << " Cut edges:            " << fCutEdges << "\n"
 658             << " Symmertrice:          " << fSymmetrice << "\n"
 659             << " Use TH2::ProjectionX: " << fUseROOTProj << "\n"
 660             << " Correct for empty:    " << fCorrEmpty << "\n"
 661             << " Normalization scheme: " << (fSchemeString ?
 662                                              fSchemeString->GetTitle() :
 663                                              "none") <<"\n"
 664             << " Trigger efficiency:   " << fTriggerEff << "\n"
 665             << " Shape correction:     " << (fShapeCorr ?
 666                                            fShapeCorr->GetName() :
 667                                            "none") << "\n"
 668             << " sqrt(s_NN):           " << (fSNNString ?
 669                                              fSNNString->GetTitle() :
 670                                            "unknown") << "\n"
 671             << " Collision system:     " << (fSysString ?
 672                                              fSysString->GetTitle() :
 673                                            "unknown") << "\n"
 674             << " Centrality bins:      " << (fCentAxis ? "" : "none");
 675   if (fCentAxis) {
 676     Int_t           nBins = fCentAxis->GetNbins();
 677     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
 678     for (Int_t i = 0; i <= nBins; i++)
 679       std::cout << (i==0 ? "" : "-") << bins[i];
 680   }
 681   std::cout << std::noboolalpha << std::endl;
 682
 683 }
 684
 685 //________________________________________________________________________
 686 TH1D*
 687 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
 688 {
 689   //
 690   // Make a copy of the input histogram and rebin that histogram
 691   //
 692   // Parameters:
 693   //    h  Histogram to rebin
 694   //
 695   // Return:
 696   //    New (rebinned) histogram
 697   //
 698   if (rebin <= 1) return 0;
 699
 700   Int_t nBins = h->GetNbinsX();
 701   if(nBins % rebin != 0) {
 702     AliWarningGeneral("AliBasedNdetaTask",
 703                       Form("Rebin factor %d is not a devisor of current number "
 704                            "of bins %d in the histogram %s",
 705                            rebin, nBins, h->GetName()));
 706     return 0;
 707   }
 708
 709   // Make a copy
 710   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
 711                                                h->GetName(), rebin)));
 712   tmp->Rebin(rebin);
 713   tmp->SetDirectory(0);
 714
 715   // The new number of bins
 716   Int_t nBinsNew = nBins / rebin;
 717   for(Int_t i = 1;i<= nBinsNew; i++) {
 718     Double_t content = 0;
 719     Double_t sumw    = 0;
 720     Double_t wsum    = 0;
 721     Int_t    nbins   = 0;
 722     for(Int_t j = 1; j<=rebin;j++) {
 723       Int_t    bin = (i-1)*rebin + j;
 724       Double_t c   =  h->GetBinContent(bin);
 725       if (c <= 0) continue;
 726
 727       if (cutEdges) {
 728         if (h->GetBinContent(bin+1)<=0 ||
 729             h->GetBinContent(bin-1)<=0) {
 730 #if 0
 731           AliWarningGeneral("AliBasedNdetaTask",
 732                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
 733                                  bin, c, h->GetName(),
 734                                  bin+1, h->GetBinContent(bin+1),
 735                                  bin-1, h->GetBinContent(bin-1)));
 736 #endif
 737           continue;
 738         }
 739       }
 740       Double_t e =  h->GetBinError(bin);
 741       Double_t w =  1 / (e*e); // 1/c/c
 742       content    += c;
 743       sumw       += w;
 744       wsum       += w * c;
 745       nbins++;
 746     }
 747
 748     if(content > 0 && nbins > 0) {
 749       tmp->SetBinContent(i, wsum / sumw);
 750       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 751     }
 752   }
 753
 754   return tmp;
 755 }
 756
 757 //__________________________________________________________________
 758 TH1*
 759 AliBasedNdetaTask::Symmetrice(const TH1* h)
 760 {
 761   //
 762   // Make an extension of @a h to make it symmetric about 0
 763   //
 764   // Parameters:
 765   //    h Histogram to symmertrice
 766   //
 767   // Return:
 768   //    Symmetric extension of @a h
 769   //
 770   Int_t nBins = h->GetNbinsX();
 771   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 772   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 773   s->Reset();
 774   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 775   s->SetMarkerColor(h->GetMarkerColor());
 776   s->SetMarkerSize(h->GetMarkerSize());
 777   s->SetMarkerStyle(h->GetMarkerStyle()+4);
 778   s->SetFillColor(h->GetFillColor());
 779   s->SetFillStyle(h->GetFillStyle());
 780   s->SetDirectory(0);
 781
 782   // Find the first and last bin with data
 783   Int_t first = nBins+1;
 784   Int_t last  = 0;
 785   for (Int_t i = 1; i <= nBins; i++) {
 786     if (h->GetBinContent(i) <= 0) continue;
 787     first = TMath::Min(first, i);
 788     last  = TMath::Max(last,  i);
 789   }
 790
 791   Double_t xfirst = h->GetBinCenter(first-1);
 792   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 793   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 794   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 795     s->SetBinContent(j, h->GetBinContent(i));
 796     s->SetBinError(j, h->GetBinError(i));
 797   }
 798   // Fill in overlap bin
 799   s->SetBinContent(l2+1, h->GetBinContent(first));
 800   s->SetBinError(l2+1, h->GetBinError(first));
 801   return s;
 802 }
 803
 804 //====================================================================
 805 AliBasedNdetaTask::CentralityBin::CentralityBin()
 806   : TNamed("", ""),
 807     fSums(0),
 808     fOutput(0),
 809     fSum(0),
 810     fSumMC(0),
 811     fTriggers(0),
 812     fLow(0),
 813     fHigh(0)
 814 {
 815   //
 816   // Constructor
 817   //
 818 }
 819 //____________________________________________________________________
 820 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
 821                                                 Short_t low, Short_t high)
 822   : TNamed(name, ""),
 823     fSums(0),
 824     fOutput(0),
 825     fSum(0),
 826     fSumMC(0),
 827     fTriggers(0),
 828     fLow(low),
 829     fHigh(high)
 830 {
 831   //
 832   // Constructor
 833   //
 834   // Parameters:
 835   //    name Name used for histograms (e.g., Forward)
 836   //    low  Lower centrality cut in percent
 837   //    high Upper centrality cut in percent
 838   //
 839   if (low <= 0 && high <= 0) {
 840     fLow  = 0;
 841     fHigh = 0;
 842     SetTitle("All centralities");
 843   }
 844   else {
 845     fLow  = low;
 846     fHigh = high;
 847     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
 848   }
 849 }
 850 //____________________________________________________________________
 851 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
 852   : TNamed(o),
 853     fSums(o.fSums),
 854     fOutput(o.fOutput),
 855     fSum(o.fSum),
 856     fSumMC(o.fSumMC),
 857     fTriggers(o.fTriggers),
 858     fLow(o.fLow),
 859     fHigh(o.fHigh)
 860 {
 861   //
 862   // Copy constructor
 863   //
 864   // Parameters:
 865   //    other Object to copy from
 866   //
 867 }
 868 //____________________________________________________________________
 869 AliBasedNdetaTask::CentralityBin::~CentralityBin()
 870 {
 871   //
 872   // Destructor
 873   //
 874   if (fSums) fSums->Delete();
 875   if (fOutput) fOutput->Delete();
 876 }
 877
 878 //____________________________________________________________________
 879 AliBasedNdetaTask::CentralityBin&
 880 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
 881 {
 882   //
 883   // Assignment operator
 884   //
 885   // Parameters:
 886   //    other Object to assign from
 887   //
 888   // Return:
 889   //    Reference to this
 890   //
 891   SetName(o.GetName());
 892   SetTitle(o.GetTitle());
 893   fSums      = o.fSums;
 894   fOutput    = o.fOutput;
 895   fSum       = o.fSum;
 896   fSumMC     = o.fSumMC;
 897   fTriggers  = o.fTriggers;
 898   fLow       = o.fLow;
 899   fHigh      = o.fHigh;
 900
 901   return *this;
 902 }
 903 //____________________________________________________________________
 904 const char*
 905 AliBasedNdetaTask::CentralityBin::GetListName() const
 906 {
 907   //
 908   // Get the list name
 909   //
 910   // Return:
 911   //    List Name
 912   //
 913   if (IsAllBin()) return "all";
 914   return Form("cent%03d_%03d", fLow, fHigh);
 915 }
 916 //____________________________________________________________________
 917 void
 918 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
 919 {
 920   //
 921   // Create output objects
 922   //
 923   // Parameters:
 924   //    dir   Parent list
 925   //
 926   fSums = new TList;
 927   fSums->SetName(GetListName());
 928   fSums->SetOwner();
 929   dir->Add(fSums);
 930
 931   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
 932   fTriggers->SetDirectory(0);
 933   fSums->Add(fTriggers);
 934 }
 935 //____________________________________________________________________
 936 void
 937 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
 938 {
 939   //
 940   // Create sum histogram
 941   //
 942   // Parameters:
 943   //    data  Data histogram to clone
 944   //    mc    (optional) MC histogram to clone
 945   //
 946   fSum = static_cast<TH2D*>(data->Clone(GetName()));
 947   fSum->SetDirectory(0);
 948   fSum->Reset();
 949   fSums->Add(fSum);
 950
 951   // If no MC data is given, then do not create MC sum histogram
 952   if (!mc) return;
 953
 954   fSumMC = static_cast<TH2D*>(mc->Clone(Form("%sMC", GetName())));
 955   fSumMC->SetDirectory(0);
 956   fSumMC->Reset();
 957   fSums->Add(fSumMC);
 958 }
 959
 960 //____________________________________________________________________
 961 Bool_t
 962 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
 963                                              Int_t triggerMask,
 964                                              Double_t vzMin, Double_t vzMax)
 965 {
 966   //
 967   // Check the trigger, vertex, and centrality
 968   //
 969   // Parameters:
 970   //    aod Event input
 971   //
 972   // Return:
 973   //    true if the event is to be used
 974   //
 975   if (!forward) return false;
 976
 977   // We do not check for centrality here - it's already done
 978   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
 979 }
 980
 981
 982 //____________________________________________________________________
 983 void
 984 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
 985                                                Int_t triggerMask,
 986                                                Double_t vzMin, Double_t vzMax,
 987                                                const TH2D* data, const TH2D* mc)
 988 {
 989   //
 990   // Process an event
 991   //
 992   // Parameters:
 993   //    forward     Forward data (for trigger, vertex, & centrality)
 994   //    triggerMask Trigger mask
 995   //    vzMin       Minimum IP z coordinate
 996   //    vzMax       Maximum IP z coordinate
 997   //    data        Data histogram
 998   //    mc          MC histogram
 999   //
1000   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1001   if (!data) return;
1002   if (!fSum) CreateSums(data, mc);
1003
1004   fSum->Add(data);
1005   if (mc) fSumMC->Add(mc);
1006 }
1007
1008 //________________________________________________________________________
1009 Double_t
1010 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1011                                                 UShort_t    scheme,
1012                                                 Double_t    trigEff,
1013                                                 Double_t&   ntotal) const
1014 {
1015   //
1016   // Calculate normalization
1017   //
1018   // Parameters:
1019   //    t       Trigger histogram
1020   //    scheme  Normaliztion scheme
1021   //    trigEff From MC
1022   //    ntotal  On return, contains the number of events.
1023   //
1024   ntotal               = 0;
1025   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1026   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1027   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1028   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1029   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1030   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1031   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1032   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1033   Double_t nAccepted   = t.GetBinContent(AliAODForwardMult::kAccepted);
1034
1035   if (nTriggered <= 0.1) {
1036     AliError("Number of triggered events <= 0");
1037     return -1;
1038   }
1039   if (nWithVertex <= 0.1) {
1040     AliError("Number of events with vertex <= 0");
1041     return -1;
1042   }
1043   ntotal          = nAccepted;
1044   Double_t vtxEff = nWithVertex / nTriggered;
1045   Double_t scaler = 1;
1046   Double_t beta   = nA + nC - 2*nE;
1047
1048   if (scheme & kEventLevel) {
1049     ntotal = nAccepted / vtxEff;
1050     scaler = vtxEff;
1051     AliInfo(Form("Calculating event normalisation as\n"
1052                  " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1053                  Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1054                  ntotal, scaler));
1055
1056     if (scheme & kBackground) {
1057       //          1            E_V             E_V
1058       //   s = --------- = ------------- = ------------
1059       //        1 - beta   1 - beta E_V    1 - beta N_V
1060       //       ---  ----       --------        ---- ---
1061       //       E_V  N_V          N_V           N_V  N_T
1062       //
1063       //          E_V
1064       //     = ------------
1065       //        1 - beta
1066       //            ----
1067       //             N_T
1068       //
1069       ntotal -= nAccepted * beta / nWithVertex;
1070       // This one is direct and correct.
1071       // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1072       // A simpler expresion
1073       scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1074       AliInfo(Form("Calculating event normalisation as\n"
1075                    " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1076                    " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1077                    Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1078                    nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1079                    Int_t(nWithVertex), ntotal, scaler));
1080     }
1081   }
1082   if (scheme & kTriggerEfficiency) {
1083     ntotal /= trigEff;
1084     scaler *= trigEff;
1085     AliInfo(Form("Correcting for trigger efficiency:\n"
1086                  " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1087                  trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1088   }
1089
1090   AliInfo(Form("\n"
1091                " Total of        %9d events for %s\n"
1092                "  of these       %9d have an offline trigger\n"
1093                "  of these N_T = %9d has the selected trigger\n"
1094                "  of these N_V = %9d has a vertex\n"
1095                "  of these N_A = %9d were in the selected range\n"
1096                "  Triggers by hardware type:\n"
1097                "    N_b   =  %9d\n"
1098                "    N_ac  =  %9d (%d+%d)\n"
1099                "    N_e   =  %9d\n"
1100                "  Vertex efficiency:          %f\n"
1101                "  Trigger efficiency:         %f\n"
1102                "  Total number of events: N = %f\n"
1103                "  Scaler (N_A/N):             %f",
1104                Int_t(nAll), GetTitle(), Int_t(nOffline),
1105                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1106                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1107                vtxEff, trigEff, ntotal, scaler));
1108   return scaler;
1109 }
1110
1111 //________________________________________________________________________
1112 void
1113 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1114                                              const char* postfix,
1115                                              bool        rootProj,
1116                                              bool        corrEmpty,
1117                                              const TH1*  shapeCorr,
1118                                              Double_t    scaler,
1119                                              bool        symmetrice,
1120                                              Int_t       rebin,
1121                                              bool        cutEdges,
1122                                              Int_t       color,
1123                                              Int_t       marker)
1124 {
1125   //
1126   // Generate the dN/deta result from input
1127   //
1128   // Parameters:
1129   //     sum        Sum of 2D hists
1130   //     postfix    Post fix on names
1131   //     rootProj   Whether to use ROOT TH2::ProjectionX
1132   //     corrEmpty  Correct for empty bins
1133   //     shapeCorr  Shape correction to use
1134   //     scaler     Event-level normalization scaler
1135   //     symmetrice Whether to make symmetric extensions
1136   //     rebin      Whether to rebin
1137   //     cutEdges   Whether to cut edges when rebinning
1138   //
1139   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1140                                                      GetName(), postfix)));
1141   TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0,
1142                            rootProj, corrEmpty, false);
1143   accNorm->SetDirectory(0);
1144
1145   // ---- Scale by shape correction ----------------------------------
1146   if (shapeCorr) copy->Divide(shapeCorr);
1147   else AliInfo("No shape correction specified, or disabled");
1148
1149   // Normalize to the acceptance -
1150   // dndeta->Divide(accNorm);
1151   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
1152     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
1153       Double_t c = copy->GetBinContent(i, j);
1154       Double_t e = copy->GetBinError(i, j);
1155       Double_t a = accNorm->GetBinContent(i);
1156       copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
1157       copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
1158     }
1159   }
1160   // --- Event-level normalization -----------------------------------
1161   copy->Scale(scaler);
1162
1163   // --- Project on X axis -------------------------------------------
1164   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1165                           1, sum->GetNbinsY(), rootProj, corrEmpty);
1166   dndeta->SetDirectory(0);
1167   // Event-level normalization
1168   dndeta->Scale(1., "width");
1169   copy->Scale(1., "width");
1170
1171   // --- Set some histogram attributes -------------------------------
1172   SetHistogramAttributes(dndeta,  color, marker, Form("ALICE %s", GetName()));
1173   SetHistogramAttributes(accNorm, color, marker, Form("ALICE %s normalisation",
1174                                                       GetName()));
1175
1176   // --- Make symmetric extensions and rebinnings --------------------
1177   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
1178   fOutput->Add(dndeta);
1179   fOutput->Add(accNorm);
1180   fOutput->Add(copy);
1181   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1182   if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1183 }
1184
1185 //________________________________________________________________________
1186 void
1187 AliBasedNdetaTask::CentralityBin::End(TList*      sums,
1188                                       TList*      results,
1189                                       UShort_t    scheme,
1190                                       const TH1*  shapeCorr,
1191                                       Double_t    trigEff,
1192                                       Bool_t      symmetrice,
1193                                       Int_t       rebin,
1194                                       Bool_t      rootProj,
1195                                       Bool_t      corrEmpty,
1196                                       Bool_t      cutEdges,
1197                                       Int_t       triggerMask,
1198                                       Int_t       color,
1199                                       Int_t       marker)
1200 {
1201   //
1202   // End of processing
1203   //
1204   // Parameters:
1205   //    sums        List of sums
1206   //    results     Output list of results
1207   //    shapeCorr   Shape correction or nil
1208   //    trigEff     Trigger efficiency
1209   //    symmetrice  Whether to symmetrice the results
1210   //    rebin       Whether to rebin the results
1211   //    corrEmpty   Whether to correct for empty bins
1212   //    cutEdges    Whether to cut edges when rebinning
1213   //    triggerMask Trigger mask
1214   //
1215
1216   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1217   if(!fSums) {
1218     AliError("Could not retrieve TList fSums");
1219     return;
1220   }
1221
1222   fOutput = new TList;
1223   fOutput->SetName(GetListName());
1224   fOutput->SetOwner();
1225   results->Add(fOutput);
1226
1227   fSum      = static_cast<TH2D*>(fSums->FindObject(GetName()));
1228   fSumMC    = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
1229   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1230
1231   if (!fTriggers) {
1232     AliError("Couldn't find histogram 'triggers' in list");
1233     return;
1234   }
1235   if (!fSum) {
1236     AliError(Form("Couldn't find histogram '%s' in list", GetName()));
1237     return;
1238   }
1239
1240   // --- Get normalization scaler ------------------------------------
1241   Double_t ntotal = 0;
1242   Double_t epsilonT = trigEff;
1243   // TEMPORARY FIX
1244   if (triggerMask == AliAODForwardMult::kNSD) {
1245     // This is a local change
1246     epsilonT = 0.92;
1247     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1248   }
1249   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1250   if (scaler < 0) {
1251     AliError("Failed to calculate normalization - bailing out");
1252     return;
1253   }
1254   fOutput->Add(fTriggers->Clone());
1255
1256   // --- Make result and store ---------------------------------------
1257   MakeResult(fSum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1258              scaler, symmetrice, rebin, cutEdges, color, marker);
1259
1260   // --- Process result from TrackRefs -------------------------------
1261   if (fSumMC) {
1262     MakeResult(fSumMC, "MC", rootProj, corrEmpty,
1263                (scheme & kShape) ? shapeCorr : 0,
1264                scaler, symmetrice, rebin, cutEdges, color+2, marker);
1265   }
1266
1267   // Temporary stuff
1268   if (!IsAllBin()) return;
1269
1270   TFile* forward = TFile::Open("forward.root", "READ");
1271   if (!forward)  {
1272     AliWarning(Form("No forward.root file found"));
1273     return;
1274   }
1275
1276   TH1D* shapeCorrProj = 0;
1277   if (shapeCorr) {
1278     shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
1279     shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
1280     shapeCorrProj->SetDirectory(0);
1281     fOutput->Add(shapeCorrProj);
1282   }
1283
1284   TList* official = static_cast<TList*>(forward->Get("official"));
1285   if (official) {
1286     TH1F*  histEta  = static_cast<TH1F*>(official->FindObject("fHistEta"));
1287     if (histEta)  {
1288       TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
1289                             histEta->GetXaxis()->GetXmin(),
1290                             histEta->GetXaxis()->GetXmax());
1291       for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1292         oEta->SetBinContent(i, histEta->GetBinContent(i));
1293         oEta->SetBinError(i, histEta->GetBinError(i));
1294       }
1295       if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1296       oEta->Scale(1./ntotal, "width");
1297       oEta->SetDirectory(0);
1298       oEta->SetMarkerStyle(marker+4);
1299       oEta->SetMarkerColor(color+5);
1300       fOutput->Add(oEta);
1301       fOutput->Add(Rebin(oEta, rebin, false));
1302     }
1303     else
1304       AliWarning(Form("Couldn't find histogram fHistEta in list %s",
1305                       official->GetName()));
1306   }
1307   else
1308     AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
1309
1310   TList* tracks = static_cast<TList*>(forward->Get("tracks"));
1311   if (tracks) {
1312     TH1F*  histEta  = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
1313     if (histEta)  {
1314       TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
1315                             histEta->GetXaxis()->GetXmin(),
1316                             histEta->GetXaxis()->GetXmax());
1317       for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1318         oEta->SetBinContent(i, histEta->GetBinContent(i));
1319         oEta->SetBinError(i, histEta->GetBinError(i));
1320       }
1321       if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1322       oEta->Scale(1./ntotal, "width");
1323       oEta->SetDirectory(0);
1324       oEta->SetMarkerStyle(marker);
1325       oEta->SetMarkerColor(color+5);
1326       fOutput->Add(oEta);
1327       fOutput->Add(Rebin(oEta, rebin, false));
1328     }
1329     else
1330       AliWarning(Form("Couldn't find histogram fHistEta in list %s",
1331                       tracks->GetName()));
1332   }
1333   else
1334     AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
1335
1336   forward->Close();
1337 }
1338
1339 //
1340 // EOF
1341 //