PWGLF/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 <AliAnalysisManager.h>
   9 #include <AliAODEvent.h>
  10 #include <AliAODHandler.h>
  11 #include <AliAODInputHandler.h>
  12 #include "AliForwardUtil.h"
  13 #include "AliAODForwardMult.h"
  14 #include <TFile.h>
  15 #include <TStyle.h>
  16
  17 //____________________________________________________________________
  18 AliBasedNdetaTask::AliBasedNdetaTask()
  19   : AliAnalysisTaskSE(),
  20     fSums(0),           // Container of sums
  21     fOutput(0),         // Container of output
  22     fVtxMin(0),         // Minimum v_z
  23     fVtxMax(0),         // Maximum v_z
  24     fTriggerMask(0),    // Trigger mask
  25     fRebin(0),          // Rebinning factor
  26     fCutEdges(false),
  27     fSymmetrice(true),
  28     fCorrEmpty(true),
  29     fUseROOTProj(false),
  30     fTriggerEff(1),
  31     fTriggerEff0(1),
  32     fShapeCorr(0),
  33     fListOfCentralities(0),
  34     fSNNString(0),
  35     fSysString(0),
  36     fCent(0),
  37     fCentAxis(0),
  38     fVtx(0),
  39     fNormalizationScheme(kFull),
  40     fSchemeString(0),
  41     fTriggerString(0),
  42     fFinalMCCorrFile(""),
  43     fSatelliteVertices(0),
  44     fglobalempiricalcorrection(0),
  45     fmeabsignalvscentr(0)
  46 {
  47   //
  48   // Constructor
  49   //
  50   DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask");
  51 }
  52
  53 //____________________________________________________________________
  54 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  55   : AliAnalysisTaskSE(name),
  56     fSums(0),           // Container of sums
  57     fOutput(0),         // Container of output
  58     fVtxMin(-10),       // Minimum v_z
  59     fVtxMax(10),        // Maximum v_z
  60     fTriggerMask(AliAODForwardMult::kInel),
  61     fRebin(5),          // Rebinning factor
  62     fCutEdges(false),
  63     fSymmetrice(true),
  64     fCorrEmpty(true),
  65     fUseROOTProj(false),
  66     fTriggerEff(1),
  67     fTriggerEff0(1),
  68     fShapeCorr(0),
  69     fListOfCentralities(0),
  70     fSNNString(0),
  71     fSysString(0),
  72     fCent(0),
  73     fCentAxis(0),
  74     fVtx(0),
  75     fNormalizationScheme(kFull),
  76     fSchemeString(0),
  77     fTriggerString(0),
  78     fFinalMCCorrFile(""),
  79     fSatelliteVertices(0),
  80     fglobalempiricalcorrection(0),
  81     fmeabsignalvscentr(0)
  82 {
  83   //
  84   // Constructor
  85   //
  86   DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name);
  87   fListOfCentralities = new TObjArray(1);
  88
  89   // Set the normalisation scheme
  90   SetNormalizationScheme(kFull);
  91
  92   // Set the trigger mask
  93   SetTriggerMask(AliAODForwardMult::kInel);
  94
  95   // Output slot #1 writes into a TH1 container
  96   DefineOutput(1, TList::Class());
  97   DefineOutput(2, TList::Class());
  98 }
  99
 100 //____________________________________________________________________
 101 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
 102   : AliAnalysisTaskSE(o),
 103     fSums(o.fSums),             // TList* - Container of sums
 104     fOutput(o.fOutput),         // Container of output
 105     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
 106     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
 107     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
 108     fRebin(o.fRebin),           // Int_t - Rebinning factor
 109     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
 110     fSymmetrice(o.fSymmetrice),
 111     fCorrEmpty(o.fCorrEmpty),
 112     fUseROOTProj(o.fUseROOTProj),
 113     fTriggerEff(o.fTriggerEff),
 114     fTriggerEff0(o.fTriggerEff0),
 115     fShapeCorr(o.fShapeCorr),
 116     fListOfCentralities(o.fListOfCentralities),
 117     fSNNString(o.fSNNString),
 118     fSysString(o.fSysString),
 119     fCent(o.fCent),
 120     fCentAxis(o.fCentAxis),
 121     fVtx(o.fVtx),
 122     fNormalizationScheme(o.fNormalizationScheme),
 123     fSchemeString(o.fSchemeString),
 124     fTriggerString(o.fTriggerString),
 125     fFinalMCCorrFile(o.fFinalMCCorrFile),
 126     fSatelliteVertices(o.fSatelliteVertices),
 127     fglobalempiricalcorrection(o.fglobalempiricalcorrection),
 128   fmeabsignalvscentr(o.fmeabsignalvscentr)
 129 {
 130   DGUARD(fDebug, 3,"Copy CTOR of AliBasedNdetaTask");
 131 }
 132
 133 //____________________________________________________________________
 134 AliBasedNdetaTask::~AliBasedNdetaTask()
 135 {
 136   //
 137   // Destructor
 138   //
 139   DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
 140   if (fSums) {
 141     fSums->Delete();
 142     delete fSums;
 143     fSums = 0;
 144   }
 145   if (fOutput) {
 146     fOutput->Delete();
 147     delete fOutput;
 148     fOutput = 0;
 149   }
 150 }
 151
 152 //________________________________________________________________________
 153 void
 154 AliBasedNdetaTask::SetDebugLevel(Int_t lvl)
 155 {
 156   AliAnalysisTaskSE::SetDebugLevel(lvl);
 157   for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) {
 158     CentralityBin* bin =
 159       static_cast<CentralityBin*>(fListOfCentralities->At(i));
 160     bin->SetDebugLevel(lvl);
 161   }
 162 }
 163
 164 //________________________________________________________________________
 165 void
 166 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
 167 {
 168   DGUARD(fDebug,3,"Set centrality axis, %d bins", n);
 169   if (!fCentAxis) {
 170     fCentAxis = new TAxis();
 171     fCentAxis->SetName("centAxis");
 172     fCentAxis->SetTitle("Centrality [%]");
 173   }
 174   TArrayD dbins(n+1);
 175   for (UShort_t i = 0; i <= n; i++)
 176     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
 177   fCentAxis->Set(n, dbins.GetArray());
 178 }
 179
 180 //________________________________________________________________________
 181 void
 182 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
 183 {
 184   //
 185   // Add a centrality bin
 186   //
 187   // Parameters:
 188   //    low  Low cut
 189   //    high High cut
 190   //
 191   DGUARD(fDebug,3,"Add a centrality bin [%d,%d] @ %d", low, high, at);
 192   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
 193   if (!bin) {
 194     Error("AddCentralityBin",
 195           "Failed to create centrality bin for %s [%d,%d] @ %d",
 196           GetName(), low, high, at);
 197     return;
 198   }
 199   bin->SetSatelliteVertices(fSatelliteVertices);
 200   bin->SetDebugLevel(fDebug);
 201   fListOfCentralities->AddAtAndExpand(bin, at);
 202 }
 203
 204 //________________________________________________________________________
 205 AliBasedNdetaTask::CentralityBin*
 206 AliBasedNdetaTask::MakeCentralityBin(const char* name,
 207                                      Short_t low, Short_t high) const
 208 {
 209   //
 210   // Make a centrality bin
 211   //
 212   // Parameters:
 213   //    name  Name used for histograms
 214   //    low   Low cut in percent
 215   //    high  High cut in percent
 216   //
 217   // Return:
 218   //    A newly created centrality bin
 219   //
 220   DGUARD(fDebug,3,"Make a centrality bin %s [%d,%d]", name, low, high);
 221   return new CentralityBin(name, low, high);
 222 }
 223
 224 #define TESTAPPEND(SCHEME,BIT,STRING) \
 225   do { if (!(SCHEME & BIT)) break;                                      \
 226     if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false)
 227
 228 //________________________________________________________________________
 229 const Char_t*
 230 AliBasedNdetaTask::NormalizationSchemeString(UShort_t scheme)
 231 {
 232   // Create a string from normalization scheme bits
 233   static TString s;
 234   s = "";
 235
 236   if (scheme == kNone)
 237     return s.Data();
 238   if (scheme == kFull) {
 239     s = "FULL";
 240     return s.Data();
 241   }
 242   TESTAPPEND(scheme, kEventLevel,        "EVENT");
 243   TESTAPPEND(scheme, kShape,             "SHAPE");
 244   TESTAPPEND(scheme, kBackground,        "BACKGROUND");
 245   TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER");
 246   TESTAPPEND(scheme, kZeroBin,           "ZEROBIN");
 247
 248   return s.Data();
 249 }
 250 //________________________________________________________________________
 251 UShort_t
 252 AliBasedNdetaTask::ParseNormalizationScheme(const char* what)
 253 {
 254   UShort_t    scheme = 0;
 255   TString     twhat(what);
 256   twhat.ToUpper();
 257   TObjString* opt;
 258   TObjArray* token = twhat.Tokenize(" ,|");
 259   TIter       next(token);
 260   while ((opt = static_cast<TObjString*>(next()))) {
 261     TString s(opt->GetString());
 262     if      (s.IsNull()) continue;
 263     Bool_t add = true;
 264     switch (s[0]) {
 265     case '-': add = false; // Fall through
 266     case '+': s.Remove(0,1);  // Remove character
 267     }
 268     UShort_t bit = 0;
 269     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
 270     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
 271     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
 272     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
 273     else if (s.CompareTo("FULL")      == 0) bit = kFull;
 274     else if (s.CompareTo("NONE")      == 0) bit = kNone;
 275     else if (s.CompareTo("ZEROBIN")   == 0) bit = kZeroBin;
 276     else
 277       ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
 278     if (add) scheme |= bit;
 279     else     scheme ^= bit;
 280   }
 281   delete token;
 282   return scheme;
 283 }
 284 //________________________________________________________________________
 285 void
 286 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
 287 {
 288   //
 289   // Set normalisation scheme
 290   //
 291   DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
 292   SetNormalizationScheme(ParseNormalizationScheme(what));
 293 }
 294 //________________________________________________________________________
 295 void
 296 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
 297 {
 298   DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
 299   fNormalizationScheme = scheme;
 300   if (fSchemeString) delete fSchemeString;
 301   fSchemeString = AliForwardUtil::MakeParameter("scheme", scheme);
 302 }
 303 //________________________________________________________________________
 304 void
 305 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 306 {
 307   //
 308   // Set the trigger maskl
 309   //
 310   // Parameters:
 311   //    mask Trigger mask
 312   //
 313   DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
 314   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
 315 }
 316 //________________________________________________________________________
 317 void
 318 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
 319 {
 320   DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
 321   fTriggerMask = mask;
 322   if (fTriggerString) delete fTriggerString;
 323   fTriggerString = AliForwardUtil::MakeParameter("trigger", fTriggerMask);
 324 }
 325
 326 //________________________________________________________________________
 327 void
 328 AliBasedNdetaTask::SetShapeCorrection(const TH2F* c)
 329 {
 330   //
 331   // Set the shape correction (a.k.a., track correction) for selected
 332   // trigger(s)
 333   //
 334   // Parameters:
 335   //    h Correction
 336   //
 337   DGUARD(fDebug,3,"Set the shape correction: %p", c);
 338   if (!c) return;
 339   fShapeCorr = static_cast<TH2F*>(c->Clone());
 340   fShapeCorr->SetDirectory(0);
 341 }
 342 //________________________________________________________________________
 343 void
 344 AliBasedNdetaTask::InitializeCentBins()
 345 {
 346   if (fListOfCentralities->GetEntries() > 0) return;
 347
 348   // Automatically add 'all' centrality bin if nothing has been defined.
 349   AddCentralityBin(0, 0, 0);
 350   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
 351     const TArrayD* bins = fCentAxis->GetXbins();
 352     Int_t          nbin = fCentAxis->GetNbins();
 353     for (Int_t i = 0; i < nbin; i++)
 354       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
 355   }
 356 }
 357
 358 //________________________________________________________________________
 359 void
 360 AliBasedNdetaTask::UserCreateOutputObjects()
 361 {
 362   //
 363   // Create output objects.
 364   //
 365   // This is called once per slave process
 366   //
 367   DGUARD(fDebug,1,"Create user ouput object");
 368   fSums = new TList;
 369   fSums->SetName(Form("%s_sums", GetName()));
 370   fSums->SetOwner();
 371
 372   InitializeCentBins();
 373   if (fCentAxis) fSums->Add(fCentAxis);
 374
 375   fSums->Add(AliForwardUtil::MakeParameter("alirootRev",
 376                                            AliForwardUtil::AliROOTRevision()));
 377   fSums->Add(AliForwardUtil::MakeParameter("alirootBranch",
 378                                            AliForwardUtil::AliROOTBranch()));
 379   fSums->Add(AliForwardUtil::MakeParameter("empirical",
 380                                            fglobalempiricalcorrection != 0));
 381   // Centrality histogram
 382   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
 383   fCent->SetDirectory(0);
 384   fCent->SetXTitle(0);
 385   fSums->Add(fCent);
 386
 387   // Custom vertex axis that will include satellite vertices
 388   // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10
 389   // Nominal vertices are usually in -10 to 10 and we should have
 390   // 10 bins in that range.  That gives us a total of
 391   //
 392   //   10+10+10=30 bins
 393   //
 394   // or 31 bin boundaries
 395   const Int_t   nCenter = 20; // Max possible is 56
 396   const Int_t   nSat    = 10;
 397   const Int_t   nBins   = 2*nSat + nCenter;
 398   Float_t       bins[nBins+1];
 399   for (Int_t i = 0; i < nSat;i++) {
 400     bins[i] = (i-nSat-.5) * 37.5;
 401     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 402   }
 403   for (Int_t i = nSat; i < nSat+nCenter+1; i++) {
 404     bins[i] = -nCenter + (i-nSat) * 2;
 405     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 406   }
 407   for (Int_t i = nSat+nCenter+1; i < 2*nSat+nCenter; i++) {
 408     bins[i] = (i-nSat-nCenter +.5) * 37.5;
 409     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 410   }
 411   bins[nBins] = (nSat + .5) * 37.5;
 412
 413   fVtx           = new TH1D("vtx", "Vertex dist", nBins, bins);
 414   fVtx->SetDirectory(0);
 415   fVtx->SetFillColor(kRed+1);
 416   fVtx->SetFillStyle(3001);
 417   fVtx->SetXTitle("IP_{z} [cm]");
 418   fVtx->SetYTitle("Events");
 419   fSums->Add(fVtx);
 420 #if 0
 421   TAxis* a = fVtx->GetXaxis();
 422   for (Int_t i = 1; i <= nBins; i++)
 423     printf("%2d/%2d: %+6.2f - %+6.2f: %+6.2f\n",
 424            i,nBins,a->GetBinLowEdge(i),a->GetBinUpEdge(i),a->GetBinCenter(i));
 425 #endif
 426
 427   // Loop over centrality bins
 428   TIter next(fListOfCentralities);
 429   CentralityBin* bin = 0;
 430   while ((bin = static_cast<CentralityBin*>(next())))
 431     bin->CreateOutputObjects(fSums, fTriggerMask);
 432
 433   fmeabsignalvscentr=new TH2D("meanAbsSignalVsCentr",
 434                               "Mean absolute signal versus centrality",
 435                               400, 0, 20, 100, 0, 100);
 436   fSums->Add(fmeabsignalvscentr);
 437
 438   Print();
 439   // Post data for ALL output slots >0 here, to get at least an empty
 440   // histogram
 441   PostData(1, fSums);
 442 }
 443
 444 //____________________________________________________________________
 445 void
 446 AliBasedNdetaTask::UserExec(Option_t *)
 447 {
 448   //
 449   // Process a single event
 450   //
 451   // Parameters:
 452   //    option Not used
 453   //
 454   // Main loop
 455   DGUARD(fDebug,1,"Analyse the AOD event");
 456   AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
 457   if (!aod) return;
 458
 459
 460   TObject* obj = aod->FindListObject("Forward");
 461   if (!obj) {
 462     AliWarning("No forward object found");
 463     return;
 464   }
 465   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 466
 467   // Fill centrality histogram
 468
 469   Double_t vtx = forward->GetIpZ();
 470   Float_t cent = forward->GetCentrality();
 471
 472   // Get the histogram(s)
 473   TH2D* data   = GetHistogram(aod, false);
 474   TH2D* dataMC = GetHistogram(aod, true);
 475   if (!data) return;
 476
 477   CheckEventData(vtx, data, dataMC);
 478
 479   if (!ApplyEmpiricalCorrection(forward,data))
 480     return;
 481
 482 #if 0
 483   // Code disabled - breaks execution
 484   Int_t notemptybins=0;
 485   Double_t sum=0.0;
 486   for (Int_t ix=1;ix<=data->GetXaxis()->GetNbins();ix++)
 487     {
 488       Double_t sumy=0.0;
 489       for(Int_t iy=1;iy<=data->GetYaxis()->GetNbins();iy++)
 490         {
 491           if(data->GetBinContent(ix,iy)>0.0)
 492             {
 493               sumy+=data->GetBinContent(ix,iy);
 494               notemptybins++;
 495             }
 496
 497         }
 498       sum+=sumy;
 499     }
 500
 501   if(notemptybins>0)
 502     {
 503       sum=sum/((Double_t)notemptybins);
 504     }
 505   else
 506     sum=-1.0;
 507   fmeabsignalvscentr->Fill(sum,cent);
 508 #endif
 509
 510   Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
 511                    !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
 512   Bool_t taken = false;
 513
 514   // Loop over centrality bins
 515   CentralityBin* allBin =
 516     static_cast<CentralityBin*>(fListOfCentralities->At(0));
 517   if (allBin->ProcessEvent(forward, fTriggerMask, isZero,
 518                            fVtxMin, fVtxMax, data, dataMC)) taken = true;
 519
 520   // Find this centrality bin
 521   if (fCentAxis && fCentAxis->GetNbins() > 0) {
 522     Int_t          icent   = fCentAxis->FindBin(cent);
 523     CentralityBin* thisBin = 0;
 524     if (icent >= 1 && icent <= fCentAxis->GetNbins())
 525       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
 526     if (thisBin)
 527       if (thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin,
 528                                 fVtxMax, data, dataMC)) taken = true;
 529   }
 530
 531   // Fill diagnostics only if we took this event
 532   if (taken) {
 533     fCent->Fill(cent);
 534     fVtx->Fill(vtx);
 535   }
 536
 537   // Here, we get the update
 538   if (!fSNNString) {
 539     fSNNString = AliForwardUtil::MakeParameter("sNN", forward->GetSNN());
 540     fSysString = AliForwardUtil::MakeParameter("sys", forward->GetSystem());
 541
 542     fSums->Add(fSNNString);
 543     fSums->Add(fSysString);
 544     fSums->Add(fSchemeString);
 545     fSums->Add(fTriggerString);
 546
 547     // Show stuff on first event
 548     Print();
 549   }
 550   PostData(1, fSums);
 551 }
 552
 553 //________________________________________________________________________
 554 void AliBasedNdetaTask::CheckEventData(Double_t,
 555                                        TH2*,
 556                                        TH2*)
 557 {
 558 }
 559
 560 //________________________________________________________________________
 561 void
 562 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 563                                           const char* title, const char* ytitle)
 564 {
 565   //
 566   // Set histogram graphical options, etc.
 567   //
 568   // Parameters:
 569   //    h       Histogram to modify
 570   //    colour  Marker color
 571   //    marker  Marker style
 572   //    title   Title of histogram
 573   //    ytitle  Title on y-axis.
 574   //
 575   h->SetTitle(title);
 576   h->SetMarkerColor(colour);
 577   h->SetMarkerStyle(marker);
 578   h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
 579   h->SetFillStyle(0);
 580   h->SetYTitle(ytitle);
 581   h->GetXaxis()->SetTitleFont(132);
 582   h->GetXaxis()->SetLabelFont(132);
 583   h->GetXaxis()->SetNdivisions(10);
 584   h->GetYaxis()->SetTitleFont(132);
 585   h->GetYaxis()->SetLabelFont(132);
 586   h->GetYaxis()->SetNdivisions(10);
 587   h->GetYaxis()->SetDecimals();
 588   h->SetStats(0);
 589 }
 590
 591 //________________________________________________________________________
 592 void
 593 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm)
 594 {
 595   // Normalize to the acceptance -
 596   // dndeta->Divide(accNorm);
 597   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 598     Double_t a = norm->GetBinContent(i);
 599     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
 600       if (a <= 0) {
 601         copy->SetBinContent(i,j,0);
 602         copy->SetBinError(i,j,0);
 603         continue;
 604       }
 605       Double_t c = copy->GetBinContent(i, j);
 606       Double_t e = copy->GetBinError(i, j);
 607       copy->SetBinContent(i, j, c / a);
 608       copy->SetBinError(i, j, e / a);
 609     }
 610   }
 611 }
 612 //________________________________________________________________________
 613 void
 614 AliBasedNdetaTask::ScaleToCoverage(TH1D* copy, const TH1D* norm)
 615 {
 616   // Normalize to the acceptance -
 617   // dndeta->Divide(accNorm);
 618   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 619     Double_t a = norm->GetBinContent(i);
 620     if (a <= 0) {
 621       copy->SetBinContent(i,0);
 622       copy->SetBinError(i,0);
 623       continue;
 624     }
 625     Double_t c = copy->GetBinContent(i);
 626     Double_t e = copy->GetBinError(i);
 627     copy->SetBinContent(i, c / a);
 628     copy->SetBinError(i, e / a);
 629   }
 630 }
 631
 632 //________________________________________________________________________
 633 TH1D*
 634 AliBasedNdetaTask::ProjectX(const TH2D* h,
 635                             const char* name,
 636                             Int_t firstbin,
 637                             Int_t lastbin,
 638                             bool  useRoot,
 639                             bool  corr,
 640                             bool  error)
 641 {
 642   //
 643   // Project onto the X axis
 644   //
 645   // Parameters:
 646   //    h         2D histogram
 647   //    name      New name
 648   //    firstbin  First bin to use
 649   //    lastbin   Last bin to use
 650   //    error     Whether to calculate errors
 651   //
 652   // Return:
 653   //    Newly created histogram or null
 654   //
 655   if (!h) return 0;
 656   if (useRoot)
 657     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 658
 659   TAxis* xaxis = h->GetXaxis();
 660   TAxis* yaxis = h->GetYaxis();
 661   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 662                           xaxis->GetXmin(), xaxis->GetXmax());
 663   static_cast<const TAttLine*>(h)->Copy(*ret);
 664   static_cast<const TAttFill*>(h)->Copy(*ret);
 665   static_cast<const TAttMarker*>(h)->Copy(*ret);
 666   ret->GetXaxis()->ImportAttributes(xaxis);
 667
 668   Int_t first = firstbin;
 669   Int_t last  = lastbin;
 670   if      (first < 0)                    first = 1;
 671   else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1;
 672   if      (last  < 0)                    last  = yaxis->GetNbins();
 673   else if (last  >= yaxis->GetNbins()+2) last  = yaxis->GetNbins()+1;
 674   if (last-first < 0) {
 675     AliWarningGeneral("AliBasedNdetaTask",
 676                       Form("Nothing to project [%d,%d]", first, last));
 677     return 0;
 678
 679   }
 680
 681   // Loop over X bins
 682   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 683   Int_t ybins = (last-first+1);
 684   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 685     Double_t content = 0;
 686     Double_t error2  = 0;
 687     Int_t    nbins   = 0;
 688
 689
 690     for (Int_t ybin = first; ybin <= last; ybin++) {
 691       Double_t c1 = h->GetCellContent(xbin, ybin);
 692       Double_t e1 = h->GetCellError(xbin, ybin);
 693
 694       // Ignore empty bins
 695       if (c1 < 1e-12) continue;
 696       if (e1 < 1e-12) {
 697         if (error) continue;
 698         e1 = 1;
 699       }
 700
 701       content    += c1;
 702       error2     += e1*e1;
 703       nbins++;
 704     } // for (ybin)
 705     if(content > 0 && nbins > 0) {
 706       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 707 #if 0
 708       AliWarningGeneral(ret->GetName(),
 709                         Form("factor @ %d is %d/%d -> %f",
 710                              xbin, ybins, nbins, factor));
 711 #endif
 712       if (error) {
 713         // Calculate weighted average
 714         ret->SetBinContent(xbin, content * factor);
 715         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 716       }
 717       else
 718         ret->SetBinContent(xbin, factor * content);
 719     }
 720   } // for (xbin)
 721
 722   return ret;
 723 }
 724
 725 //________________________________________________________________________
 726 void
 727 AliBasedNdetaTask::Terminate(Option_t *)
 728 {
 729   //
 730   // Called at end of event processing..
 731   //
 732   // This is called once in the master
 733   //
 734   // Parameters:
 735   //    option Not used
 736   //
 737   // Draw result to screen, or perform fitting, normalizations Called
 738   // once at the end of the query
 739   DGUARD(fDebug,1,"Process final merged results");
 740
 741   fSums = dynamic_cast<TList*> (GetOutputData(1));
 742   if(!fSums) {
 743     AliError("Could not retrieve TList fSums");
 744     return;
 745   }
 746
 747   fOutput = new TList;
 748   fOutput->SetName(Form("%s_result", GetName()));
 749   fOutput->SetOwner();
 750
 751   fSNNString     = fSums->FindObject("sNN");
 752   fSysString     = fSums->FindObject("sys");
 753   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
 754   fSchemeString  = fSums->FindObject("scheme");
 755   fTriggerString = fSums->FindObject("trigger");
 756
 757   if(fSysString && fSNNString &&
 758      fSysString->GetUniqueID() == AliForwardUtil::kPP)
 759     LoadNormalizationData(fSysString->GetUniqueID(),
 760                           fSNNString->GetUniqueID());
 761
 762   InitializeCentBins();
 763
 764   // Print before we loop
 765   Print();
 766
 767   // Loop over centrality bins
 768   TIter next(fListOfCentralities);
 769   CentralityBin* bin = 0;
 770   gStyle->SetPalette(1);
 771   THStack* dndetaStack        = new THStack("dndeta", "dN/d#eta");
 772   THStack* dndetaStackRebin   = new THStack(Form("dndeta_rebin%02d", fRebin),
 773                                             "dN_{ch}/d#eta");
 774   THStack* dndetaMCStack      = new THStack("dndetaMC", "dN_{ch}/d#eta");
 775   THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
 776                                             "dN_{ch}/d#eta");
 777
 778   TList* mclist = 0;
 779   TList* truthlist = 0;
 780
 781   if (fFinalMCCorrFile.Contains(".root")) {
 782     TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
 783     if(ftest) {
 784       mclist    = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
 785       truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
 786     }
 787     else
 788       AliWarning("MC analysis file invalid - no final MC correction possible");
 789   }
 790   Int_t style = GetMarker();
 791   Int_t color = GetColor();
 792
 793   AliInfo(Form("Marker style=%d, color=%d", style, color));
 794   while ((bin = static_cast<CentralityBin*>(next()))) {
 795
 796     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr,
 797              fTriggerEff, fTriggerEff0,
 798              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
 799              fTriggerMask, style, color, mclist, truthlist);
 800     if (fCentAxis && bin->IsAllBin()) continue;
 801     TH1* dndeta      = bin->GetResult(0, false, "");
 802     TH1* dndetaSym   = bin->GetResult(0, true,  "");
 803     TH1* dndetaMC    = bin->GetResult(0, false, "MC");
 804     TH1* dndetaMCSym = bin->GetResult(0, true,  "MC");
 805     if (dndeta)      dndetaStack->Add(dndeta);
 806     if (dndetaSym)   dndetaStack->Add(dndetaSym);
 807     if (dndetaMC)    dndetaMCStack->Add(dndetaMC);
 808     if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
 809     if (fRebin > 1) {
 810       dndeta      = bin->GetResult(fRebin, false, "");
 811       dndetaSym   = bin->GetResult(fRebin, true,  "");
 812       dndetaMC    = bin->GetResult(fRebin, false, "MC");
 813       dndetaMCSym = bin->GetResult(fRebin, true,  "MC");
 814       if (dndeta)      dndetaStackRebin->Add(dndeta);
 815       if (dndetaSym)   dndetaStackRebin->Add(dndetaSym);
 816       if (dndetaMC)    dndetaMCStackRebin->Add(dndetaMC);
 817       if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
 818     }
 819   }
 820   // Output the stack
 821   fOutput->Add(dndetaStack);
 822
 823   // If available output rebinned stack
 824   if (!dndetaStackRebin->GetHists() ||
 825       dndetaStackRebin->GetHists()->GetEntries() <= 0) {
 826     AliWarning("No rebinned histograms found");
 827     delete dndetaStackRebin;
 828     dndetaStackRebin = 0;
 829   }
 830   if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
 831
 832   // If available, output track-ref stack
 833   if (!dndetaMCStack->GetHists() ||
 834       dndetaMCStack->GetHists()->GetEntries() <= 0) {
 835     AliWarning("No MC histograms found");
 836     delete dndetaMCStack;
 837     dndetaMCStack = 0;
 838   }
 839   if (dndetaMCStack) fOutput->Add(dndetaMCStack);
 840
 841   // If available, output rebinned track-ref stack
 842   if (!dndetaMCStackRebin->GetHists() ||
 843       dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
 844     AliWarning("No rebinned MC histograms found");
 845     delete dndetaMCStackRebin;
 846     dndetaMCStackRebin = 0;
 847   }
 848   if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
 849
 850   // Output collision energy string
 851   if (fSNNString) {
 852     UShort_t sNN = fSNNString->GetUniqueID();
 853     TNamed* sNNObj = new TNamed(fSNNString->GetName(),
 854                                 AliForwardUtil::CenterOfMassEnergyString(sNN));
 855     sNNObj->SetUniqueID(sNN);
 856     fOutput->Add(sNNObj); // fSNNString->Clone());
 857   }
 858
 859   // Output collision system string
 860   if (fSysString) {
 861     UShort_t sys = fSysString->GetUniqueID();
 862     TNamed* sysObj = new TNamed(fSysString->GetName(),
 863                                 AliForwardUtil::CollisionSystemString(sys));
 864     sysObj->SetUniqueID(sys);
 865     fOutput->Add(sysObj); // fSysString->Clone());
 866   }
 867
 868   // Output centrality axis
 869   if (fCentAxis) fOutput->Add(fCentAxis);
 870
 871   // Output trigger string
 872   if (fTriggerString) {
 873     UShort_t mask = fTriggerString->GetUniqueID();
 874     TNamed* maskObj = new TNamed(fTriggerString->GetName(),
 875                                  AliAODForwardMult::GetTriggerString(mask));
 876     maskObj->SetUniqueID(mask);
 877     fOutput->Add(maskObj); // fTriggerString->Clone());
 878   }
 879
 880   // Normalization string
 881   if (fSchemeString) {
 882     UShort_t scheme = fSchemeString->GetUniqueID();
 883     TNamed* schemeObj = new TNamed(fSchemeString->GetName(),
 884                                    NormalizationSchemeString(scheme));
 885     schemeObj->SetUniqueID(scheme);
 886     fOutput->Add(schemeObj); // fSchemeString->Clone());
 887   }
 888
 889   // Output vertex axis
 890   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 891   vtxAxis->SetName("vtxAxis");
 892   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 893   fOutput->Add(vtxAxis);
 894
 895   // Output trigger efficiency and shape correction
 896   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff));
 897   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0));
 898   if (fShapeCorr) fOutput->Add(fShapeCorr);
 899
 900   TNamed* options = new TNamed("options","");
 901   TString str;
 902   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
 903   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
 904   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
 905   options->SetTitle(str);
 906   fOutput->Add(options);
 907
 908   PostData(2, fOutput);
 909 }
 910 //________________________________________________________________________
 911 void
 912 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 913 {
 914   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
 915   DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
 916          sys, energy);
 917   TString type("pp");
 918   TString snn("900");
 919   if(energy == 7000) snn.Form("7000");
 920   if(energy == 2750) snn.Form("2750");
 921
 922   if(fShapeCorr &&  (fTriggerEff != 1)) {
 923     AliInfo("Objects already set for normalization - no action taken");
 924     return;
 925   }
 926
 927   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
 928                                 "Normalization/normalizationHists_%s_%s.root",
 929                                 type.Data(),snn.Data()));
 930   if(!fin) {
 931     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
 932     return;
 933   }
 934
 935   // Shape correction
 936   if ((fNormalizationScheme & kShape) && !fShapeCorr) {
 937     TString trigName("All");
 938     if (fTriggerMask == AliAODForwardMult::kInel ||
 939         fTriggerMask == AliAODForwardMult::kNClusterGt0)
 940       trigName = "Inel";
 941     else if (fTriggerMask == AliAODForwardMult::kNSD)
 942       trigName = "NSD";
 943     else if (fTriggerMask == AliAODForwardMult::kInelGt0)
 944       trigName = "InelGt0";
 945     else {
 946       AliWarning(Form("Normalization for trigger %s not known, using all",
 947                       AliAODForwardMult::GetTriggerString(fTriggerMask)));
 948     }
 949
 950     TString shapeCorName(Form("h%sNormalization", trigName.Data()));
 951     TH2F*   shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
 952     if (shapeCor) SetShapeCorrection(shapeCor);
 953     else {
 954       AliWarning(Form("No shape correction found for %s", trigName.Data()));
 955     }
 956   }
 957
 958   // Trigger efficiency
 959   TString effName(Form("%sTriggerEff",
 960                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 961                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 962                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
 963                        "inelgt0" : "all"));
 964
 965   Double_t trigEff = 1;
 966   if (fNormalizationScheme & kTriggerEfficiency) {
 967     TObject* eff = fin->Get(effName);
 968     if (eff) AliForwardUtil::GetParameter(eff, trigEff);
 969   }
 970   if (fTriggerEff != 1) SetTriggerEff(trigEff);
 971   if (fTriggerEff < 0)  fTriggerEff = 1;
 972
 973   // Trigger efficiency
 974   TString eff0Name(Form("%sTriggerEff0",
 975                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 976                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 977                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
 978                        "inelgt0" : "all"));
 979
 980   Double_t trigEff0 = 1;
 981   if (fNormalizationScheme & kTriggerEfficiency) {
 982     TObject* eff = fin->Get(eff0Name);
 983     if (eff) AliForwardUtil::GetParameter(eff, trigEff0);
 984   }
 985   if (fTriggerEff0 != 1) SetTriggerEff0(trigEff0);
 986   if (fTriggerEff0 < 0)  fTriggerEff0 = 1;
 987
 988   // TEMPORARY FIX
 989   // Rescale the shape correction by the trigger efficiency
 990   if (fShapeCorr) {
 991     AliWarning(Form("Rescaling shape correction by trigger efficiency: "
 992                     "1/E_X=1/%f", trigEff));
 993     fShapeCorr->Scale(1. / trigEff);
 994   }
 995
 996   // Print - out
 997   if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
 998 }
 999
