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