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