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