PWGLF/FORWARD/analysis2/AliForwardMultDists.cxx

   1 #include "AliForwardMultDists.h"
   2 #include <AliForwardUtil.h>
   3 #include <AliAODForwardMult.h>
   4 #include <AliAODCentralMult.h>
   5 #include <AliAODEvent.h>
   6 #include <AliLog.h>
   7 #include <TH1.h>
   8 #include <TH2.h>
   9 #include <TVector2.h>
  10 #include <THStack.h>
  11 #include <TParameter.h>
  12 #include <TError.h>
  13 #include <iostream>
  14 #include <iomanip>
  15
  16 //====================================================================
  17 namespace {
  18   /**
  19    * Marker styles
  20    */
  21   enum {
  22     kSolid        = 0x000,
  23     kHollow       = 0x001,
  24     kCircle       = 0x002,
  25     kSquare       = 0x004,
  26     kUpTriangle   = 0x006,
  27     kDownTriangle = 0x008,
  28     kDiamond      = 0x00a,
  29     kCross        = 0x00c,
  30     kStar         = 0x00e
  31   };
  32   /**
  33    * Get a ROOT marker style from bit options
  34    *
  35    * @param bits Bit mask of options
  36    *
  37    * @return ROOT marker style
  38    */
  39   Int_t GetMarkerStyle(UShort_t bits)
  40   {
  41     Int_t  base   = bits & (0xFE);
  42     Bool_t hollow = bits & kHollow;
  43     switch (base) {
  44     case kCircle:       return (hollow ? 24 : 20);
  45     case kSquare:       return (hollow ? 25 : 21);
  46     case kUpTriangle:   return (hollow ? 26 : 22);
  47     case kDownTriangle: return (hollow ? 32 : 23);
  48     case kDiamond:      return (hollow ? 27 : 33);
  49     case kCross:        return (hollow ? 28 : 34);
  50     case kStar:         return (hollow ? 30 : 29);
  51     }
  52     return 1;
  53   }
  54 #if 0
  55   /**
  56    * Get the marker option bits from a ROOT style
  57    *
  58    * @param style ROOT style
  59    *
  60    * @return Pattern of marker options
  61    */
  62   UShort_t GetMarkerBits(Int_t style)
  63   {
  64     UShort_t bits = 0;
  65     switch (style) {
  66     case 24: case 25: case 26: case 27: case 28: case 30: case 32:
  67       bits |= kHollow; break;
  68     }
  69     switch (style) {
  70     case 20: case 24: bits |= kCircle;       break;
  71     case 21: case 25: bits |= kSquare;       break;
  72     case 22: case 26: bits |= kUpTriangle;   break;
  73     case 23: case 32: bits |= kDownTriangle; break;
  74     case 27: case 33: bits |= kDiamond;      break;
  75     case 28: case 34: bits |= kCross;        break;
  76     case 29: case 30: bits |= kStar;         break;
  77     }
  78     return bits;
  79   }
  80 #endif
  81   static Int_t GetIndexMarker(Int_t idx)
  82   {
  83     const UShort_t nMarkers = 7;
  84     UShort_t markers[] = {
  85       kCircle,
  86       kSquare,
  87       kUpTriangle,
  88       kDownTriangle,
  89       kDiamond,
  90       kCross,
  91       kStar
  92     };
  93
  94     return markers[idx % nMarkers];
  95   }
  96 #if 0
  97   /**
  98    * Flip the 'hollow-ness' of a marker
  99    *
 100    * @param style ROOT Style
 101    *
 102    * @return ROOT style
 103    */
 104   Int_t FlipHollowStyle(Int_t style)
 105   {
 106     UShort_t bits = GetMarkerBits(style);
 107     Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
 108     return ret;
 109   }
 110 #endif
 111 }
 112
 113 //====================================================================
 114 AliForwardMultDists::BinSpec::BinSpec(Double_t etaMin,
 115                                       Double_t etaMax,
 116                                       Double_t nchLow)
 117   : fEtaMin(etaMin), fEtaMax(etaMax), fLow(nchLow), fN(0), fD(0), fAxis(1,0,1)
 118 {
 119 }
 120 //____________________________________________________________________
 121 void
 122 AliForwardMultDists::BinSpec::Push(UShort_t n, Double_t d)
 123 {
 124   Int_t l = fN.GetSize();
 125   fN.Set(l+1); // Not terribly efficient, but that's ok because we do
 126   fD.Set(l+1); // this infrequently
 127   fN[l] = n;
 128   fD[l] = d;
 129 }
 130 //____________________________________________________________________
 131 const TAxis&
 132 AliForwardMultDists::BinSpec::Axis() const
 133 {
 134   if (fAxis.GetNbins() > 1) return fAxis;
 135   if (fN.GetSize() <= 0) return fAxis;
 136   if (fN.GetSize() == 1) {
 137     // Exactly one spec,
 138     fAxis.Set(fN[0], fLow, fN[0]*fD[0]+fLow);
 139   }
 140   else {
 141     Int_t n = Int_t(fN.GetSum());
 142     Int_t l = fN.GetSize();
 143     Int_t i = 0;
 144     TArrayD b(n+1);
 145     b[0]    = fLow;
 146     for (Int_t j = 0; j < l; j++) {
 147       for (Int_t k = 0; k < fN[j]; k++) {
 148         b[i+1] = b[i] + fD[j];
 149         i++;
 150       }
 151     }
 152     if (i != n) {
 153       ::Warning("Axis", "Assigned bins, and number of bins not equal");
 154       n = TMath::Min(i, n);
 155     }
 156     fAxis.Set(n, b.GetArray());
 157   }
 158   return fAxis;
 159 }
 160
 161 //====================================================================
 162 AliForwardMultDists::AliForwardMultDists()
 163   : AliBaseAODTask(),
 164     fBins(),
 165     fSymmetric(0),
 166     fNegative(0),
 167     fPositive(0),
 168     fMCVertex(0),
 169     fDiag(0),
 170     fForwardCache(0),
 171     fCentralCache(0),
 172     fMCCache(0),
 173     fUsePhiAcc(true),
 174     fIsSelected(false)
 175 {}
 176
 177 //____________________________________________________________________
 178 AliForwardMultDists::AliForwardMultDists(const char* name)
 179   : AliBaseAODTask(name,"AliForwardMultDists"),
 180     fBins(),
 181     fSymmetric(0),
 182     fNegative(0),
 183     fPositive(0),
 184     fMCVertex(0),
 185     fDiag(0),
 186     fForwardCache(0),
 187     fCentralCache(0),
 188     fMCCache(0),
 189     fUsePhiAcc(true),
 190     fIsSelected(false)
 191 {
 192   /*
 193    * User constructor
 194    *
 195    * @param name Name of the task
 196    */
 197 }
 198
 199
 200 //____________________________________________________________________
 201 Bool_t
 202 AliForwardMultDists::Book()
 203 {
 204   /*
 205    * Create output objects - called at start of job in slave
 206    *
 207    */
 208   return true;
 209 }
 210 //____________________________________________________________________
 211 Bool_t
 212 AliForwardMultDists::PreData()
 213 {
 214   /*
 215    * Set-up internal structures on first seen event
 216    *
 217    */
 218   AliAODEvent*       aod     = AliForwardUtil::GetAODEvent(this);
 219   AliAODForwardMult* forward = GetForward(*aod);
 220   const TH2D&        hist    = forward->GetHistogram();
 221   Bool_t             useMC   = GetPrimary(*aod) != 0;
 222
 223   fSums->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
 224   fSums->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
 225
 226   const TAxis* xaxis = hist.GetXaxis();
 227   if (!xaxis->GetXbins() || xaxis->GetXbins()->GetSize() <= 0) {
 228     fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
 229                              xaxis->GetNbins(), xaxis->GetXmin(),
 230                              xaxis->GetXmax());
 231     fCentralCache = new TH1D("centralCache", "Projection of Central data",
 232                              xaxis->GetNbins(), xaxis->GetXmin(),
 233                              xaxis->GetXmax());
 234     if (useMC)
 235       fMCCache = new TH1D("mcCache", "Projection of Mc data",
 236                           xaxis->GetNbins(), xaxis->GetXmin(),
 237                           xaxis->GetXmax());
 238   }
 239   else {
 240     fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
 241                              xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 242     fCentralCache = new TH1D("centralCache", "Projection of Central data",
 243                              xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 244     if (useMC)
 245       fMCCache = new TH1D("mcCache", "Projection of Mc data",
 246                           xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 247   }
 248   fForwardCache->SetDirectory(0);
 249   fForwardCache->GetXaxis()->SetTitle("#eta");
 250   fForwardCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 251   fForwardCache->Sumw2();
 252   fCentralCache->SetDirectory(0);
 253   fCentralCache->GetXaxis()->SetTitle("#eta");
 254   fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 255   fCentralCache->Sumw2();
 256
 257   if (useMC) {
 258     fMCCache->SetDirectory(0);
 259     fMCCache->GetXaxis()->SetTitle("#eta");
 260     fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 261     fMCCache->Sumw2();
 262
 263     fMCVertex = static_cast<TH1*>(fVertex->Clone("mcVertex"));
 264     fMCVertex->SetTitle("Vertex distribution from MC");
 265     fMCVertex->SetDirectory(0);
 266     fMCVertex->SetFillColor(kBlue+2);
 267     fSums->Add(fMCVertex);
 268   }
 269
 270   UShort_t xBase = kTrigger-1;
 271   UShort_t yBase = kAnalysis-1;
 272   fDiag = new TH2I("diagnostics", "Event selection",
 273                    kTriggerVertex-xBase, kTrigger-.5,  kTriggerVertex+.5,
 274                    kTriggerVertex-yBase, kAnalysis-.5, kTriggerVertex+.5);
 275   fDiag->SetDirectory(0);
 276   fDiag->GetXaxis()->SetTitle("Selection");
 277   fDiag->GetXaxis()->SetBinLabel(kTrigger      -xBase, "Trigger");
 278   fDiag->GetXaxis()->SetBinLabel(kVertex       -xBase, "Vertex");
 279   fDiag->GetXaxis()->SetBinLabel(kTriggerVertex-xBase, "Trigger+Vertex");
 280   fDiag->GetYaxis()->SetTitle("Type/MC selection");
 281   fDiag->GetYaxis()->SetBinLabel(kAnalysis     -yBase, "Analysis");
 282   fDiag->GetYaxis()->SetBinLabel(kMC           -yBase, "MC");
 283   fDiag->GetYaxis()->SetBinLabel(kTrigger      -yBase, "MC Trigger");
 284   fDiag->GetYaxis()->SetBinLabel(kVertex       -yBase, "MC Vertex");
 285   fDiag->GetYaxis()->SetBinLabel(kTriggerVertex-yBase, "MC Trigger+Vertex");
 286   fDiag->SetMarkerSize(3);
 287   fDiag->SetMarkerColor(kWhite);
 288   fSums->Add(fDiag);
 289
 290   TIter   next(&fBins);
 291   EtaBin* bin = 0;
 292   while ((bin = static_cast<EtaBin*>(next()))) {
 293     bin->SetupForData(fSums, hist, useMC);
 294   }
 295   return true;
 296 }
 297 //____________________________________________________________________
 298 Bool_t
 299 AliForwardMultDists::Event(AliAODEvent& aod)
 300 {
 301   /*
 302    * Analyse a single event
 303    *
 304    * @param aod Input event
 305    */
 306   AliAODForwardMult* forward = GetForward(aod, false, true);
 307   if (!forward) return false;
 308
 309   AliAODCentralMult* central = GetCentral(aod, false, true);
 310   if (!central) return false;
 311
 312   TH2* primary = GetPrimary(aod);
 313
 314   const TH2& forwardData = forward->GetHistogram();
 315   const TH2& centralData = central->GetHistogram();
 316
 317   Double_t vz    = forward->GetIpZ();
 318   Bool_t   trg   = forward->IsTriggerBits(fTriggerMask);
 319   Bool_t   vtx   = (vz <= fMaxIpZ && vz >= fMinIpZ);
 320   Bool_t   ok    = true;  // Should bins process this event?
 321   Bool_t   mcTrg = false;
 322   Bool_t   mcVtx = false;
 323   // If we have MC inpunt
 324   if (primary) {
 325     // For MC, we need to check if we should process the event for
 326     // MC information or not.
 327     if ((fTriggerMask & (AliAODForwardMult::kV0AND|AliAODForwardMult::kNSD))
 328         && !forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
 329       // Bail out if this event is not MC NSD event
 330       ok = false;
 331     else
 332       mcTrg = true;
 333     Double_t mcVz = primary->GetBinContent(0,0);
 334     fMCVertex->Fill(mcVz);
 335     if (mcVz > fMaxIpZ || mcVz < fMinIpZ)
 336       // Bail out if this event was not in the valid range
 337       ok = false;
 338     else
 339       mcVtx = true;
 340   }
 341   // Fill diagnostics
 342   if (trg) {
 343     fDiag->Fill(kTrigger, kAnalysis);
 344     if (mcTrg)          fDiag->Fill(kTrigger,       kTrigger);
 345     if (mcVtx)          fDiag->Fill(kTrigger,       kVertex);
 346     if (mcTrg && mcVtx) fDiag->Fill(kTrigger,       kTriggerVertex);
 347   }
 348   if (vtx) {
 349     fDiag->Fill(kVertex, kAnalysis);
 350     if (mcTrg)          fDiag->Fill(kVertex,        kTrigger);
 351     if (mcVtx)          fDiag->Fill(kVertex,        kVertex);
 352     if (mcTrg && mcVtx) fDiag->Fill(kVertex,        kTriggerVertex);
 353   }
 354   if (trg && vtx) {
 355     fDiag->Fill(kTriggerVertex, kAnalysis);
 356     if (mcTrg)          fDiag->Fill(kTriggerVertex, kTrigger);
 357     if (mcVtx)          fDiag->Fill(kTriggerVertex, kVertex);
 358     if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kTriggerVertex);
 359   }
 360   if (primary) {
 361     if (mcTrg)          fDiag->Fill(kTrigger,       kMC);
 362     if (mcVtx)          fDiag->Fill(kVertex,        kMC);
 363     if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kMC);
 364   }
 365   if (!fIsSelected && !ok) {
 366     // Printf("Event is neither accepted by analysis nor by MC");
 367     return false;
 368   }
 369
 370   // forward->Print();
 371   // Info("", "Event vz=%f", forward->GetIpZ());
 372
 373   ProjectX(forwardData, *fForwardCache, fUsePhiAcc);
 374   ProjectX(centralData, *fCentralCache, fUsePhiAcc);
 375   ProjectX(primary, fMCCache);
 376
 377   TIter   next(&fBins);
 378   EtaBin* bin = 0;
 379   while ((bin = static_cast<EtaBin*>(next()))) {
 380     bin->Process(*fForwardCache, *fCentralCache,
 381                  forwardData,    centralData,
 382                  fIsSelected,    fMCCache);
 383   }
 384
 385   return true;
 386 }
 387
 388 //_____________________________________________________________________
 389 Bool_t
 390 AliForwardMultDists::CheckEvent(const AliAODForwardMult& fwd)
 391 {
 392   fIsSelected = AliBaseAODTask::CheckEvent(fwd);
 393   // We always return true, so that we can process MC truth;
 394   return true;
 395 }
 396
 397 //____________________________________________________________________
 398 Bool_t
 399 AliForwardMultDists::Finalize()
 400 {
 401   /*
 402    * Called at the end of the final processing of the job on the
 403    * full data set (merged data)
 404    */
 405   fResults->Add(fSums->FindObject("triggers")->Clone());
 406   fResults->Add(fSums->FindObject("status")->Clone());
 407   fResults->Add(fSums->FindObject("diagnostics")->Clone());
 408
 409   THStack* sym    = new THStack("all", "All distributions");
 410   THStack* pos    = new THStack("all", "All distributions");
 411   THStack* neg    = new THStack("all", "All distributions");
 412   THStack* oth    = new THStack("all", "All distributions");
 413   THStack* sta    = 0;
 414   EtaBin*  bin    = 0;
 415   TIter    next(&fBins);
 416   while ((bin = static_cast<EtaBin*>(next()))) {
 417     bin->Terminate(fSums, fResults);
 418
 419     sta = oth;
 420     if      (bin->IsSymmetric()) sta = sym;
 421     else if (bin->IsNegative())  sta = neg;
 422     else if (bin->IsPositive())  sta = pos;
 423
 424     TH1*  hh[] = { bin->fSum, bin->fTruth, bin->fTruthAccepted, 0 };
 425     TH1** ph   = hh;
 426
 427     Int_t idx     = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
 428     if (bin->fTruth) idx /= 2;
 429
 430     Int_t mkrBits = GetIndexMarker(idx);
 431     while (*ph) {
 432       Int_t factor = Int_t(TMath::Power(10, idx));
 433       TH1* h = static_cast<TH1*>((*ph)->Clone());
 434       h->SetDirectory(0);
 435       h->Scale(factor);
 436       h->SetTitle(Form("%s (#times%d)", h->GetTitle(), Int_t(factor)));
 437       h->SetMarkerStyle(GetMarkerStyle(mkrBits));
 438       mkrBits ^= kHollow;
 439
 440       sta->Add(h, "p");
 441       ph++;
 442     }
 443   }
 444   TList* lsym = static_cast<TList*>(fResults->FindObject("symmetric"));
 445   TList* lneg = static_cast<TList*>(fResults->FindObject("negative"));
 446   TList* lpos = static_cast<TList*>(fResults->FindObject("positive"));
 447   TList* loth = static_cast<TList*>(fResults->FindObject("other"));
 448
 449   if (lsym) lsym->Add(sym);
 450   if (lneg) lneg->Add(neg);
 451   if (lpos) lpos->Add(pos);
 452   if (loth) loth->Add(oth);
 453
 454   return true;
 455 }
 456
 457 //____________________________________________________________________
 458 void
 459 AliForwardMultDists::ProjectX(const TH2& input, TH1& cache, Bool_t usePhiAcc)
 460 {
 461   /*
 462    * Project a 2D histogram into a 1D histogram taking care to use
 463    * either the @f$\phi@f$ acceptance stored in the overflow bins, or
 464    * the @f$\eta@f$ coverage stored in the underflow bins.
 465    *
 466    * @param input      2D histogram to project
 467    * @param cache      1D histogram to project into
 468    * @param usePhiAcc  If true, use the @f$\phi@f$ acceptance stored in
 469    * the overflow bins, or if false the @f$\eta@f$ coverage stored in
 470    * the underflow bins.
 471    */
 472   cache.Reset();
 473
 474   Int_t nPhi = input.GetNbinsY();
 475   Int_t nEta = input.GetNbinsX();
 476
 477   for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 478     Double_t phiAcc = input.GetBinContent(iEta, nPhi+1);
 479     Double_t etaCov = input.GetBinContent(iEta, 0);
 480     Double_t factor = usePhiAcc ? phiAcc : etaCov;
 481
 482     if (factor < 1e-3) continue;
 483     Double_t sum    = 0;
 484     Double_t e2sum  = 0;
 485     for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
 486       Double_t c = input.GetBinContent(iEta, iPhi);
 487       Double_t e = input.GetBinError(iEta, iPhi);
 488
 489       sum   += c;
 490       e2sum += e * e;
 491     }
 492     cache.SetBinContent(iEta, factor * sum);
 493     cache.SetBinError(iEta, factor * TMath::Sqrt(e2sum));
 494   }
 495 }
 496 //____________________________________________________________________
 497 void
 498 AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
 499 {
 500   /*
 501    * Project on @f$\eta@f$ axis.  If any of the pointers passed is
 502    * zero, do nothing.
 503    *
 504    * @param input
 505    * @param cache
 506    */
 507   if (!input || !cache) return;
 508   cache->Reset();
 509
 510   Int_t nPhi = input->GetNbinsY();
 511   Int_t nEta = input->GetNbinsX();
 512
 513   for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 514     Double_t sum    = 0;
 515     Double_t e2sum  = 0;
 516     for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
 517       Double_t c = input->GetBinContent(iEta, iPhi);
 518       Double_t e = input->GetBinError(iEta, iPhi);
 519
 520       sum   += c;
 521       e2sum += e * e;
 522     }
 523     cache->SetBinContent(iEta, sum);
 524     cache->SetBinError(iEta, TMath::Sqrt(e2sum));
 525   }
 526 }
 527 //____________________________________________________________________
 528 void
 529 AliForwardMultDists::AddBin(const BinSpec& spec)
 530 {
 531   AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
 532 }
 533
 534 //____________________________________________________________________
 535 void
 536 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
 537                             const TAxis& mAxis)
 538 {
 539   /*
 540    * Add an @f$\eta@f$ bin
 541    *
 542    * @param etaLow Low cut on @f$\eta@f$
 543    * @param etaMax High cut on @f$\eta@f$
 544    */
 545   // Symmetric bin
 546   if (etaMax >= kInvalidEta) {
 547     etaLow = -TMath::Abs(etaLow);
 548     etaMax = +TMath::Abs(etaLow);
 549   }
 550   EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
 551   // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
 552   fBins.Add(bin);
 553 }
 554 //____________________________________________________________________
 555 void
 556 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
 557                             UShort_t nMax, UShort_t nDiv)
 558 {
 559   /*
 560    * Add an @f$\eta@f$ bin
 561    *
 562    * @param etaLow Low cut on @f$\eta@f$
 563    * @param etaMax High cut on @f$\eta@f$
 564    */
 565   // Symmetric bin
 566   if (etaMax >= kInvalidEta) {
 567     etaLow = -TMath::Abs(etaLow);
 568     etaMax = +TMath::Abs(etaLow);
 569   }
 570   TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
 571   EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
 572   // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
 573   fBins.Add(bin);
 574 }
 575 #define PFV(N,VALUE)                                    \
 576   do {                                                  \
 577     AliForwardUtil::PrintName(N);                       \
 578     std::cout << (VALUE) << std::endl; } while(false)
 579 #define PFB(N,FLAG)                             \
 580   do {                                                                  \
 581     AliForwardUtil::PrintName(N);                                       \
 582     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 583   } while(false)
 584
 585 //____________________________________________________________________
 586 void
 587 AliForwardMultDists::Print(Option_t* option) const
 588 {
 589   /*
 590    * Print this task
 591    *
 592    * @param option Not used
 593    */
 594   AliBaseAODTask::Print(option);
 595   gROOT->IncreaseDirLevel();
 596   PFB("Use phi acceptance", fUsePhiAcc);
 597   // AliForwardUtil::PrintName("Bins");
 598
 599   Int_t first = true;
 600   TIter   next(&fBins);
 601   EtaBin* bin = 0;
 602   while ((bin = static_cast<EtaBin*>(next()))) {
 603     PFV(first ? "Bins" : "", bin->GetName());
 604     first = false;
 605   }
 606   gROOT->DecreaseDirLevel();
 607 }
 608
 609 //====================================================================
 610 AliForwardMultDists::EtaBin::EtaBin()
 611   : fName(""),
 612     fMAxis(1,0,0),
 613     fTAxis(1,0,0),
 614     fMinEta(0),
 615     fMaxEta(0),
 616     fMinBin(0),
 617     fMaxBin(0),
 618     fSum(0),
 619     fCorr(0),
 620     fResponse(0),
 621     fTruth(0),
 622     fTruthAccepted(0),
 623     fCoverage(0)
 624 {
 625   /*
 626    * I/O constructor
 627    */
 628 }
 629
 630 //____________________________________________________________________
 631 AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta,
 632                                     const TAxis& mAxis)
 633   : fName(""),
 634     fMAxis(1,0,0),
 635     fTAxis(1,0,0),
 636     fMinEta(minEta),
 637     fMaxEta(maxEta),
 638     fMinBin(0),
 639     fMaxBin(0),
 640     fSum(0),
 641     fCorr(0),
 642     fResponse(0),
 643     fTruth(0),
 644     fTruthAccepted(0),
 645     fCoverage(0)
 646 {
 647   /*
 648    * User constructor
 649    *
 650    * @param minEta Least @f$\eta@f$ to consider
 651    * @param maxEta Largest @f$\eta@f$ to consider
 652    */
 653   fName = TString::Format("%+05.2f_%+05.2f", fMinEta, fMaxEta);
 654   fName.ReplaceAll("-", "m");
 655   fName.ReplaceAll("+", "p");
 656   fName.ReplaceAll(".", "d");
 657
 658   // Copy to other our object
 659   mAxis.Copy(fMAxis);
 660
 661   if (mAxis.GetXbins() && mAxis.GetXbins()->GetArray()) {
 662     const TArrayD& mA  = *(mAxis.GetXbins());
 663     TArrayD        tA(mA.GetSize());
 664     Int_t          j   = 0;
 665     Double_t       min = mA[0];
 666     tA[0]              = min;
 667     for (Int_t i = 1; i < mA.GetSize(); i++) {
 668       Double_t d = mA[i] - min;
 669       if (d < 1)
 670         // Not full integer bin
 671         continue;
 672       tA[j+1] = tA[j] + Int_t(d);
 673       min     = tA[j+1];
 674       j++;
 675     }
 676     fTAxis.Set(j, tA.GetArray());
 677   }
 678   else {
 679     // Rounded down maximum and minimum
 680     Int_t max = Int_t(mAxis.GetXmax());
 681     Int_t min = Int_t(mAxis.GetXmin());
 682     fTAxis.Set((max-min)+1, min-.5, max+.5);
 683   }
 684 }
 685 //____________________________________________________________________
 686 AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
 687   : TObject(o),
 688     fName(o.fName),
 689     fMAxis(o.fMAxis),
 690     fTAxis(o.fTAxis),
 691     fMinEta(o.fMinEta),
 692     fMaxEta(o.fMaxEta),
 693     fMinBin(o.fMinBin),
 694     fMaxBin(o.fMaxBin),
 695     fSum(o.fSum),
 696     fCorr(o.fCorr),
 697     fResponse(o.fResponse),
 698     fTruth(o.fTruth),
 699     fTruthAccepted(o.fTruthAccepted),
 700     fCoverage(o.fCoverage)
 701 {}
 702 //____________________________________________________________________
 703 AliForwardMultDists::EtaBin&
 704 AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
 705 {
 706   if (&o == this) return *this;
 707
 708   fName          = o.fName;
 709   fMAxis         = o.fMAxis;
 710   fTAxis         = o.fTAxis;
 711   fMinEta        = o.fMinEta;
 712   fMaxEta        = o.fMaxEta;
 713   fMinBin        = o.fMinBin;
 714   fMaxBin        = o.fMaxBin;
 715   fSum           = o.fSum;
 716   fCorr          = o.fCorr;
 717   fResponse      = o.fResponse;
 718   fTruth         = o.fTruth;
 719   fTruthAccepted = o.fTruthAccepted;
 720   fCoverage      = o.fCoverage;
 721
 722   return *this;
 723 }
 724 //____________________________________________________________________
 725 Bool_t
 726 AliForwardMultDists::EtaBin::IsSymmetric() const
 727 {
 728   return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
 729 }
 730 //____________________________________________________________________
 731 Bool_t
 732 AliForwardMultDists::EtaBin::IsNegative() const
 733 {
 734   return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
 735 }
 736 //____________________________________________________________________
 737 Bool_t
 738 AliForwardMultDists::EtaBin::IsPositive() const
 739 {
 740   return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
 741 }
 742 //____________________________________________________________________
 743 const char*
 744 AliForwardMultDists::EtaBin::ParentName() const
 745 {
 746   if      (IsSymmetric()) return "symmetric";
 747   else if (IsNegative())  return "negative";
 748   else if (IsPositive())  return "positive";
 749   return "other";
 750 }
 751 //____________________________________________________________________
 752 TList*
 753 AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
 754 {
 755   const char* parent = ParentName();
 756   TObject*    op     = l->FindObject(parent);
 757
 758   if (op) return static_cast<TList*>(op);
 759   if (!create) return 0;
 760
 761   // Info("FindParent", "Parent %s not found in %s, creating and adding",
 762   //      parent, l->GetName());
 763   TList* p = new TList;
 764   p->SetName(parent);
 765   p->SetOwner();
 766   l->Add(p);
 767
 768   TList* lowEdges = new TList;
 769   lowEdges->SetName("lowEdges");
 770   lowEdges->SetOwner();
 771   p->Add(lowEdges);
 772
 773   TList* highEdges = new TList;
 774   highEdges->SetName("highEdges");
 775   highEdges->SetOwner();
 776   p->Add(highEdges);
 777
 778   return p;
 779 }
 780
 781 //____________________________________________________________________
 782 TH1*
 783 AliForwardMultDists::EtaBin::CreateH1(const char*  name,
 784                                       const char*  title,
 785                                       const TAxis& xAxis)
 786 {
 787   TH1* ret = 0;
 788
 789   if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
 790     ret = new TH1D(name,title,xAxis.GetNbins(), xAxis.GetXbins()->GetArray());
 791   else
 792     ret = new TH1D(name,title,xAxis.GetNbins(),xAxis.GetXmin(),xAxis.GetXmax());
 793   return ret;
 794 }
 795 //____________________________________________________________________
 796 TH2*
 797 AliForwardMultDists::EtaBin::CreateH2(const char*  name,
 798                                       const char*  title,
 799                                       const TAxis& xAxis,
 800                                       const TAxis& yAxis)
 801 {
 802   TH2* ret = 0;
 803
 804   if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray()) {
 805     if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
 806       // Variable variable
 807       ret = new TH2D(name,title,
 808                      xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
 809                      yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
 810     else
 811       // variable fixed
 812       ret = new TH2D(name,title,
 813                      xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
 814                      yAxis.GetNbins(),yAxis.GetXmin(),yAxis.GetXmax());
 815   }
 816   else {
 817     if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
 818       // fixed variable
 819       ret = new TH2D(name,title,
 820                      xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
 821                      yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
 822     else
 823       // fixed fixed
 824       ret = new TH2D(name,title,
 825                      xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
 826                      yAxis.GetNbins(), yAxis.GetXmin(), yAxis.GetXmax());
 827   }
 828   return ret;
 829 }
 830
 831 //____________________________________________________________________
 832 void
 833 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
 834                                           Bool_t useMC)
 835 {
 836   /*
 837    * Set-up internal structures on first event.
 838    *
 839    * @param list  List to add information to
 840    * @param hist  Template histogram
 841    * @param max   Maximum number of particles
 842    * @param useMC Whether to set-up for MC input
 843    */
 844   TList* p = FindParent(list, true);
 845   TList* l = new TList;
 846   l->SetName(GetName());
 847   l->SetOwner();
 848   p->Add(l);
 849   TList* le = static_cast<TList*>(p->FindObject("lowEdges"));
 850   TList* he = static_cast<TList*>(p->FindObject("highEdges"));
 851   if (!le || !he) {
 852     AliError("Failed to get bin edge lists from parent");
 853     return;
 854   }
 855   else {
 856     Int_t n = le->GetEntries();
 857     TParameter<double>* lp =
 858       new TParameter<double>(Form("minEta%02d", n), fMinEta);
 859     TParameter<double>* hp =
 860       new TParameter<double>(Form("maxEta%02d", n), fMaxEta);
 861     lp->SetMergeMode('f');
 862     hp->SetMergeMode('f');
 863     le->Add(lp);
 864     he->Add(hp);
 865   }
 866
 867   const TAxis * xax = hist.GetXaxis();
 868   fMinBin = xax->FindFixBin(fMinEta);
 869   fMaxBin = xax->FindFixBin(fMaxEta-.00001);
 870
 871   TString t(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
 872   fSum = CreateH1("rawDist",Form("Raw P(#it{N}_{ch}) in %s",t.Data()),fMAxis);
 873   fSum->SetDirectory(0);
 874   fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
 875   fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
 876   fSum->SetFillColor(kRed+1);
 877   fSum->SetFillStyle(0);
 878   // fSum->SetOption("hist e");
 879   fSum->SetMarkerStyle(20);
 880   fSum->SetMarkerColor(kRed+1);
 881   fSum->SetLineColor(kBlack);
 882   fSum->Sumw2();
 883   l->Add(fSum);
 884
 885   fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
 886   fCorr->SetDirectory(0);
 887   fCorr->GetXaxis()->SetTitle("Forward");
 888   fCorr->GetYaxis()->SetTitle("Central");
 889   fCorr->SetOption("colz");
 890   l->Add(fCorr);
 891
 892   fCoverage = new TH1D("coverage", "Fraction of bins with coverage",
 893                        101, -.5, 100.5);
 894   fCoverage->SetDirectory(0);
 895   fCoverage->SetXTitle("Fraction of bins [%]");
 896   fCoverage->SetFillStyle(3001);
 897   fCoverage->SetFillColor(kGreen+1);
 898   l->Add(fCoverage);
 899
 900   if (!useMC) return;
 901   fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
 902                        fMAxis, fTAxis);
 903   fResponse->SetDirectory(0);
 904   fResponse->SetYTitle("MC");
 905   fResponse->SetXTitle("Signal");
 906   fResponse->SetOption("colz");
 907   l->Add(fResponse);
 908
 909   fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
 910   fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
 911   fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
 912   fTruth->SetLineColor(kBlack);
 913   fTruth->SetFillColor(kBlue+1);
 914   fTruth->SetFillStyle(0);
 915   fTruth->SetDirectory(0);
 916   /// fTruth->SetOption("");
 917   fTruth->SetMarkerColor(kBlue+1);
 918   fTruth->SetMarkerStyle(24);
 919   fTruth->Sumw2();
 920   l->Add(fTruth);
 921
 922   fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
 923   fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s",
 924                                 t.Data()));
 925   fTruthAccepted->SetLineColor(kGray+2);
 926   fTruthAccepted->SetFillColor(kOrange+2);
 927   fTruthAccepted->SetDirectory(0);
 928   /// fTruth->SetOption("");
 929   fTruthAccepted->SetMarkerColor(kOrange+2);
 930   l->Add(fTruthAccepted);
 931 }
 932 //____________________________________________________________________
 933 void
 934 AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
 935                                      const TH1& sumCentral,
 936                                      const TH2& forward,
 937                                      const TH2& central,
 938                                      Bool_t     accepted,
 939                                      const TH1* mc)
 940 {
 941   /*
 942    * Process a single event
 943    *
 944    * @param sumForward  Projection of forward data
 945    * @param sumCentral  Projection of the central data
 946    * @param forward     The original forward data
 947    * @param central     The original central data
 948    */
 949   if (!mc && !accepted)
 950     // If we're not looking at MC data, and the event wasn't selected,
 951     // we bail out - this check is already done, but we do it again
 952     // for safety.
 953     return;
 954
 955   Double_t sum        = 0;
 956   Double_t e2sum      = 0;
 957   Int_t    covered    = 0;
 958   Double_t fsum       = -1;
 959   Double_t csum       = -1;
 960   Double_t mcsum      = 0;
 961   Double_t mce2sum    = 0;
 962   for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) {
 963     if (mc) {
 964       Double_t cM = mc->GetBinContent(iEta);
 965       Double_t eM = mc->GetBinError(iEta);
 966       mcsum   += cM;
 967       mce2sum += eM * eM;
 968     }
 969     if (!accepted)
 970       // Event wasn't selected, but we still need to get the rest of
 971       // the MC data.
 972       continue;
 973
 974     Double_t cF = sumForward.GetBinContent(iEta);
 975     Double_t eF = sumForward.GetBinError(iEta);
 976     Bool_t   bF = forward.GetBinContent(iEta, 0) > 0;
 977     Double_t cC = sumCentral.GetBinContent(iEta);
 978     Double_t eC = sumCentral.GetBinError(iEta);
 979     Bool_t   bC = central.GetBinContent(iEta, 0) > 0;
 980     Double_t c  = 0;
 981     Double_t e  = 0;
 982
 983     // If we have an overlap - as given by the eta-coverage,
 984     // calculate the mean
 985     if (bF & bC) {
 986       c = (cF + cC) / 2;
 987       e = TMath::Sqrt(eF*eF + eC*eC);
 988       covered++;
 989     }
 990     // Else, if we have coverage from forward, use that
 991     else if (bF) {
 992       c = cF;
 993       e = eF;
 994       covered++;
 995     }
 996     // Else, if we have covrage from central, use that
 997     else if (bC) {
 998       c = cC;
 999       e = eC;
1000       covered++;
1001     }
1002     // Else, we have incomplete coverage
1003
1004     if (bF) {
1005       if (fsum < 0) fsum = 0;
1006       fsum += cF;
1007     }
1008     if (bC) {
1009       if (csum < 0) csum = 0;
1010       csum += cC;
1011     }
1012
1013     sum   += c;
1014     e2sum += e*e;
1015   }
1016
1017   if (accepted) {
1018     // Only update the histograms if the event was triggered.
1019     // Fill with weight
1020     Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
1021     fSum->Fill(sum, w);
1022     fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
1023
1024     Int_t nTotal = fMaxBin - fMinBin + 1;
1025     fCoverage->Fill(100*float(covered) / nTotal);
1026   }
1027
1028   if (mc) {
1029     Double_t w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
1030     if (fTruth) {
1031       fTruth->Fill(mcsum, w);
1032     }
1033     if (accepted) {
1034       if (fResponse)
1035         // Only fill response matrix for accepted events
1036         fResponse->Fill(sum, mcsum);
1037       if (fTruthAccepted)
1038         fTruthAccepted->Fill(mcsum, w);
1039     }
1040   }
1041 }
1042 //____________________________________________________________________
1043 void
1044 AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
1045 {
1046   /*
1047    * Called at the end of the final processing of the job on the
1048    * full data set (merged data)
1049    *
1050    * @param in    Input list
1051    * @param out   Output list
1052    * @param maxN  Maximum number of @f$N_{ch}@f$ to consider
1053    */
1054   TList* ip = FindParent(in, false);
1055   if (!ip) {
1056     AliErrorF("Parent folder %s not found in input", ParentName());
1057     return;
1058   }
1059
1060   TList* i = dynamic_cast<TList*>(ip->FindObject(GetName()));
1061   if (!i) {
1062     AliErrorF("List %s not found in input", GetName());
1063     return;
1064   }
1065
1066   TList* op = FindParent(out, true);
1067   TList* o  = static_cast<TList*>(i->Clone());
1068   o->SetOwner();
1069   op->Add(o);
1070
1071   fSum           = static_cast<TH1*>(o->FindObject("rawDist"));
1072   fTruth         = static_cast<TH1*>(o->FindObject("truth"));
1073   fTruthAccepted = static_cast<TH1*>(o->FindObject("truthAccepted"));
1074
1075   TH1*  hists[] = { fSum, fTruth, fTruthAccepted, 0 };
1076   TH1** phist   = hists;
1077   while (*phist) {
1078     TH1* h = *phist;
1079     if (h) {
1080       Int_t    maxN = h->GetNbinsX();
1081       Double_t intg = h->Integral(1, maxN);
1082       h->Scale(1. / intg, "width");
1083     }
1084     phist++;
1085   }
1086
1087   if (fTruth && fTruthAccepted) {
1088     TH1*  trgVtx  = static_cast<TH1*>(fTruthAccepted->Clone("triggerVertex"));
1089     TString tit(fTruth->GetTitle());
1090     tit.ReplaceAll("True P(#it{N}_{ch})", "C_{trigger,vertex}");
1091     trgVtx->SetTitle(tit);
1092     trgVtx->SetYTitle("P_{MC}(#it{N}_{ch})/P_{MC,accepted}(#it{N}_{ch})");
1093     trgVtx->Divide(fTruth);
1094     trgVtx->SetDirectory(0);
1095     o->Add(trgVtx);
1096   }
1097
1098 }
1099 //====================================================================
1100 //
1101 // EOF
1102 //