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