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     fUseShapeCorr(true),
  33     fSNNString(0),
  34     fSysString(0),
  35     fCent(0)
  36 {
  37   //
  38   // Constructor
  39   //
  40 }
  41
  42 //____________________________________________________________________
  43 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  44   : AliAnalysisTaskSE(name),
  45     fSums(0),           // Container of sums
  46     fOutput(0),         // Container of output
  47     fVtxMin(-10),       // Minimum v_z
  48     fVtxMax(10),        // Maximum v_z
  49     fTriggerMask(AliAODForwardMult::kInel),
  50     fRebin(5),          // Rebinning factor
  51     fCutEdges(false),
  52     fSymmetrice(true),
  53     fCorrEmpty(true),
  54     fTriggerEff(1),
  55     fShapeCorr(0),
  56     fListOfCentralities(0),
  57     fUseShapeCorr(true),
  58     fSNNString(0),
  59     fSysString(0),
  60     fCent(0)
  61 {
  62   //
  63   // Constructor
  64   //
  65   fListOfCentralities = new TList;
  66
  67   // Output slot #1 writes into a TH1 container
  68   DefineOutput(1, TList::Class());
  69   DefineOutput(2, TList::Class());
  70 }
  71
  72 //____________________________________________________________________
  73 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
  74   : AliAnalysisTaskSE(o),
  75     fSums(o.fSums),             // TList* - Container of sums
  76     fOutput(o.fOutput),         // Container of output
  77     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
  78     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
  79     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
  80     fRebin(o.fRebin),           // Int_t - Rebinning factor
  81     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
  82     fSymmetrice(true),
  83     fCorrEmpty(true),
  84     fTriggerEff(o.fTriggerEff),
  85     fShapeCorr(o.fShapeCorr),
  86     fListOfCentralities(o.fListOfCentralities),
  87     fUseShapeCorr(o.fUseShapeCorr),
  88     fSNNString(o.fSNNString),
  89     fSysString(o.fSysString),
  90     fCent(o.fCent)
  91 {}
  92
  93 //____________________________________________________________________
  94 AliBasedNdetaTask::~AliBasedNdetaTask()
  95 {
  96   //
  97   // Destructor
  98   //
  99   if (fSums) {
 100     fSums->Delete();
 101     delete fSums;
 102     fSums = 0;
 103   }
 104   if (fOutput) {
 105     fOutput->Delete();
 106     delete fOutput;
 107     fOutput = 0;
 108   }
 109 }
 110
 111 //________________________________________________________________________
 112 void
 113 AliBasedNdetaTask::AddCentralityBin(Short_t low, Short_t high)
 114 {
 115   //
 116   // Add a centrality bin
 117   //
 118   // Parameters:
 119   //    low  Low cut
 120   //    high High cut
 121   //
 122   fListOfCentralities->Add(MakeCentralityBin(GetName(), low, high));
 123 }
 124
 125 //________________________________________________________________________
 126 AliBasedNdetaTask::CentralityBin*
 127 AliBasedNdetaTask::MakeCentralityBin(const char* name,
 128                                      Short_t low, Short_t high) const
 129 {
 130   //
 131   // Make a centrality bin
 132   //
 133   // Parameters:
 134   //    name  Name used for histograms
 135   //    low   Low cut in percent
 136   //    high  High cut in percent
 137   //
 138   // Return:
 139   //    A newly created centrality bin
 140   //
 141   return new CentralityBin(name, low, high);
 142 }
 143 //________________________________________________________________________
 144 void
 145 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 146 {
 147   //
 148   // Set the trigger maskl
 149   //
 150   // Parameters:
 151   //    mask Trigger mask
 152   //
 153   UShort_t    trgMask = 0;
 154   TString     trgs(mask);
 155   trgs.ToUpper();
 156   TObjString* trg;
 157   TIter       next(trgs.Tokenize(" ,|"));
 158   while ((trg = static_cast<TObjString*>(next()))) {
 159     TString s(trg->GetString());
 160     if      (s.IsNull()) continue;
 161     if      (s.CompareTo("INEL")  == 0) trgMask = AliAODForwardMult::kInel;
 162     else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
 163     else if (s.CompareTo("NSD")   == 0) trgMask = AliAODForwardMult::kNSD;
 164     else
 165       Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
 166   }
 167   if (trgMask == 0) trgMask = 1;
 168   SetTriggerMask(trgMask);
 169 }
 170
 171 //________________________________________________________________________
 172 void
 173 AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
 174 {
 175   //
 176   // Set the shape correction (a.k.a., track correction) for selected
 177   // trigger(s)
 178   //
 179   // Parameters:
 180   //    h Correction
 181   //
 182   if (!c) return;
 183   fShapeCorr = static_cast<TH1*>(c->Clone());
 184   fShapeCorr->SetDirectory(0);
 185 }
 186
 187 //________________________________________________________________________
 188 void
 189 AliBasedNdetaTask::UserCreateOutputObjects()
 190 {
 191   //
 192   // Create output objects.
 193   //
 194   // This is called once per slave process
 195   //
 196   fSums = new TList;
 197   fSums->SetName(Form("%s_sums", GetName()));
 198   fSums->SetOwner();
 199
 200   // Automatically add 'all' centrality bin if nothing has been defined.
 201   if (fListOfCentralities->GetEntries() <= 0) AddCentralityBin(0,0);
 202
 203   // Centrality histogram
 204   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
 205   fCent->SetDirectory(0);
 206   fCent->SetXTitle(0);
 207   fSums->Add(fCent);
 208
 209   // Loop over centrality bins
 210   TIter next(fListOfCentralities);
 211   CentralityBin* bin = 0;
 212   while ((bin = static_cast<CentralityBin*>(next())))
 213     bin->CreateOutputObjects(fSums);
 214
 215   // Check that we have an AOD input handler
 216   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 217   AliAODInputHandler* ah =
 218     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
 219   if (!ah) AliFatal("No AOD input handler set in analysis manager");
 220
 221   // Post data for ALL output slots >0 here, to get at least an empty histogram
 222   PostData(1, fSums);
 223 }
 224 //____________________________________________________________________
 225 void
 226 AliBasedNdetaTask::UserExec(Option_t *)
 227 {
 228   //
 229   // Process a single event
 230   //
 231   // Parameters:
 232   //    option Not used
 233   //
 234   // Main loop
 235   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 236   if (!aod) {
 237     AliError("Cannot get the AOD event");
 238     return;
 239   }
 240
 241   TObject* obj = aod->FindListObject("Forward");
 242   if (!obj) {
 243     AliWarning("No forward object found");
 244     return;
 245   }
 246   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 247
 248   // Fill centrality histogram
 249   fCent->Fill(forward->GetCentrality());
 250
 251   // Get the histogram(s)
 252   TH2D* data   = GetHistogram(aod, false);
 253   TH2D* dataMC = GetHistogram(aod, true);
 254
 255   // Loop over centrality bins
 256   TIter next(fListOfCentralities);
 257   CentralityBin* bin = 0;
 258   while ((bin = static_cast<CentralityBin*>(next())))
 259     bin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
 260
 261   // Here, we get the update
 262   if (!fSNNString) {
 263     UShort_t sNN = forward->GetSNN();
 264     fSNNString = new TNamed("sNN", "");
 265     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
 266     fSNNString->SetUniqueID(sNN);
 267     fSums->Add(fSNNString);
 268
 269     UShort_t sys = forward->GetSystem();
 270     fSysString = new TNamed("sys", "");
 271     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
 272     fSysString->SetUniqueID(sys);
 273     fSums->Add(fSysString);
 274   }
 275
 276   PostData(1, fSums);
 277 }
 278
 279 //________________________________________________________________________
 280 void
 281 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 282                                           const char* title, const char* ytitle)
 283 {
 284   //
 285   // Set histogram graphical options, etc.
 286   //
 287   // Parameters:
 288   //    h       Histogram to modify
 289   //    colour  Marker color
 290   //    marker  Marker style
 291   //    title   Title of histogram
 292   //    ytitle  Title on y-axis.
 293   //
 294   h->SetTitle(title);
 295   h->SetMarkerColor(colour);
 296   h->SetMarkerStyle(marker);
 297   h->SetMarkerSize(1);
 298   h->SetFillStyle(0);
 299   h->SetYTitle(ytitle);
 300   h->GetXaxis()->SetTitleFont(132);
 301   h->GetXaxis()->SetLabelFont(132);
 302   h->GetXaxis()->SetNdivisions(10);
 303   h->GetYaxis()->SetTitleFont(132);
 304   h->GetYaxis()->SetLabelFont(132);
 305   h->GetYaxis()->SetNdivisions(10);
 306   h->GetYaxis()->SetDecimals();
 307   h->SetStats(0);
 308 }
 309
 310 //________________________________________________________________________
 311 TH1D*
 312 AliBasedNdetaTask::ProjectX(const TH2D* h,
 313                             const char* name,
 314                             Int_t firstbin,
 315                             Int_t lastbin,
 316                             bool  corr,
 317                             bool  error)
 318 {
 319   //
 320   // Project onto the X axis
 321   //
 322   // Parameters:
 323   //    h         2D histogram
 324   //    name      New name
 325   //    firstbin  First bin to use
 326   //    lastbin   Last bin to use
 327   //    error     Whether to calculate errors
 328   //
 329   // Return:
 330   //    Newly created histogram or null
 331   //
 332   if (!h) return 0;
 333 #if USE_ROOT_PROJECT
 334   return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 335 #endif
 336
 337   TAxis* xaxis = h->GetXaxis();
 338   TAxis* yaxis = h->GetYaxis();
 339   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 340                           xaxis->GetXmin(), xaxis->GetXmax());
 341   static_cast<const TAttLine*>(h)->Copy(*ret);
 342   static_cast<const TAttFill*>(h)->Copy(*ret);
 343   static_cast<const TAttMarker*>(h)->Copy(*ret);
 344   ret->GetXaxis()->ImportAttributes(xaxis);
 345
 346   Int_t first = firstbin;
 347   Int_t last  = lastbin;
 348   if (first < 0)                         first = 0;
 349   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
 350   if      (last  < 0)                    last  = yaxis->GetNbins();
 351   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
 352   if (last-first < 0) {
 353     AliWarningGeneral("AliBasedNdetaTask",
 354                       Form("Nothing to project [%d,%d]", first, last));
 355     return 0;
 356
 357   }
 358
 359   // Loop over X bins
 360   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 361   Int_t ybins = (last-first+1);
 362   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 363     Double_t content = 0;
 364     Double_t error2  = 0;
 365     Int_t    nbins   = 0;
 366
 367
 368     for (Int_t ybin = first; ybin <= last; ybin++) {
 369       Double_t c1 = h->GetCellContent(xbin, ybin);
 370       Double_t e1 = h->GetCellError(xbin, ybin);
 371
 372       // Ignore empty bins
 373       if (c1 < 1e-12) continue;
 374       if (e1 < 1e-12) {
 375         if (error) continue;
 376         e1 = 1;
 377       }
 378
 379       content    += c1;
 380       error2     += e1*e1;
 381       nbins++;
 382     } // for (ybin)
 383     if(content > 0 && nbins > 0) {
 384       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 385       if (error) {
 386         // Calculate weighted average
 387         ret->SetBinContent(xbin, content * factor);
 388         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 389       }
 390       else
 391         ret->SetBinContent(xbin, factor * content);
 392     }
 393   } // for (xbin)
 394
 395   return ret;
 396 }
 397
 398 //________________________________________________________________________
 399 void
 400 AliBasedNdetaTask::Terminate(Option_t *)
 401 {
 402   //
 403   // Called at end of event processing..
 404   //
 405   // This is called once in the master
 406   //
 407   // Parameters:
 408   //    option Not used
 409   //
 410   // Draw result to screen, or perform fitting, normalizations Called
 411   // once at the end of the query
 412
 413   fSums = dynamic_cast<TList*> (GetOutputData(1));
 414   if(!fSums) {
 415     AliError("Could not retrieve TList fSums");
 416     return;
 417   }
 418
 419   fOutput = new TList;
 420   fOutput->SetName(Form("%s_result", GetName()));
 421   fOutput->SetOwner();
 422
 423   fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
 424   fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
 425
 426   if(fSysString && fSNNString &&
 427      fSysString->GetUniqueID() == AliForwardUtil::kPP)
 428     LoadNormalizationData(fSysString->GetUniqueID(),
 429                           fSNNString->GetUniqueID());
 430
 431   // Loop over centrality bins
 432   TIter next(fListOfCentralities);
 433   CentralityBin* bin = 0;
 434   while ((bin = static_cast<CentralityBin*>(next())))
 435     bin->End(fSums, fOutput, fUseShapeCorr ? fShapeCorr : 0, fTriggerEff,
 436              fSymmetrice, fRebin, fCorrEmpty, fCutEdges,
 437              fVtxMin, fVtxMax, fTriggerMask);
 438
 439   // Output collision energy string
 440   if (fSNNString) fOutput->Add(fSNNString->Clone());
 441
 442   // Output collision system string
 443   if (fSysString) fOutput->Add(fSysString->Clone());
 444
 445   // Output trigger string
 446   TNamed* trigString =
 447     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
 448   trigString->SetUniqueID(fTriggerMask);
 449   fOutput->Add(trigString);
 450
 451   // Output vertex axis
 452   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 453   vtxAxis->SetName("vtxAxis");
 454   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 455   fOutput->Add(vtxAxis);
 456
 457   // Output trigger efficiency and shape correction
 458   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
 459   if (fShapeCorr) fOutput->Add(fShapeCorr);
 460
 461   PostData(2, fOutput);
 462 }
 463 //________________________________________________________________________
 464 void
 465 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 466 {
 467   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
 468   TString type("pp");
 469   TString snn("900");
 470   if(energy == 7000) snn.Form("7000");
 471   if(energy == 2750) snn.Form("2750");
 472
 473   if(fShapeCorr && (fTriggerEff != 1)) {
 474     AliInfo("Objects already set for normalization - no action taken");
 475     return;
 476   }
 477
 478   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/"
 479                                 "Normalization/normalizationHists_%s_%s.root",
 480                                 type.Data(),snn.Data()));
 481   if(!fin) {
 482     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
 483     return;
 484   }
 485
 486   // Shape correction
 487   TString shapeCorName(Form("h%sNormalization",
 488                             fTriggerMask == AliAODForwardMult::kInel ? "Inel" :
 489                             fTriggerMask == AliAODForwardMult::kNSD ? "NSD" :
 490                             fTriggerMask == AliAODForwardMult::kInelGt0 ?
 491                             "InelGt0" : "All"));
 492   TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
 493   if (shapeCor) SetShapeCorrection(shapeCor);
 494
 495   // Trigger efficiency
 496   TString effName(Form("%sTriggerEff",
 497                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 498                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 499                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
 500                        "inelgt0" : "all"));
 501   TParameter<float>* eff = static_cast<TParameter<float>*>(fin->Get(effName));
 502   if (eff) SetTriggerEff(eff->GetVal());
 503
 504   // Print - out
 505   if (shapeCor && eff) AliInfo("Loaded objects for normalization.");
 506 }
 507 //________________________________________________________________________
 508 TH1D*
 509 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
 510 {
 511   //
 512   // Make a copy of the input histogram and rebin that histogram
 513   //
 514   // Parameters:
 515   //    h  Histogram to rebin
 516   //
 517   // Return:
 518   //    New (rebinned) histogram
 519   //
 520   if (rebin <= 1) return 0;
 521
 522   Int_t nBins = h->GetNbinsX();
 523   if(nBins % rebin != 0) {
 524     AliWarningGeneral("AliBasedNdetaTask",
 525                       Form("Rebin factor %d is not a devisor of current number "
 526                            "of bins %d in the histogram %s",
 527                            rebin, nBins, h->GetName()));
 528     return 0;
 529   }
 530
 531   // Make a copy
 532   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
 533                                                h->GetName(), rebin)));
 534   tmp->Rebin(rebin);
 535   tmp->SetDirectory(0);
 536
 537   // The new number of bins
 538   Int_t nBinsNew = nBins / rebin;
 539   for(Int_t i = 1;i<= nBinsNew; i++) {
 540     Double_t content = 0;
 541     Double_t sumw    = 0;
 542     Double_t wsum    = 0;
 543     Int_t    nbins   = 0;
 544     for(Int_t j = 1; j<=rebin;j++) {
 545       Int_t    bin = (i-1)*rebin + j;
 546       Double_t c   =  h->GetBinContent(bin);
 547       if (c <= 0) continue;
 548
 549       if (cutEdges) {
 550         if (h->GetBinContent(bin+1)<=0 ||
 551             h->GetBinContent(bin-1)<=0) {
 552 #if 0
 553           AliWarningGeneral("AliBasedNdetaTask",
 554                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
 555                                  bin, c, h->GetName(),
 556                                  bin+1, h->GetBinContent(bin+1),
 557                                  bin-1, h->GetBinContent(bin-1)));
 558 #endif
 559           continue;
 560         }
 561       }
 562       Double_t e =  h->GetBinError(bin);
 563       Double_t w =  1 / (e*e); // 1/c/c
 564       content    += c;
 565       sumw       += w;
 566       wsum       += w * c;
 567       nbins++;
 568     }
 569
 570     if(content > 0 && nbins > 0) {
 571       tmp->SetBinContent(i, wsum / sumw);
 572       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 573     }
 574   }
 575
 576   return tmp;
 577 }
 578
 579 //__________________________________________________________________
 580 TH1*
 581 AliBasedNdetaTask::Symmetrice(const TH1* h)
 582 {
 583   //
 584   // Make an extension of @a h to make it symmetric about 0
 585   //
 586   // Parameters:
 587   //    h Histogram to symmertrice
 588   //
 589   // Return:
 590   //    Symmetric extension of @a h
 591   //
 592   Int_t nBins = h->GetNbinsX();
 593   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 594   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 595   s->Reset();
 596   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 597   s->SetMarkerColor(h->GetMarkerColor());
 598   s->SetMarkerSize(h->GetMarkerSize());
 599   s->SetMarkerStyle(h->GetMarkerStyle()+4);
 600   s->SetFillColor(h->GetFillColor());
 601   s->SetFillStyle(h->GetFillStyle());
 602   s->SetDirectory(0);
 603
 604   // Find the first and last bin with data
 605   Int_t first = nBins+1;
 606   Int_t last  = 0;
 607   for (Int_t i = 1; i <= nBins; i++) {
 608     if (h->GetBinContent(i) <= 0) continue;
 609     first = TMath::Min(first, i);
 610     last  = TMath::Max(last,  i);
 611   }
 612
 613   Double_t xfirst = h->GetBinCenter(first-1);
 614   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 615   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 616   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 617     s->SetBinContent(j, h->GetBinContent(i));
 618     s->SetBinError(j, h->GetBinError(i));
 619   }
 620   // Fill in overlap bin
 621   s->SetBinContent(l2+1, h->GetBinContent(first));
 622   s->SetBinError(l2+1, h->GetBinError(first));
 623   return s;
 624 }
 625
 626 //====================================================================
 627 AliBasedNdetaTask::CentralityBin::CentralityBin()
 628   : TNamed("", ""),
 629     fSums(0),
 630     fOutput(0),
 631     fSum(0),
 632     fSumMC(0),
 633     fTriggers(0),
 634     fLow(0),
 635     fHigh(0)
 636 {
 637   //
 638   // Constructor
 639   //
 640 }
 641 //____________________________________________________________________
 642 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
 643                                                 Short_t low, Short_t high)
 644   : TNamed(name, ""),
 645     fSums(0),
 646     fOutput(0),
 647     fSum(0),
 648     fSumMC(0),
 649     fTriggers(0),
 650     fLow(low),
 651     fHigh(high)
 652 {
 653   //
 654   // Constructor
 655   //
 656   // Parameters:
 657   //    name Name used for histograms (e.g., Forward)
 658   //    low  Lower centrality cut in percent
 659   //    high Upper centrality cut in percent
 660   //
 661   if (low <= 0 && high <= 0) {
 662     fLow  = 0;
 663     fHigh = 0;
 664     SetTitle("All centralities");
 665   }
 666   else {
 667     fLow  = low;
 668     fHigh = high;
 669     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
 670   }
 671 }
 672 //____________________________________________________________________
 673 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
 674   : TNamed(o),
 675     fSums(o.fSums),
 676     fOutput(o.fOutput),
 677     fSum(o.fSum),
 678     fSumMC(o.fSumMC),
 679     fTriggers(o.fTriggers),
 680     fLow(o.fLow),
 681     fHigh(o.fHigh)
 682 {
 683   //
 684   // Copy constructor
 685   //
 686   // Parameters:
 687   //    other Object to copy from
 688   //
 689 }
 690 //____________________________________________________________________
 691 AliBasedNdetaTask::CentralityBin::~CentralityBin()
 692 {
 693   //
 694   // Destructor
 695   //
 696   if (fSums) fSums->Delete();
 697   if (fOutput) fOutput->Delete();
 698 }
 699
 700 //____________________________________________________________________
 701 AliBasedNdetaTask::CentralityBin&
 702 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
 703 {
 704   //
 705   // Assignment operator
 706   //
 707   // Parameters:
 708   //    other Object to assign from
 709   //
 710   // Return:
 711   //    Reference to this
 712   //
 713   SetName(o.GetName());
 714   SetTitle(o.GetTitle());
 715   fSums      = o.fSums;
 716   fOutput    = o.fOutput;
 717   fSum       = o.fSum;
 718   fSumMC     = o.fSumMC;
 719   fTriggers  = o.fTriggers;
 720   fLow       = o.fLow;
 721   fHigh      = o.fHigh;
 722
 723   return *this;
 724 }
 725 //____________________________________________________________________
 726 const char*
 727 AliBasedNdetaTask::CentralityBin::GetListName() const
 728 {
 729   //
 730   // Get the list name
 731   //
 732   // Return:
 733   //    List Name
 734   //
 735   if (IsAllBin()) return "all";
 736   return Form("cent%03d_%03d", fLow, fHigh);
 737 }
 738 //____________________________________________________________________
 739 void
 740 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
 741 {
 742   //
 743   // Create output objects
 744   //
 745   // Parameters:
 746   //    dir   Parent list
 747   //
 748   fSums = new TList;
 749   fSums->SetName(GetListName());
 750   fSums->SetOwner();
 751   dir->Add(fSums);
 752
 753   fTriggers = new TH1D("triggers", "Number of triggers",
 754                        kAccepted, 1, kAccepted);
 755   fTriggers->SetYTitle("# of events");
 756   fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
 757   fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
 758   fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
 759   fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
 760   fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
 761   fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
 762   fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
 763   fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
 764   fTriggers->GetXaxis()->SetBinLabel(kPileUp,      "w/Pileup");
 765   fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
 766   fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
 767   fTriggers->SetFillColor(kRed+1);
 768   fTriggers->SetFillStyle(3001);
 769   fTriggers->SetStats(0);
 770   fSums->Add(fTriggers);
 771 }
 772 //____________________________________________________________________
 773 void
 774 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
 775 {
 776   //
 777   // Create sum histogram
 778   //
 779   // Parameters:
 780   //    data  Data histogram to clone
 781   //    mc    (optional) MC histogram to clone
 782   //
 783   fSum = static_cast<TH2D*>(data->Clone(GetName()));
 784   fSum->SetDirectory(0);
 785   fSum->Reset();
 786   fSums->Add(fSum);
 787
 788   // If no MC data is given, then do not create MC sum histogram
 789   if (!mc) return;
 790
 791   fSumMC = static_cast<TH2D*>(mc->Clone(Form("%sMC", GetName())));
 792   fSumMC->SetDirectory(0);
 793   fSumMC->Reset();
 794   fSums->Add(fSumMC);
 795 }
 796
 797 //____________________________________________________________________
 798 Bool_t
 799 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
 800                                              Int_t triggerMask,
 801                                              Double_t vzMin, Double_t vzMax)
 802 {
 803   //
 804   // Check the trigger, vertex, and centrality
 805   //
 806   // Parameters:
 807   //    aod Event input
 808   //
 809   // Return:
 810   //    true if the event is to be used
 811   //
 812   if (!forward) return false;
 813
 814   // Check the centrality class unless this is the 'all' bin
 815   if (!IsAllBin()) {
 816     Double_t centrality = forward->GetCentrality();
 817     if (centrality < fLow || centrality >= fHigh) return false;
 818   }
 819
 820   fTriggers->AddBinContent(kAll);
 821   if (forward->IsTriggerBits(AliAODForwardMult::kB))
 822     fTriggers->AddBinContent(kB);
 823   if (forward->IsTriggerBits(AliAODForwardMult::kA))
 824     fTriggers->AddBinContent(kA);
 825   if (forward->IsTriggerBits(AliAODForwardMult::kC))
 826     fTriggers->AddBinContent(kC);
 827   if (forward->IsTriggerBits(AliAODForwardMult::kE))
 828     fTriggers->AddBinContent(kE);
 829   if (forward->IsTriggerBits(AliAODForwardMult::kInel))
 830     fTriggers->AddBinContent(kMB);
 831   if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
 832     fTriggers->AddBinContent(kMCNSD);
 833   if (forward->IsTriggerBits(AliAODForwardMult::kPileUp))
 834     fTriggers->AddBinContent(kPileUp);
 835
 836   // Check if we have an event of interest.
 837   if (!forward->IsTriggerBits(triggerMask)) return false;
 838
 839   //Check for pileup
 840   if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return false;
 841
 842   fTriggers->AddBinContent(kWithTrigger);
 843
 844   // Check that we have a valid vertex
 845   if (!forward->HasIpZ()) return false;
 846   fTriggers->AddBinContent(kWithVertex);
 847
 848   // Check that vertex is within cuts
 849   if (!forward->InRange(vzMin, vzMax)) return false;
 850   fTriggers->AddBinContent(kAccepted);
 851
 852   return true;
 853 }
 854
 855
 856 //____________________________________________________________________
 857 void
 858 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
 859                                                Int_t triggerMask,
 860                                                Double_t vzMin, Double_t vzMax,
 861                                                const TH2D* data, const TH2D* mc)
 862 {
 863   //
 864   // Process an event
 865   //
 866   // Parameters:
 867   //    forward     Forward data (for trigger, vertex, & centrality)
 868   //    triggerMask Trigger mask
 869   //    vzMin       Minimum IP z coordinate
 870   //    vzMax       Maximum IP z coordinate
 871   //    data        Data histogram
 872   //    mc          MC histogram
 873   //
 874   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
 875   if (!data) return;
 876   if (!fSum) CreateSums(data, mc);
 877
 878   fSum->Add(data);
 879   if (mc) fSumMC->Add(mc);
 880 }
 881
 882 //________________________________________________________________________
 883 void
 884 AliBasedNdetaTask::CentralityBin::End(TList*      sums,
 885                                       TList*      results,
 886                                       const TH1*  shapeCorr,
 887                                       Double_t    trigEff,
 888                                       Bool_t      symmetrice,
 889                                       Int_t       rebin,
 890                                       Bool_t      corrEmpty,
 891                                       Bool_t      cutEdges,
 892                                       Double_t    vzMin,
 893                                       Double_t    vzMax,
 894                                       Int_t       triggerMask)
 895 {
 896   //
 897   // End of processing
 898   //
 899   // Parameters:
 900   //    sums        List of sums
 901   //    results     Output list of results
 902   //    shapeCorr   Shape correction or nil
 903   //    trigEff     Trigger efficiency
 904   //    symmetrice  Whether to symmetrice the results
 905   //    rebin       Whether to rebin the results
 906   //    corrEmpty   Whether to correct for empty bins
 907   //    cutEdges    Whether to cut edges when rebinning
 908   //    vzMin       Minimum IP z coordinate
 909   //    vzMax     Maximum IP z coordinate
 910   //    triggerMask Trigger mask
 911   //
 912
 913   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
 914   if(!fSums) {
 915     AliError("Could not retrieve TList fSums");
 916     return;
 917   }
 918
 919   fOutput = new TList;
 920   fOutput->SetName(GetListName());
 921   fOutput->SetOwner();
 922   results->Add(fOutput);
 923
 924   fSum      = static_cast<TH2D*>(fSums->FindObject(GetName()));
 925   fSumMC    = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
 926   fTriggers = static_cast<TH1D*>(fSums->FindObject("triggers"));
 927
 928   if (!fTriggers) {
 929     AliError("Couldn't find histogram 'triggers' in list");
 930     return;
 931   }
 932   if (!fSum) {
 933     AliError("Couldn't find histogram 'base' in list");
 934     return;
 935   }
 936
 937   Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
 938   Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
 939   Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
 940   Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
 941   Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
 942   Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
 943   Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
 944   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
 945   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
 946   Int_t nGood       = nB - nA - nC + 2 * nE;
 947
 948   Double_t alpha    = Double_t(nAccepted) / nWithVertex;
 949   Double_t vNorm     = nAccepted + alpha*(nTriggered - nWithVertex);
 950   Double_t vtxEff   = Double_t(nAccepted) / vNorm;
 951   vtxEff            = vtxEff * trigEff;
 952
 953   if (nTriggered <= 0) {
 954     AliError("Number of triggered events <= 0");
 955     return;
 956   }
 957   if (nGood <= 0) {
 958     AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
 959                     nGood, nB, nA, nC, nE));
 960     nGood = nMB;
 961   }
 962   AliInfo(Form("Total of %9d events for %s\n"
 963                "                   of these %9d are minimum bias\n"
 964                "                   of these %9d has a %s trigger\n"
 965                "                   of these %9d has a vertex\n"
 966                "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
 967                "                   Triggers by type:\n"
 968                "                     B   = %9d\n"
 969                "                     A|C = %9d (%9d+%-9d)\n"
 970                "                     E   = %9d\n"
 971                "                   Implies %9d good triggers\n"
 972                "                   Vertex efficiency: %f (%f)",
 973                nAll, GetTitle(), nMB, nTriggered,
 974                AliAODForwardMult::GetTriggerString(triggerMask),
 975                nWithVertex, nAccepted,
 976                vzMin, vzMax,
 977                nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
 978
 979   const char* name = GetName();
 980
 981   // Get acceptance normalisation from underflow bins
 982   TH1D* norm = ProjectX(fSum, Form("norm%s",name), 0, 0, corrEmpty, false);
 983   norm->SetDirectory(0);
 984
 985   // Scale by shape correction
 986   if(shapeCorr) fSum->Divide(shapeCorr);
 987   else AliInfo("No shape correction specified, or disabled");
 988
 989   // Project on X axis
 990   TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
 991                           corrEmpty);
 992   dndeta->SetDirectory(0);
 993
 994   // Normalize to the acceptance -
 995   dndeta->Divide(norm);
 996
 997   // Scale by the vertex efficiency
 998   dndeta->Scale(vtxEff, "width");
 999
1000   SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
1001   SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name));
1002
1003   fOutput->Add(fTriggers->Clone());
1004   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
1005   fOutput->Add(dndeta);
1006   fOutput->Add(norm);
1007   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1008   if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1009
1010   if (fSumMC) {
1011     // Get acceptance normalisation from underflow bins
1012     TH1D* normMC   = ProjectX(fSumMC,Form("norm%sMC", name), 0, 1,
1013                               corrEmpty, false);
1014     // Project onto eta axis - _ignoring_underflow_bins_!
1015     TH1D* dndetaMC = ProjectX(fSumMC,Form("dndeta%sMC", name),1,
1016                               fSum->GetNbinsY(), corrEmpty);
1017     // Normalize to the acceptance
1018     dndetaMC->Divide(normMC);
1019     // Scale by the vertex efficiency
1020     dndetaMC->Scale(vNorm, "width");
1021     normMC->Scale(1. / nAccepted);
1022
1023     SetHistogramAttributes(dndetaMC, kRed+3, 21, Form("ALICE %s (MC)",name));
1024     SetHistogramAttributes(normMC, kRed+3, 21, Form("ALICE %s (MC) "
1025                                                     "normalisation",name));
1026
1027     fOutput->Add(dndetaMC);
1028     if (symmetrice)   fOutput->Add(Symmetrice(dndetaMC));
1029
1030     fOutput->Add(normMC);
1031     fOutput->Add(Rebin(dndetaMC, rebin, cutEdges));
1032
1033     if (symmetrice)
1034       fOutput->Add(Symmetrice(Rebin(dndetaMC, rebin, cutEdges)));
1035   }
1036 }
1037
1038 //
1039 // EOF
1040 //