1000
1001 //________________________________________________________________________
1002 void
1003 AliBasedNdetaTask::Print(Option_t*) const
1004 {
1005   //
1006   // Print information
1007   //
1008   TString trigString("none");
1009   TString schemeString("none");
1010   TString sysString("unknown");
1011   TString sNNString("unknown");
1012   if (fTriggerString)
1013     trigString = AliAODForwardMult::GetTriggerString(fTriggerString->
1014                                                      GetUniqueID());
1015   if (fSchemeString)
1016     schemeString = NormalizationSchemeString(fSchemeString->GetUniqueID());
1017   if (fSysString)
1018     sysString = AliForwardUtil::CollisionSystemString(fSysString->
1019                                                       GetUniqueID());
1020   if (fSNNString)
1021     sNNString = AliForwardUtil::CenterOfMassEnergyString(fSNNString->
1022                                                          GetUniqueID());
1023
1024
1025   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
1026             << std::boolalpha
1027             << " Trigger:                    " << trigString << "\n"
1028             << " Vertex range:               [" << fVtxMin << ":"
1029             << fVtxMax << "]\n"
1030             << " Rebin factor:               " << fRebin << "\n"
1031             << " Cut edges:                  " << fCutEdges << "\n"
1032             << " Symmertrice:                " << fSymmetrice << "\n"
1033             << " Use TH2::ProjectionX:       " << fUseROOTProj << "\n"
1034             << " Correct for empty:          " << fCorrEmpty << "\n"
1035             << " Normalization scheme:       " << schemeString <<"\n"
1036             << " Trigger efficiency:         " << fTriggerEff << "\n"
1037             << " Bin-0 Trigger efficiency:   " << fTriggerEff0 << "\n"
1038             << " Shape correction:           " << (fShapeCorr ?
1039                                                    fShapeCorr->GetName() :
1040                                                    "none") << "\n"
1041             << " sqrt(s_NN):                 " << sNNString << "\n"
1042             << " Collision system:           " << sysString << "\n"
1043             << " Centrality bins:            " << (fCentAxis ? "" : "none");
1044   if (fCentAxis) {
1045     Int_t           nBins = fCentAxis->GetNbins();
1046     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
1047     for (Int_t i = 0; i <= nBins; i++)
1048       std::cout << (i==0 ? "" : "-") << bins[i];
1049   }
1050   std::cout << std::noboolalpha << std::endl;
1051
1052 }
1053
1054 //________________________________________________________________________
1055 TH1D*
1056 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
1057 {
1058   //
1059   // Make a copy of the input histogram and rebin that histogram
1060   //
1061   // Parameters:
1062   //    h  Histogram to rebin
1063   //
1064   // Return:
1065   //    New (rebinned) histogram
1066   //
1067   if (rebin <= 1) return 0;
1068
1069   Int_t nBins = h->GetNbinsX();
1070   if(nBins % rebin != 0) {
1071     AliWarningGeneral("AliBasedNdetaTask",
1072                       Form("Rebin factor %d is not a devisor of current number "
1073                            "of bins %d in the histogram %s",
1074                            rebin, nBins, h->GetName()));
1075     return 0;
1076   }
1077
1078   // Make a copy
1079   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
1080                                                h->GetName(), rebin)));
1081   tmp->Rebin(rebin);
1082   // Hist should be reset, as it otherwise messes with the cutEdges option
1083   tmp->Reset();
1084   tmp->SetDirectory(0);
1085
1086   // The new number of bins
1087   Int_t nBinsNew = nBins / rebin;
1088   for(Int_t i = 1;i<= nBinsNew; i++) {
1089     Double_t content = 0;
1090     Double_t sumw    = 0;
1091     Double_t wsum    = 0;
1092     Int_t    nbins   = 0;
1093     for(Int_t j = 1; j<=rebin;j++) {
1094       Int_t    bin = (i-1)*rebin + j;
1095       Double_t c   =  h->GetBinContent(bin);
1096       if (c <= 0)  {
1097         continue; // old TODO: check
1098         //content = -1; // new
1099         //break; // also new
1100       }
1101
1102       if (cutEdges) {
1103         if (h->GetBinContent(bin+1)<=0 ||
1104             h->GetBinContent(bin-1)<=0) {
1105 #if 0
1106           AliWarningGeneral("AliBasedNdetaTask",
1107                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
1108                                  bin, c, h->GetName(),
1109                                  bin+1, h->GetBinContent(bin+1),
1110                                  bin-1, h->GetBinContent(bin-1)));
1111 #endif
1112           continue;
1113         }
1114       }
1115       Double_t e =  h->GetBinError(bin);
1116       Double_t w =  1 / (e*e); // 1/c/c
1117       content    += c;
1118       sumw       += w;
1119       wsum       += w * c;
1120       nbins++;
1121     }
1122
1123     if(content > 0 && nbins > 0) {
1124       tmp->SetBinContent(i, wsum / sumw);
1125       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1126     }
1127   }
1128
1129   return tmp;
1130 }
1131
1132 //__________________________________________________________________
1133 TH1*
1134 AliBasedNdetaTask::Symmetrice(const TH1* h)
1135 {
1136   //
1137   // Make an extension of @a h to make it symmetric about 0
1138   //
1139   // Parameters:
1140   //    h Histogram to symmertrice
1141   //
1142   // Return:
1143   //    Symmetric extension of @a h
1144   //
1145   Int_t nBins = h->GetNbinsX();
1146   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1147   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1148   s->Reset();
1149   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1150   s->SetMarkerColor(h->GetMarkerColor());
1151   s->SetMarkerSize(h->GetMarkerSize());
1152   s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
1153   s->SetFillColor(h->GetFillColor());
1154   s->SetFillStyle(h->GetFillStyle());
1155   s->SetDirectory(0);
1156
1157   // Find the first and last bin with data
1158   Int_t first = nBins+1;
1159   Int_t last  = 0;
1160   for (Int_t i = 1; i <= nBins; i++) {
1161     if (h->GetBinContent(i) <= 0) continue;
1162     first = TMath::Min(first, i);
1163     last  = TMath::Max(last,  i);
1164   }
1165
1166   Double_t xfirst = h->GetBinCenter(first-1);
1167   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1168   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1169   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1170     s->SetBinContent(j, h->GetBinContent(i));
1171     s->SetBinError(j, h->GetBinError(i));
1172   }
1173   // Fill in overlap bin
1174   s->SetBinContent(l2+1, h->GetBinContent(first));
1175   s->SetBinError(l2+1, h->GetBinError(first));
1176   return s;
1177 }
1178
1179 //__________________________________________________________________
1180 Int_t
1181 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
1182 {
1183   Int_t  base   = bits & (0xFE);
1184   Bool_t hollow = bits & kHollow;
1185   switch (base) {
1186   case kCircle:       return (hollow ? 24 : 20);
1187   case kSquare:       return (hollow ? 25 : 21);
1188   case kUpTriangle:   return (hollow ? 26 : 22);
1189   case kDownTriangle: return (hollow ? 32 : 23);
1190   case kDiamond:      return (hollow ? 27 : 33);
1191   case kCross:        return (hollow ? 28 : 34);
1192   case kStar:         return (hollow ? 30 : 29);
1193   }
1194   return 1;
1195 }
1196 //__________________________________________________________________
1197 UShort_t
1198 AliBasedNdetaTask::GetMarkerBits(Int_t style)
1199 {
1200   UShort_t bits = 0;
1201   switch (style) {
1202   case 24: case 25: case 26: case 27: case 28: case 30: case 32:
1203     bits |= kHollow; break;
1204   }
1205   switch (style) {
1206   case 20: case 24: bits |= kCircle;       break;
1207   case 21: case 25: bits |= kSquare;       break;
1208   case 22: case 26: bits |= kUpTriangle;   break;
1209   case 23: case 32: bits |= kDownTriangle; break;
1210   case 27: case 33: bits |= kDiamond;      break;
1211   case 28: case 34: bits |= kCross;        break;
1212   case 29: case 30: bits |= kStar;         break;
1213   }
1214   return bits;
1215 }
1216 //__________________________________________________________________
1217 Int_t
1218 AliBasedNdetaTask::FlipHollowStyle(Int_t style)
1219 {
1220   UShort_t bits = GetMarkerBits(style);
1221   Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
1222   return ret;
1223 }
1224
1225 //====================================================================
1226 void
1227 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1228 {
1229   DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1230   TString n(GetHistName(0));
1231   TString n0(GetHistName(1));
1232   const char* postfix = GetTitle();
1233
1234   fSum = static_cast<TH2D*>(data->Clone(n));
1235   if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1236   fSum->SetDirectory(0);
1237   fSum->SetMarkerColor(col);
1238   fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1239   fSum->Reset();
1240   list->Add(fSum);
1241
1242   fSum0 = static_cast<TH2D*>(data->Clone(n0));
1243   if (postfix)
1244     fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1245   else
1246     fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1247   fSum0->SetDirectory(0);
1248   fSum0->SetMarkerColor(col);
1249   fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1250   fSum0->Reset();
1251   list->Add(fSum0);
1252
1253   fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1254   fEvents->SetDirectory(0);
1255   fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1256   fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1257   list->Add(fEvents);
1258 }
1259
1260 //____________________________________________________________________
1261 TString
1262 AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post)
1263 {
1264   TString n(name);
1265   if      (what == 1) n.Append("0");
1266   else if (what == 2) n.Append("Events");
1267   if (post && post[0] != '\0')  n.Append(post);
1268   return n;
1269 }
1270
1271 //____________________________________________________________________
1272 TString
1273 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1274 {
1275   return GetHistName(GetName(), what, GetTitle());
1276 }
1277
1278 //____________________________________________________________________
1279 void
1280 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1281 {
1282   DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1283   if (isZero) fSum0->Add(data);
1284   else        fSum->Add(data);
1285   fEvents->Fill(isZero ? 1 : 0);
1286 }
1287
1288 //____________________________________________________________________
1289 TH2D*
1290 AliBasedNdetaTask::Sum::CalcSum(TList*       output,
1291                                 Double_t&    ntotal,
1292                                 Double_t     epsilon0,
1293                                 Double_t     epsilon,
1294                                 Int_t        marker,
1295                                 Bool_t       rootProj,
1296                                 Bool_t       corrEmpty) const
1297 {
1298   DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1299
1300   // The return value `ret' is not scaled in anyway
1301   TH2D* ret      = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1302   ret->SetDirectory(0);
1303   ret->Reset();
1304   Int_t n        = Int_t(fEvents->GetBinContent(1));
1305   Int_t n0       = Int_t(fEvents->GetBinContent(2));
1306
1307   AliInfoF("Adding histograms %s(%d) and %s(%d) with weights %f and %f resp.",
1308            fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon, 1./epsilon0);
1309   DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.",
1310        fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1311   // Generate merged histogram
1312   ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon);
1313   ntotal = n / epsilon + n0 / epsilon0;
1314
1315   TList* out = new TList;
1316   out->SetOwner();
1317   const char* postfix = GetTitle();
1318   if (!postfix) postfix = "";
1319   out->SetName(Form("partial%s", postfix));
1320   output->Add(out);
1321
1322   // Now make copies, normalize them, and store in output list
1323   // Note, these are the only ones normalized here
1324   // These are mainly for diagnostics
1325   TH2D* sumCopy  = static_cast<TH2D*>(fSum->Clone("sum"));
1326   TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1327   TH2D* retCopy  = static_cast<TH2D*>(ret->Clone("sumAll"));
1328   sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1329   sumCopy->SetDirectory(0);
1330   sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1331   sum0Copy->SetDirectory(0);
1332   retCopy->SetMarkerStyle(marker);
1333   retCopy->SetDirectory(0);
1334
1335   Int_t nY      = fSum->GetNbinsY();
1336   Int_t o       = 0; // nY+1;
1337   TH1D* norm    = ProjectX(fSum,  "norm",    o, o, rootProj, corrEmpty, false);
1338   TH1D* norm0   = ProjectX(fSum0, "norm0",   o, o, rootProj, corrEmpty, false);
1339   TH1D* normAll = ProjectX(ret,   "normAll", o, o, rootProj, corrEmpty, false);
1340   norm->SetTitle("#eta coverage - >0-bin");
1341   norm0->SetTitle("#eta coverage - 0-bin");
1342   normAll->SetTitle("#eta coverage");
1343   norm->SetDirectory(0);
1344   norm0->SetDirectory(0);
1345   normAll->SetDirectory(0);
1346
1347   TH1D* sumCopyPx  = ProjectX(sumCopy,  "average",    1, nY,rootProj,corrEmpty);
1348   TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0",   1, nY,rootProj,corrEmpty);
1349   TH1D* retCopyPx  = ProjectX(retCopy,  "averageAll", 1, nY,rootProj,corrEmpty);
1350   sumCopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sumCopy->GetTitle()));
1351   sum0CopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sum0Copy->GetTitle()));
1352   retCopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", retCopy->GetTitle()));
1353   sumCopyPx->SetDirectory(0);
1354   sum0CopyPx->SetDirectory(0);
1355   retCopyPx->SetDirectory(0);
1356
1357   TH1D* phi    = ProjectX(fSum,  "phi",    nY+1, nY+1,rootProj,corrEmpty,false);
1358   TH1D* phi0   = ProjectX(fSum0, "phi0",   nY+1, nY+1,rootProj,corrEmpty,false);
1359   TH1D* phiAll = ProjectX(ret,   "phiAll", nY+1, nY+1,rootProj,corrEmpty,false);
1360   phi->SetTitle("#phi acceptance from dead strips - >0-bin");
1361   phi0->SetTitle("#phi acceptance from dead strips - 0-bin");
1362   phiAll->SetTitle("#phi acceptance from dead strips");
1363   phi->SetDirectory(0);
1364   phi0->SetDirectory(0);
1365   phiAll->SetDirectory(0);
1366
1367   const TH1D* cov    = (corrEmpty ? norm    : phi);
1368   const TH1D* cov0   = (corrEmpty ? norm0   : phi0);
1369   const TH1D* covAll = (corrEmpty ? normAll : phiAll);
1370
1371   // Here, we scale to the coverage (or phi acceptance)
1372   ScaleToCoverage(sumCopy,  cov);
1373   ScaleToCoverage(sum0Copy, cov0);
1374   ScaleToCoverage(retCopy,  covAll);
1375
1376   // Scale our 1D histograms
1377   sumCopyPx->Scale(1., "width");
1378   sum0CopyPx->Scale(1., "width");
1379   retCopyPx->Scale(1., "width");
1380
1381   AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1382                sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1383
1384   // Scale the normalization - they should be 1 at the maximum
1385   norm->Scale(n > 0   ? 1. / n  : 1);
1386   norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1387   normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1388
1389   // Scale the normalization - they should be 1 at the maximum
1390   phi->Scale(n > 0   ? 1. / n  : 1);
1391   phi0->Scale(n0 > 0 ? 1. / n0 : 1);
1392   phiAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1393
1394   out->Add(sumCopy);
1395   out->Add(sum0Copy);
1396   out->Add(retCopy);
1397   out->Add(sumCopyPx);
1398   out->Add(sum0CopyPx);
1399   out->Add(retCopyPx);
1400   out->Add(norm);
1401   out->Add(norm0);
1402   out->Add(normAll);
1403   out->Add(phi);
1404   out->Add(phi0);
1405   out->Add(phiAll);
1406
1407   AliInfoF("Returning  (1/%f * %s + 1/%f * %s), "
1408            "1/%f * %d + 1/%f * %d = %d",
1409            epsilon0, fSum0->GetName(), epsilon, fSum->GetName(),
1410            epsilon0, n0, epsilon, n, int(ntotal));
1411 #if 0
1412   for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1413     Double_t nc  = sum->GetBinContent(i, 0);
1414     Double_t nc0 = sum0->GetBinContent(i, 0);
1415     ret->SetBinContent(i, 0, nc + nc0); // Just count events
1416   }
1417 #endif
1418
1419   return ret;
1420 }
1421
1422 //====================================================================
1423 AliBasedNdetaTask::CentralityBin::CentralityBin()
1424   : TNamed("", ""),
1425     fSums(0),
1426     fOutput(0),
1427     fSum(0),
1428     fSumMC(0),
1429     fTriggers(0),
1430     fStatus(0),
1431     fLow(0),
1432     fHigh(0),
1433     fDoFinalMCCorrection(false),
1434     fSatelliteVertices(false),
1435     fDebug(0)
1436 {
1437   //
1438   // Constructor
1439   //
1440   DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin");
1441 }
1442 //____________________________________________________________________
1443 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
1444                                                 Short_t low, Short_t high)
1445   : TNamed(name, ""),
1446     fSums(0),
1447     fOutput(0),
1448     fSum(0),
1449     fSumMC(0),
1450     fTriggers(0),
1451     fStatus(0),
1452     fLow(low),
1453     fHigh(high),
1454     fDoFinalMCCorrection(false),
1455     fSatelliteVertices(false),
1456     fDebug(0)
1457 {
1458   //
1459   // Constructor
1460   //
1461   // Parameters:
1462   //    name Name used for histograms (e.g., Forward)
1463   //    low  Lower centrality cut in percent
1464   //    high Upper centrality cut in percent
1465   //
1466   DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: %s [%3d,%3d]",
1467          name,low,high);
1468   if (low <= 0 && high <= 0) {
1469     fLow  = 0;
1470     fHigh = 0;
1471     SetTitle("All centralities");
1472   }
1473   else {
1474     fLow  = low;
1475     fHigh = high;
1476     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1477   }
1478 }
1479 //____________________________________________________________________
1480 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1481   : TNamed(o),
1482     fSums(o.fSums),
1483     fOutput(o.fOutput),
1484     fSum(o.fSum),
1485     fSumMC(o.fSumMC),
1486     fTriggers(o.fTriggers),
1487     fStatus(o.fStatus),
1488     fLow(o.fLow),
1489     fHigh(o.fHigh),
1490     fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1491     fSatelliteVertices(o.fSatelliteVertices),
1492     fDebug(o.fDebug)
1493 {
1494   //
1495   // Copy constructor
1496   //
1497   // Parameters:
1498   //    other Object to copy from
1499   //
1500   DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin");
1501 }
1502 //____________________________________________________________________
1503 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1504 {
1505   //
1506   // Destructor
1507   //
1508   DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin");
1509
1510   // if (fSums) fSums->Delete();
1511   // if (fOutput) fOutput->Delete();
1512 }
1513
1514 //____________________________________________________________________
1515 AliBasedNdetaTask::CentralityBin&
1516 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1517 {
1518   //
1519   // Assignment operator
1520   //
1521   // Parameters:
1522   //    other Object to assign from
1523   //
1524   // Return:
1525   //    Reference to this
1526   //
1527   DGUARD(fDebug,3,"Centrality bin assignment");
1528   if (&o == this) return *this;
1529   SetName(o.GetName());
1530   SetTitle(o.GetTitle());
1531   fSums      = o.fSums;
1532   fOutput    = o.fOutput;
1533   fSum       = o.fSum;
1534   fSumMC     = o.fSumMC;
1535   fTriggers  = o.fTriggers;
1536   fStatus    = o.fStatus;
1537   fLow       = o.fLow;
1538   fHigh      = o.fHigh;
1539   fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1540   fSatelliteVertices = o.fSatelliteVertices;
1541
1542   return *this;
1543 }
1544 //____________________________________________________________________
1545 Int_t
1546 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const
1547 {
1548   if (IsAllBin()) return fallback;
1549   Float_t  fc       = (fLow+double(fHigh-fLow)/2) / 100;
1550   Int_t    nCol     = gStyle->GetNumberOfColors();
1551   Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
1552   Int_t    col      = gStyle->GetColorPalette(icol);
1553   return col;
1554 }
1555 //____________________________________________________________________
1556 const char*
1557 AliBasedNdetaTask::CentralityBin::GetListName() const
1558 {
1559   //
1560   // Get the list name
1561   //
1562   // Return:
1563   //    List Name
1564   //
1565   if (IsAllBin()) return "all";
1566   return Form("cent%03d_%03d", fLow, fHigh);
1567 }
1568 //____________________________________________________________________
1569 void
1570 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir, Int_t mask)
1571 {
1572   //
1573   // Create output objects
1574   //
1575   // Parameters:
1576   //    dir   Parent list
1577   //
1578   DGUARD(fDebug,1,"Create centrality bin output objects");
1579   fSums = new TList;
1580   fSums->SetName(GetListName());
1581   fSums->SetOwner();
1582   dir->Add(fSums);
1583
1584   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers", mask);
1585   fTriggers->SetDirectory(0);
1586
1587   fStatus = AliAODForwardMult::MakeStatusHistogram("status");
1588   fStatus->SetDirectory(0);
1589
1590   fSums->Add(fTriggers);
1591   fSums->Add(fStatus);
1592 }
1593 //____________________________________________________________________
1594 void
1595 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1596 {
1597   fDebug = lvl;
1598   if (fSum) fSum->fDebug = lvl;
1599   if (fSumMC) fSumMC->fDebug = lvl;
1600 }
1601
1602 //____________________________________________________________________
1603 Bool_t
1604 AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc)
1605 {
1606   const char* post = (mc ? "MC" : "");
1607   TString     sn   = Sum::GetHistName(GetName(),0,post);
1608   TString     sn0  = Sum::GetHistName(GetName(),1,post);
1609   TString     ev   = Sum::GetHistName(GetName(),2,post);
1610   TH2D* sum      = static_cast<TH2D*>(list->FindObject(sn));
1611   TH2D* sum0     = static_cast<TH2D*>(list->FindObject(sn0));
1612   TH1I* events   = static_cast<TH1I*>(list->FindObject(ev));
1613   if (!sum || !sum0 || !events) {
1614     AliWarningF("Failed to find one or more histograms: "
1615                 "%s (%p) %s (%p) %s (%p)",
1616                 sn.Data(), sum,
1617                 sn0.Data(), sum0,
1618                 ev.Data(), events);
1619     return false;
1620   }
1621   Sum* ret     = new Sum(GetName(), post);
1622   ret->fSum    = sum;
1623   ret->fSum0   = sum0;
1624   ret->fEvents = events;
1625   ret->fDebug  = fDebug;
1626   if (mc) fSumMC = ret;
1627   else    fSum   = ret;
1628
1629   return true;
1630 }
1631
1632 //____________________________________________________________________
1633 void
1634 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1635 {
1636   //
1637   // Create sum histogram
1638   //
1639   // Parameters:
1640   //    data  Data histogram to clone
1641   //    mc    (optional) MC histogram to clone
1642   //
1643   DGUARD(fDebug,1,"Create centrality bin sums from %s",
1644          data ? data->GetName() : "(null)");
1645   if (data) {
1646     fSum = new Sum(GetName(),"");
1647     fSum->Init(fSums, data, GetColor());
1648     fSum->fDebug = fDebug;
1649   }
1650
1651   // If no MC data is given, then do not create MC sum histogram
1652   if (!mc) return;
1653
1654   fSumMC = new Sum(GetName(), "MC");
1655   fSumMC->Init(fSums, mc, GetColor());
1656   fSumMC->fDebug = fDebug;
1657 }
1658
1659 //____________________________________________________________________
1660 Bool_t
1661 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1662                                              Int_t triggerMask,
1663                                              Double_t vzMin, Double_t vzMax)
1664 {
1665   //
1666   // Check the trigger, vertex, and centrality
1667   //
1668   // Parameters:
1669   //    aod Event input
1670   //
1671   // Return:
1672   //    true if the event is to be used
1673   //
1674   if (!forward) return false;
1675
1676   DGUARD(fDebug,2,"Check the event");
1677   // We do not check for centrality here - it's already done
1678   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0,
1679                              fTriggers, fStatus);
1680 }
1681
1682
1683 //____________________________________________________________________
1684 Bool_t
1685 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1686                                                Int_t triggerMask, Bool_t isZero,
1687                                                Double_t vzMin, Double_t vzMax,
1688                                                const TH2D* data, const TH2D* mc)
1689 {
1690   //
1691   // Process an event
1692   //
1693   // Parameters:
1694   //    forward     Forward data (for trigger, vertex, & centrality)
1695   //    triggerMask Trigger mask
1696   //    vzMin       Minimum IP z coordinate
1697   //    vzMax       Maximum IP z coordinate
1698   //    data        Data histogram
1699   //    mc          MC histogram
1700   //
1701   DGUARD(fDebug,1,"Process one event for %s a given centrality bin",
1702          data ? data->GetName() : "(null)");
1703   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return false;
1704   if (!data) return false;
1705   if (!fSum) CreateSums(data, mc);
1706
1707   fSum->Add(data, isZero);
1708   if (mc) fSumMC->Add(mc, isZero);
1709
1710   return true;
1711 }
1712
1713 //________________________________________________________________________
1714 Double_t
1715 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1716                                                 UShort_t    scheme,
1717                                                 Double_t    trigEff,
1718                                                 Double_t&   ntotal,
1719                                                 TString*    text) const
1720 {
1721   //
1722   // Calculate normalization
1723   //
1724   // Parameters:
1725   //    t       Trigger histogram
1726   //    scheme  Normaliztion scheme
1727   //    trigEff From MC
1728   //    ntotal  On return, contains the number of events.
1729   //
1730   DGUARD(fDebug,1,"Normalize centrality bin %s [%3d-%3d%%] with %s",
1731          GetName(), fLow, fHigh, t.GetName());
1732   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1733   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1734   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1735   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1736   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1737   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1738   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1739   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1740   Double_t nAccepted   = ntotal;
1741   ntotal               = 0;
1742
1743   if (nTriggered <= 0.1) {
1744     AliError("Number of triggered events <= 0");
1745     return -1;
1746   }
1747   if (nWithVertex <= 0.1) {
1748     AliError("Number of events with vertex <= 0");
1749     return -1;
1750   }
1751   ntotal          = nAccepted;
1752   Double_t vtxEff = nWithVertex / nTriggered;
1753   Double_t scaler = 1;
1754   Double_t beta   = nA + nC - 2*nE;
1755
1756
1757   TString rhs("N = N_acc");
1758   if (!(scheme & kZeroBin)) {
1759     if (scheme & kEventLevel) {
1760       ntotal = nAccepted / vtxEff;
1761       scaler = vtxEff;
1762       AliInfoF("Calculating event normalisation as\n"
1763                " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1764                Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1765                ntotal, scaler);
1766       if (scheme & kBackground) {
1767         //          1            E_V             E_V
1768         //   s = --------- = ------------- = ------------
1769         //        1 - beta   1 - beta E_V    1 - beta N_V
1770         //       ---  ----       --------        ---- ---
1771         //       E_V  N_V          N_V           N_V  N_T
1772         //
1773         //          E_V
1774         //     = ------------
1775         //        1 - beta
1776         //            ----
1777         //             N_T
1778         //
1779         ntotal -= nAccepted * beta / nWithVertex;
1780         // This one is direct and correct.
1781         // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1782         // A simpler expresion
1783         scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1784         AliInfo(Form("Calculating event normalisation as\n"
1785                      " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1786                      " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1787                      Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1788                      nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1789                      Int_t(nWithVertex), ntotal, scaler));
1790         rhs.Append("(1/eps_V - beta/N_vtx)");
1791       } // Background
1792       else
1793         rhs.Append("/eps_V");
1794     } // Event-level
1795     if (scheme & kTriggerEfficiency) {
1796       ntotal /= trigEff;
1797       scaler *= trigEff;
1798       AliInfo(Form("Correcting for trigger efficiency:\n"
1799                    " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1800                    trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1801       rhs.Append("/eps_T");
1802     } // Trigger efficiency
1803   }
1804   else  {
1805     // Calculate as
1806     //
1807     //  N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1808     //    = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1809     //    = N_A (1 + 1/E_X (1/E_V - 1 - beta / N_V))
1810     //
1811     //  s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1812     //    = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1813     //
1814     if (!(scheme & kBackground)) beta = 0;
1815     ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1816                                          - beta / nWithVertex));
1817     scaler = nWithVertex / (nWithVertex +
1818                             1/trigEff * (nTriggered-nWithVertex-beta));
1819     AliInfo(Form("Calculating event normalisation as\n"
1820                  "  beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1821                  "  N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1822                  "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1823                  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1824                  Int_t(nAccepted), trigEff, Int_t(nTriggered),
1825                  Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1826                  ntotal, scaler));
1827     rhs.Append("(1+1/eps_T(1/eps_V-1-beta/N_vtx))");
1828   }
1829
1830   if (text) {
1831     text->Append(Form("%-40s = %d\n", "N_all",            UInt_t(nAll)));
1832     text->Append(Form("%-40s = %d\n", "N_acc",            UInt_t(nAccepted)));
1833     text->Append(Form("%-40s = %d\n", "N_trg",            UInt_t(nTriggered)));
1834     text->Append(Form("%-40s = %d\n", "N_vtx",            UInt_t(nWithVertex)));
1835     text->Append(Form("%-40s = %d\n", "N_B",              UInt_t(nB)));
1836     text->Append(Form("%-40s = %d\n", "N_A",              UInt_t(nA)));
1837     text->Append(Form("%-40s = %d\n", "N_C",              UInt_t(nC)));
1838     text->Append(Form("%-40s = %d\n", "N_E",              UInt_t(nE)));
1839     text->Append(Form("%-40s = %d\n", "beta = N_A + N_C - 2N_E",UInt_t(beta)));
1840     text->Append(Form("%-40s = %f\n", "eps_V = N_vtx/N_trg",vtxEff));
1841     text->Append(Form("%-40s = %f\n", "eps_T",            trigEff));
1842     text->Append(Form("%-40s = %f\n", rhs.Data(),         ntotal));
1843   }
1844   AliInfo(Form("\n"
1845                " Total of        %9d events for %s\n"
1846                "  of these       %9d have an offline trigger\n"
1847                "  of these N_T = %9d has the selected trigger\n"
1848                "  of these N_V = %9d has a vertex\n"
1849                "  of these N_A = %9d were in the selected range\n"
1850                "  Triggers by hardware type:\n"
1851                "    N_b   =  %9d\n"
1852                "    N_ac  =  %9d (%d+%d)\n"
1853                "    N_e   =  %9d\n"
1854                "  Vertex efficiency:          %f\n"
1855                "  Trigger efficiency:         %f\n"
1856                "  Total number of events: N = %f\n"
1857                "  Scaler (N_A/N):             %f\n"
1858                "  %25s = %f",
1859                Int_t(nAll), GetTitle(), Int_t(nOffline),
1860                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1861                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1862                vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal));
1863   return scaler;
1864 }
1865
1866 //________________________________________________________________________
1867 const char*
1868 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1869                                                 Bool_t sym,
1870                                                 const char* postfix) const
1871 {
1872   static TString n;
1873   n = Form("dndeta%s%s",GetName(), postfix);
1874   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1875   if (sym)       n.Append("_mirror");
1876   return n.Data();
1877 }
1878 //________________________________________________________________________
1879 TH1*
1880 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1881                                             Bool_t sym,
1882                                             const char* postfix) const
1883 {
1884   if (!fOutput) {
1885     AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
1886                     fLow, fHigh));
1887     return 0;
1888   }
1889   TString  n = GetResultName(rebin, sym, postfix);
1890   TObject* o = fOutput->FindObject(n.Data());
1891   if (!o) {
1892     // AliWarning(Form("Object %s not found in output list", n.Data()));
1893     return 0;
1894   }
1895   return static_cast<TH1*>(o);
1896 }
1897
1898 //________________________________________________________________________
1899 void
1900 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1901                                              const char* postfix,
1902                                              bool        rootProj,
1903                                              bool        corrEmpty,
1904                                              const TH2F* shapeCorr,
1905                                              Double_t    scaler,
1906                                              bool        symmetrice,
1907                                              Int_t       rebin,
1908                                              bool        cutEdges,
1909                                              Int_t       marker,
1910                                              Int_t       color,
1911                                              TList*      mclist,
1912                                              TList*      truthlist)
1913 {
1914   //
1915   // Generate the dN/deta result from input
1916   //
1917   // Parameters:
1918   //     sum        Sum of 2D hists
1919   //     postfix    Post fix on names
1920   //     rootProj   Whether to use ROOT TH2::ProjectionX
1921   //     corrEmpty  Correct for empty bins
1922   //     shapeCorr  Shape correction to use
1923   //     scaler     Event-level normalization scaler
1924   //     symmetrice Whether to make symmetric extensions
1925   //     rebin      Whether to rebin
1926   //     cutEdges   Whether to cut edges when rebinning
1927   //
1928   DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1929   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1930                                                      GetName(), postfix)));
1931
1932   TH1D* accNorm = 0;
1933   Int_t nY      = sum->GetNbinsY();
1934   // Hack HHD Hans test code to manually remove FMD2 dead channel (but
1935   // it is on outer?)
1936   //
1937   // cholm comment: The original hack has been moved to
1938   // AliForwarddNdetaTask::CheckEventData - this simplifies things a
1939   // great deal, and we could in principle use the new phi acceptance.
1940   //
1941   // However, since we may have filtered out the dead sectors in the
1942   // AOD production already, we can't be sure we can recalculate the
1943   // phi acceptance correctly, so for now, we rely on fCorrEmpty being set.
1944   Int_t o       = (corrEmpty ? 0 : nY+1);
1945   accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), o, o,
1946                      rootProj, corrEmpty, false);
1947   accNorm->SetDirectory(0);
1948
1949   // ---- Scale by shape correction ----------------------------------
1950   if (shapeCorr) copy->Divide(shapeCorr);
1951   else AliInfo("No shape correction specified, or disabled");
1952
1953   // --- Normalize to the coverage -----------------------------------
1954   if (corrEmpty) {
1955     ScaleToCoverage(copy, accNorm);
1956     // --- Event-level normalization ---------------------------------
1957     copy->Scale(scaler);
1958   }
1959
1960   // --- Project on X axis -------------------------------------------
1961   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1962                           1, nY, rootProj, corrEmpty);
1963   dndeta->SetDirectory(0);
1964   // Event-level normalization
1965   if (!corrEmpty) {
1966     ScaleToCoverage(dndeta, accNorm);
1967     dndeta->Scale(scaler);
1968   }
1969   dndeta->Scale(1., "width");
1970   copy->Scale(1., "width");
1971
1972   TH1D* dndetaMCCorrection = 0;
1973   TList* centlist          = 0;
1974   TH1D* dndetaMCtruth      = 0;
1975   TList* truthcentlist     = 0;
1976
1977   // --- Possible final correction to <MC analysis> / <MC truth> -----
1978   // we get the rebinned distribution for satellite to make the correction
1979   TString rebinSuf(fSatelliteVertices ? "_rebin05" : "");
1980   if(mclist)
1981     centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1982   if(centlist)
1983     dndetaMCCorrection =
1984       static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s%s",
1985                                                    GetName(), postfix,
1986                                                    rebinSuf.Data())));
1987   if (truthlist)
1988     truthcentlist =
1989       static_cast<TList*>(truthlist->FindObject(GetListName()));
1990   if (truthcentlist)
1991     // TODO here new is "dndetaTruth"
1992     dndetaMCtruth =
1993       static_cast<TH1D*>(truthcentlist->FindObject(Form("dndetaMCTruth%s",
1994                                                         rebinSuf.Data())));
1995
1996   if (dndetaMCCorrection && dndetaMCtruth) {
1997     AliInfo("Correcting with final MC correction");
1998     TString testString(dndetaMCCorrection->GetName());
1999
2000     // We take the measured MC dN/deta and divide with truth
2001     dndetaMCCorrection->Divide(dndetaMCtruth);
2002     dndetaMCCorrection->SetTitle("Final MC correction");
2003     dndetaMCCorrection->SetName("finalMCCorr");
2004     for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
2005       if(dndetaMCCorrection->GetBinContent(m) < 0.5 ||
2006          dndetaMCCorrection->GetBinContent(m) > 1.75) {
2007         dndetaMCCorrection->SetBinContent(m,1.);
2008         dndetaMCCorrection->SetBinError(m,0.1);
2009       }
2010     }
2011     // Applying the correction
2012     if (!fSatelliteVertices)
2013       // For non-satellites we took the same binning, so we do a straight
2014       // division
2015       dndeta->Divide(dndetaMCCorrection);
2016     else {
2017       // For satelitte events, we took coarser binned histograms, so
2018       // we need to do a bit more
2019       for(Int_t m = 1; m <= dndeta->GetNbinsX(); m++) {
2020         if(dndeta->GetBinContent(m) <= 0.01 ) continue;
2021
2022         Double_t eta     = dndeta->GetXaxis()->GetBinCenter(m);
2023         Int_t    bin     = dndetaMCCorrection->GetXaxis()->FindBin(eta);
2024         Double_t mccorr  = dndetaMCCorrection->GetBinContent(bin);
2025         Double_t mcerror = dndetaMCCorrection->GetBinError(bin);
2026         if (mccorr < 1e-6) {
2027           dndeta->SetBinContent(m, 0);
2028           dndeta->SetBinError(m, 0);
2029         }
2030         Double_t value   = dndeta->GetBinContent(m);
2031         Double_t error   = dndeta->GetBinError(m);
2032         Double_t sumw2   = (error   * error   * mccorr * mccorr +
2033                             mcerror * mcerror * value  * value);
2034         dndeta->SetBinContent(m,value/mccorr) ;
2035         dndeta->SetBinError(m,TMath::Sqrt(sumw2)/mccorr/mccorr);
2036       }
2037     }
2038   }
2039   else
2040     AliInfo("No final MC correction applied");
2041
2042   // --- Set some histogram attributes -------------------------------
2043   TString post;
2044   Int_t rColor = GetColor(color);
2045   if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
2046   SetHistogramAttributes(dndeta,  rColor, marker,
2047                          Form("ALICE %s%s", GetName(), post.Data()));
2048   SetHistogramAttributes(accNorm, rColor, marker,
2049                          Form("ALICE %s normalisation%s",
2050                               GetName(), post.Data()));
2051
2052   // --- Make symmetric extensions and rebinnings --------------------
2053   if (symmetrice) fOutput->Add(Symmetrice(dndeta));
2054   fOutput->Add(dndeta);
2055   fOutput->Add(accNorm);
2056   fOutput->Add(copy);
2057   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
2058   if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
2059   if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
2060
2061   // HHD Test of dN/deta in phi bins add flag later?
2062   //
2063   // cholm comment: We disable this for now
2064 #if 0
2065   for (Int_t nn=1; nn <= sum->GetNbinsY(); nn++) {
2066     TH1D* dndeta_phi = ProjectX(copy, Form("dndeta%s%s_phibin%d",
2067                                            GetName(), postfix, nn),
2068                                 nn, nn, rootProj, corrEmpty);
2069     dndeta_phi->SetDirectory(0);
2070     // Event-level normalization
2071     dndeta_phi->Scale(TMath::Pi()/10., "width");
2072
2073     if(centlist)
2074       dndetaMCCorrection =
2075         static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s_phibin%d",
2076                                                      GetName(), postfix,nn)));
2077     if(truthcentlist)
2078       dndetaMCtruth
2079         = static_cast<TH1D*>(truthcentlist->FindObject("dndetaMCTruth"));
2080
2081     if (dndetaMCCorrection && dndetaMCtruth) {
2082       AliInfo("Correcting with final MC correction");
2083       TString testString(dndetaMCCorrection->GetName());
2084       dndetaMCCorrection->Divide(dndetaMCtruth);
2085       dndetaMCCorrection->SetTitle(Form("Final_MC_correction_phibin%d",nn));
2086       dndetaMCCorrection->SetName(Form("Final_MC_correction_phibin%d",nn));
2087       for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
2088         if(dndetaMCCorrection->GetBinContent(m) < 0.25 ||
2089            dndetaMCCorrection->GetBinContent(m) > 1.75) {
2090           dndetaMCCorrection->SetBinContent(m,1.);
2091           dndetaMCCorrection->SetBinError(m,0.1);
2092         }
2093       }
2094       //Applying the correction
2095       dndeta_phi->Divide(dndetaMCCorrection);
2096     }
2097     fOutput->Add(dndeta_phi);
2098     fOutput->Add(Rebin(dndeta_phi, rebin, cutEdges));
2099     if(dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
2100   } // End of phi
2101 #endif
2102 }
2103
2104 //________________________________________________________________________
2105 void
2106 AliBasedNdetaTask::CentralityBin::End(TList*      sums,
2107                                       TList*      results,
2108                                       UShort_t    scheme,
2109                                       const TH2F* shapeCorr,
2110                                       Double_t    trigEff,
2111                                       Double_t    trigEff0,
2112                                       Bool_t      symmetrice,
2113                                       Int_t       rebin,
2114                                       Bool_t      rootProj,
2115                                       Bool_t      corrEmpty,
2116                                       Bool_t      cutEdges,
2117                                       Int_t       triggerMask,
2118                                       Int_t       marker,
2119                                       Int_t       color,
2120                                       TList*      mclist,
2121                                       TList*      truthlist)
2122 {
2123   //
2124   // End of processing
2125   //
2126   // Parameters:
2127   //    sums        List of sums
2128   //    results     Output list of results
2129   //    shapeCorr   Shape correction or nil
2130   //    trigEff     Trigger efficiency
2131   //    symmetrice  Whether to symmetrice the results
2132   //    rebin       Whether to rebin the results
2133   //    corrEmpty   Whether to correct for empty bins
2134   //    cutEdges    Whether to cut edges when rebinning
2135   //    triggerMask Trigger mask
2136   //
2137   DGUARD(fDebug,1,"End centrality bin procesing");
2138
2139   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
2140   if(!fSums) {
2141     AliError("Could not retrieve TList fSums");
2142     return;
2143   }
2144
2145   fOutput = new TList;
2146   fOutput->SetName(GetListName());
2147   fOutput->SetOwner();
2148   results->Add(fOutput);
2149
2150   if (!fSum) {
2151     if (!ReadSum(fSums, false)) {
2152       AliInfo("This task did not produce any output");
2153       return;
2154     }
2155   }
2156   if (!fSumMC) ReadSum(fSums, true);
2157
2158   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
2159
2160   if (!fTriggers) {
2161     AliError("Couldn't find histogram 'triggers' in list");
2162     return;
2163   }
2164
2165   // --- Get normalization scaler ------------------------------------
2166   Double_t epsilonT  = trigEff;
2167   Double_t epsilonT0 = trigEff0;
2168   AliInfoF("Using epsilonT=%f, epsilonT0=%f for %d",
2169            epsilonT, epsilonT0, triggerMask);
2170 #if 0
2171   // These hard-coded trigger efficiencies are not used anymore, and
2172   // are only left in the code for reference.  We should remove this
2173   // soon.
2174   if (triggerMask == AliAODForwardMult::kNSD) {
2175     // This is a local change
2176     epsilonT = 0.96; // New value since HHD code 29/08/2012, why?
2177     //epsilonT = 0.92; //First paper...
2178     //epsilonT = 0.954; //First paper...
2179     //epsilonT = 1.03; //phojet
2180     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
2181   }
2182   else if (triggerMask == AliAODForwardMult::kInel) {
2183     // This is a local change
2184     epsilonT = 0.934; // New value since HHD code 29/08/2012, why?
2185     // 900 GeV Inel eff from Martin
2186     //epsilonT = 0.916;
2187     //epsilonT = 0.97;  //phojet
2188
2189     AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
2190   }
2191   if (scheme & kZeroBin) {
2192     if (triggerMask==AliAODForwardMult::kInel)
2193       epsilonT0 = 0.785021; // 0.100240;
2194     else if (triggerMask==AliAODForwardMult::kInelGt0)
2195       epsilonT0 = 0;
2196     else if (triggerMask==AliAODForwardMult::kNSD)
2197       epsilonT0 = .706587;
2198     epsilonT = 1;
2199     AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
2200                     epsilonT0));
2201   }
2202 #endif
2203
2204   // Get our histograms
2205   Double_t nSum   = 0;
2206   TH2D*    sum    = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT,
2207                                  marker, rootProj, corrEmpty);
2208   Double_t nSumMC = 0;
2209   TH2D*    sumMC  = 0;
2210   if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC,
2211                                       epsilonT0, epsilonT, marker,
2212                                       rootProj, corrEmpty);
2213   if (!sum) {
2214     AliError("Failed to get sum from summer - bailing out");
2215     return;
2216   }
2217
2218   TString  text;
2219   Double_t ntotal = nSum;
2220   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
2221   if (scaler < 0) {
2222     AliError("Failed to calculate normalization - bailing out");
2223     return;
2224   }
2225   fOutput->Add(fTriggers->Clone());
2226   fOutput->Add(new TNamed("normCalc", text.Data()));
2227
2228   // --- Make result and store ---------------------------------------
2229   MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
2230              scaler, symmetrice, rebin, cutEdges, marker, color,
2231              mclist, truthlist);
2232
2233   // --- Process result from TrackRefs -------------------------------
2234   if (sumMC)
2235     MakeResult(sumMC, "MC", rootProj, corrEmpty,
2236                (scheme & kShape) ? shapeCorr : 0,
2237                scaler, symmetrice, rebin, cutEdges,
2238                GetMarkerStyle(GetMarkerBits(marker)+4), color,
2239                mclist, truthlist);
2240
2241   // Temporary stuff
2242   // if (!IsAllBin()) return;
2243
2244 }
2245 //____________________________________________________________________
2246 Bool_t
2247 AliBasedNdetaTask::ApplyEmpiricalCorrection(const AliAODForwardMult* aod,
2248                                             TH2D* data)
2249 {
2250   if (!fglobalempiricalcorrection || !data)
2251     return true;
2252   Float_t zvertex=aod->GetIpZ();
2253   Int_t binzvertex=fglobalempiricalcorrection->GetXaxis()->FindBin(zvertex);
2254   if(binzvertex<1||binzvertex>fglobalempiricalcorrection->GetNbinsX())
2255     return false;
2256   for (int i=1;i<=data->GetNbinsX();i++) {
2257     Int_t bincorrection=fglobalempiricalcorrection->GetYaxis()
2258       ->FindBin(data->GetXaxis()->GetBinCenter(i));
2259     if(bincorrection<1||bincorrection>fglobalempiricalcorrection->GetNbinsY())
2260       return false;
2261     Float_t correction=fglobalempiricalcorrection
2262       ->GetBinContent(binzvertex,bincorrection);
2263     if(correction<0.001) {
2264       data->SetBinContent(i,0,0);
2265       data->SetBinContent(i,data->GetNbinsY()+1,0);
2266     }
2267     for(int j=1;j<=data->GetNbinsY();j++) {
2268       if (data->GetBinContent(i,j)>0.0) {
2269         data->SetBinContent(i,j,data->GetBinContent(i,j)*correction);
2270         data->SetBinError(i,j,data->GetBinError(i,j)*correction);
2271       }
2272     }
2273   }
2274   return true;
2275 }
2276
2277 //
2278 // EOF
2279 //