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 <iostream>
  13 #include <iomanip>
  14
  15 //____________________________________________________________________
  16 AliForwardMultDists::AliForwardMultDists()
  17   : AliAnalysisTaskSE(),
  18     fBins(),
  19     fSymmetric(0),
  20     fNegative(0),
  21     fPositive(0),
  22     fList(0),
  23     fTriggers(0),
  24     fStatus(0),
  25     fTriggerMask(0),
  26     fMinIpZ(0),
  27     fMaxIpZ(0),
  28     fFirstEvent(true),
  29     fForwardCache(0),
  30     fCentralCache(0),
  31     fMCCache(0),
  32     fMaxN(200),
  33     fUsePhiAcc(true)
  34 {}
  35
  36 //____________________________________________________________________
  37 AliForwardMultDists::AliForwardMultDists(const char* name)
  38   : AliAnalysisTaskSE(name),
  39     fBins(),
  40     fSymmetric(0),
  41     fNegative(0),
  42     fPositive(0),
  43     fList(0),
  44     fTriggers(0),
  45     fStatus(0),
  46     fTriggerMask(AliAODForwardMult::kV0AND),
  47     fMinIpZ(-4),
  48     fMaxIpZ(4),
  49     fFirstEvent(true),
  50     fForwardCache(0),
  51     fCentralCache(0),
  52     fMCCache(0),
  53     fMaxN(200),
  54     fUsePhiAcc(true)
  55 {
  56   /**
  57    * User constructor
  58    *
  59    * @param name Name of the task
  60    */
  61   DefineOutput(1, TList::Class());
  62   DefineOutput(2, TList::Class());
  63 }
  64
  65 //____________________________________________________________________
  66 AliForwardMultDists::AliForwardMultDists(const AliForwardMultDists& o)
  67   : AliAnalysisTaskSE(o),
  68     fBins(),
  69     fSymmetric(o.fSymmetric),
  70     fNegative(o.fNegative),
  71     fPositive(o.fPositive),
  72     fList(o.fList),
  73     fTriggers(o.fTriggers),
  74     fStatus(o.fStatus),
  75     fTriggerMask(o.fTriggerMask),
  76     fMinIpZ(o.fMinIpZ),
  77     fMaxIpZ(o.fMaxIpZ),
  78     fFirstEvent(o.fFirstEvent),
  79     fForwardCache(o.fForwardCache),
  80     fCentralCache(o.fCentralCache),
  81     fMCCache(o.fMCCache),
  82     fMaxN(o.fMaxN),
  83     fUsePhiAcc(o.fUsePhiAcc)
  84 {}
  85 //____________________________________________________________________
  86 AliForwardMultDists&
  87 AliForwardMultDists::operator=(const AliForwardMultDists& o)
  88 {
  89   if (&o == this) return *this;
  90
  91   fSymmetric            = o.fSymmetric;
  92   fNegative             = o.fNegative;
  93   fPositive             = o.fPositive;
  94   fList                 = o.fList;
  95   fTriggers             = o.fTriggers;
  96   fStatus               = o.fStatus;
  97   fTriggerMask          = o.fTriggerMask;
  98   fMinIpZ               = o.fMinIpZ;
  99   fMaxIpZ               = o.fMaxIpZ;
 100   fFirstEvent           = o.fFirstEvent;
 101   fForwardCache         = o.fForwardCache;
 102   fCentralCache         = o.fCentralCache;
 103   fMCCache              = o.fMCCache;
 104   fMaxN                 = o.fMaxN;
 105   fUsePhiAcc            = o.fUsePhiAcc;
 106
 107   return *this;
 108 }
 109
 110 //____________________________________________________________________
 111 void
 112 AliForwardMultDists::UserCreateOutputObjects()
 113 {
 114   /**
 115    * Create output objects - called at start of job in slave
 116    *
 117    */
 118
 119   fList = new TList;
 120   fList->SetName("myMultAna");
 121   fList->SetOwner();
 122
 123   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers",
 124                                                       fTriggerMask);
 125   fTriggers->SetDirectory(0);
 126
 127   fStatus = AliAODForwardMult::MakeStatusHistogram("status");
 128   fStatus->SetDirectory(0);
 129
 130   fList->Add(fTriggers);
 131   fList->Add(fStatus);
 132
 133   PostData(1, fList);
 134 }
 135 //____________________________________________________________________
 136 void
 137 AliForwardMultDists::UserExec(Option_t* /*option=""*/)
 138 {
 139   /**
 140    * Analyse a single event
 141    *
 142    * @param option Not used
 143    */
 144   AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
 145   if (!aod) return;
 146
 147   TObject* obj = aod->FindListObject("Forward");
 148   if (!obj) {
 149     AliWarning("No forward object found");
 150     return;
 151   }
 152   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 153
 154   obj = aod->FindListObject("CentralClusters");
 155   if (!obj) {
 156     AliWarning("No central object found");
 157     return;
 158   }
 159   AliAODCentralMult* central = static_cast<AliAODCentralMult*>(obj);
 160
 161   TH2* primary = 0;
 162   obj          = aod->FindListObject("primary");
 163   // We should have a forward object at least
 164   if (obj) primary = static_cast<TH2D*>(obj);
 165
 166   const TH2& forwardData = forward->GetHistogram();
 167   const TH2& centralData = central->GetHistogram();
 168
 169   if (fFirstEvent) {
 170     StoreInformation(forward);
 171     SetupForData(forwardData, primary != 0);
 172     fFirstEvent = false;
 173   }
 174
 175   if (!forward->CheckEvent(fTriggerMask, fMinIpZ, fMaxIpZ, 0, 0,
 176                            fTriggers, fStatus)) return;
 177   // forward->Print();
 178   // Info("", "Event vz=%f", forward->GetIpZ());
 179
 180   ProjectX(forwardData, *fForwardCache, fUsePhiAcc);
 181   ProjectX(centralData, *fCentralCache, fUsePhiAcc);
 182   ProjectX(primary, fMCCache);
 183
 184   TIter   next(&fBins);
 185   EtaBin* bin = 0;
 186   while ((bin = static_cast<EtaBin*>(next()))) {
 187     bin->Process(*fForwardCache, *fCentralCache,
 188                  forwardData,    centralData,
 189                  fMCCache);
 190   }
 191
 192   PostData(1, fList);
 193 }
 194 namespace {
 195   /**
 196    * Marker styles
 197    */
 198   enum {
 199     kSolid        = 0x000,
 200     kHollow       = 0x001,
 201     kCircle       = 0x002,
 202     kSquare       = 0x004,
 203     kUpTriangle   = 0x006,
 204     kDownTriangle = 0x008,
 205     kDiamond      = 0x00a,
 206     kCross        = 0x00c,
 207     kStar         = 0x00e
 208   };
 209   /**
 210    * Get a ROOT marker style from bit options
 211    *
 212    * @param bits Bit mask of options
 213    *
 214    * @return ROOT marker style
 215    */
 216   Int_t GetMarkerStyle(UShort_t bits)
 217   {
 218     Int_t  base   = bits & (0xFE);
 219     Bool_t hollow = bits & kHollow;
 220     switch (base) {
 221     case kCircle:       return (hollow ? 24 : 20);
 222     case kSquare:       return (hollow ? 25 : 21);
 223     case kUpTriangle:   return (hollow ? 26 : 22);
 224     case kDownTriangle: return (hollow ? 32 : 23);
 225     case kDiamond:      return (hollow ? 27 : 33);
 226     case kCross:        return (hollow ? 28 : 34);
 227     case kStar:         return (hollow ? 30 : 29);
 228     }
 229     return 1;
 230   }
 231   /**
 232    * Get the marker option bits from a ROOT style
 233    *
 234    * @param style ROOT style
 235    *
 236    * @return Pattern of marker options
 237    */
 238   UShort_t GetMarkerBits(Int_t style)
 239   {
 240     UShort_t bits = 0;
 241     switch (style) {
 242     case 24: case 25: case 26: case 27: case 28: case 30: case 32:
 243       bits |= kHollow; break;
 244     }
 245     switch (style) {
 246     case 20: case 24: bits |= kCircle;       break;
 247     case 21: case 25: bits |= kSquare;       break;
 248     case 22: case 26: bits |= kUpTriangle;   break;
 249     case 23: case 32: bits |= kDownTriangle; break;
 250     case 27: case 33: bits |= kDiamond;      break;
 251     case 28: case 34: bits |= kCross;        break;
 252     case 29: case 30: bits |= kStar;         break;
 253     }
 254     return bits;
 255   }
 256   static Int_t GetIndexMarker(Int_t idx)
 257   {
 258     const UShort_t nMarkers = 7;
 259     UShort_t markers[] = {
 260       kCircle,
 261       kSquare,
 262       kUpTriangle,
 263       kDownTriangle,
 264       kDiamond,
 265       kCross,
 266       kStar
 267     };
 268
 269     return markers[idx % nMarkers];
 270   }
 271   /**
 272    * Flip the 'hollow-ness' of a marker
 273    *
 274    * @param style ROOT Style
 275    *
 276    * @return ROOT style
 277    */
 278   Int_t FlipHollowStyle(Int_t style)
 279   {
 280     UShort_t bits = GetMarkerBits(style);
 281     Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
 282     return ret;
 283   }
 284 }
 285 //____________________________________________________________________
 286 void
 287 AliForwardMultDists::Terminate(Option_t* /*option=""*/)
 288 {
 289   /**
 290    * Called at the end of the final processing of the job on the
 291    * full data set (merged data)
 292    *
 293    *
 294    * @param option Not used
 295    */
 296   TList* out = new TList;
 297   out->SetName("results");
 298   out->SetOwner();
 299
 300   TList* in = dynamic_cast<TList*>(GetOutputData(1));
 301   if (!in) {
 302     AliError("No data connected to slot 1 here");
 303     return;
 304   }
 305   out->Add(in->FindObject("triggers")->Clone());
 306   out->Add(in->FindObject("status")->Clone());
 307
 308   THStack* sym    = new THStack("all", "All distributions");
 309   THStack* pos    = new THStack("all", "All distributions");
 310   THStack* neg    = new THStack("all", "All distributions");
 311   THStack* oth    = new THStack("all", "All distributions");
 312   THStack* sta    = 0;
 313   EtaBin*  bin    = 0;
 314   TIter    next(&fBins);
 315   while ((bin = static_cast<EtaBin*>(next()))) {
 316     bin->Terminate(in, out, fMaxN);
 317
 318     sta = oth;
 319     if      (bin->IsSymmetric()) sta = sym;
 320     else if (bin->IsNegative())  sta = neg;
 321     else if (bin->IsPositive())  sta = pos;
 322
 323     TH1*  hh[] = { bin->fSum, bin->fTruth, 0 };
 324     TH1** ph   = hh;
 325
 326     Int_t idx     = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
 327     if (bin->fTruth) idx /= 2;
 328
 329     Int_t mkrBits = GetIndexMarker(idx);
 330     while (*ph) {
 331       Int_t factor = TMath::Power(10, idx);
 332       TH1* h = static_cast<TH1*>((*ph)->Clone());
 333       h->SetDirectory(0);
 334       h->Scale(factor);
 335       h->SetTitle(Form("%s (#times%d)", h->GetTitle(), Int_t(factor)));
 336       h->SetMarkerStyle(GetMarkerStyle(mkrBits));
 337       mkrBits ^= kHollow;
 338
 339       sta->Add(h, "p");
 340       ph++;
 341     }
 342   }
 343   TList* lsym = static_cast<TList*>(out->FindObject("symmetric"));
 344   TList* lneg = static_cast<TList*>(out->FindObject("negative"));
 345   TList* lpos = static_cast<TList*>(out->FindObject("positive"));
 346   TList* loth = static_cast<TList*>(out->FindObject("other"));
 347
 348   if (lsym) lsym->Add(sym);
 349   if (lneg) lneg->Add(neg);
 350   if (lpos) lpos->Add(pos);
 351   if (loth) loth->Add(oth);
 352
 353   // out->Add(stack);
 354   PostData(2,out);
 355 }
 356 //____________________________________________________________________
 357 void
 358 AliForwardMultDists::StoreInformation(const AliAODForwardMult* aod)
 359 {
 360   fList->Add(AliForwardUtil::MakeParameter("sys", aod->GetSystem()));
 361   fList->Add(AliForwardUtil::MakeParameter("snn", aod->GetSNN()));
 362   fList->Add(AliForwardUtil::MakeParameter("trigger", ULong_t(fTriggerMask)));
 363   fList->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
 364   fList->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
 365   fList->Add(AliForwardUtil::MakeParameter("maxN", UShort_t(fMaxN)));
 366 }
 367
 368 //____________________________________________________________________
 369 void
 370 AliForwardMultDists::SetupForData(const TH2& hist, Bool_t useMC)
 371 {
 372   /**
 373    * Set-up internal structures on first seen event
 374    *
 375    * @param hist Basic histogram template from AOD object
 376    */
 377   TAxis* xaxis = hist.GetXaxis();
 378   if (!xaxis->GetXbins() || xaxis->GetXbins()->GetSize() <= 0) {
 379     fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
 380                              xaxis->GetNbins(), xaxis->GetXmin(),
 381                              xaxis->GetXmax());
 382     fCentralCache = new TH1D("centralCache", "Projection of Central data",
 383                              xaxis->GetNbins(), xaxis->GetXmin(),
 384                              xaxis->GetXmax());
 385     if (useMC)
 386       fMCCache = new TH1D("mcCache", "Projection of Mc data",
 387                           xaxis->GetNbins(), xaxis->GetXmin(),
 388                           xaxis->GetXmax());
 389   }
 390   else {
 391     fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
 392                              xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 393     fCentralCache = new TH1D("centralCache", "Projection of Central data",
 394                              xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 395     if (useMC)
 396       fMCCache = new TH1D("mcCache", "Projection of Mc data",
 397                           xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
 398   }
 399   fForwardCache->SetDirectory(0);
 400   fForwardCache->GetXaxis()->SetTitle("#eta");
 401   fForwardCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 402   fForwardCache->Sumw2();
 403   fCentralCache->SetDirectory(0);
 404   fCentralCache->GetXaxis()->SetTitle("#eta");
 405   fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 406   fCentralCache->Sumw2();
 407   if (useMC) {
 408     fMCCache->SetDirectory(0);
 409     fMCCache->GetXaxis()->SetTitle("#eta");
 410     fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
 411     fMCCache->Sumw2();
 412   }
 413
 414   TIter   next(&fBins);
 415   EtaBin* bin = 0;
 416   while ((bin = static_cast<EtaBin*>(next()))) {
 417     bin->SetupForData(fList, hist, fMaxN, useMC);
 418   }
 419
 420 }
 421 //____________________________________________________________________
 422 void
 423 AliForwardMultDists::ProjectX(const TH2& input, TH1& cache, Bool_t usePhiAcc)
 424 {
 425   /**
 426    * Project a 2D histogram into a 1D histogram taking care to use
 427    * either the @f$\phi2f$ acceptance stored in the overflow bins, or
 428    * the @f$\eta@f$ coverage stored in the underflow bins.
 429    *
 430    * @param input      2D histogram to project
 431    * @param cache      1D histogram to project into
 432    * @param usePhiAcc  If true, use the @f$\phi2f$ acceptance stored in
 433    * the overflow bins, or if false the @f$\eta@f$ coverage stored in
 434    * the underflow bins.
 435    */
 436   cache.Reset();
 437
 438   Int_t nPhi = input.GetNbinsY();
 439   Int_t nEta = input.GetNbinsX();
 440
 441   for (Int_t iEta = 1; iEta <= nEta; iEta++) {
 442     Double_t phiAcc = input.GetBinContent(iEta, nPhi+1);
 443     Double_t etaCov = input.GetBinContent(iEta, 0);
 444     Double_t factor = usePhiAcc ? phiAcc : etaCov;
 445
 446     if (factor < 1e-3) continue;
 447     Double_t sum    = 0;
 448     Double_t e2sum  = 0;
 449     for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
 450       Double_t c = input.GetBinContent(iEta, iPhi);
 451       Double_t e = input.GetBinError(iEta, iPhi);
 452
 453       sum   += c;
 454       e2sum += e * e;
 455     }
 456     cache.SetBinContent(iEta, factor * sum);
 457     cache.SetBinError(iEta, factor * TMath::Sqrt(e2sum));
 458   }
 459 }
 460 //____________________________________________________________________
 461 void
 462 AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
 463 {
 464   /**
 465    * Project on @f$\eta@f$ axis.  If any of the pointers passed is
 466    * zero, do nothing.
 467    *
 468    * @param input
 469    * @param cache
 470    */
 471   if (!input || !cache) return;
 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 sum    = 0;
 479     Double_t e2sum  = 0;
 480     for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
 481       Double_t c = input->GetBinContent(iEta, iPhi);
 482       Double_t e = input->GetBinError(iEta, iPhi);
 483
 484       sum   += c;
 485       e2sum += e * e;
 486     }
 487     cache->SetBinContent(iEta, sum);
 488     cache->SetBinError(iEta, TMath::Sqrt(e2sum));
 489   }
 490 }
 491 //____________________________________________________________________
 492 void
 493 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax)
 494 {
 495   /**
 496    * Add an @f$\eta@f$ bin
 497    *
 498    * @param etaLow Low cut on @f$\eta@f$
 499    * @param etaMax High cut on @f$\eta@f$
 500    */
 501   // Symmetric bin
 502   if (etaMax >= kInvalidEta) {
 503     etaLow = -TMath::Abs(etaLow);
 504     etaMax = +TMath::Abs(etaLow);
 505   }
 506   EtaBin* bin = new EtaBin(etaLow, etaMax);
 507   // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
 508   fBins.Add(bin);
 509 }
 510 //____________________________________________________________________
 511 void
 512 AliForwardMultDists::SetTriggerMask(const char* mask)
 513 {
 514   /**
 515    * Set the trigger mask
 516    *
 517    * @param mask Mask
 518    */
 519   fTriggerMask = AliAODForwardMult::MakeTriggerMask(mask);
 520 }
 521 //____________________________________________________________________
 522 void
 523 AliForwardMultDists::Print(Option_t* /*option=""*/) const
 524 {
 525   /**
 526    * Print this task
 527    *
 528    * @param option Not used
 529    */
 530   std::cout << "Task: " << GetName() << " " << ClassName() << "\n"
 531             << "  Trigger mask:        "
 532             << AliAODForwardMult::GetTriggerString(fTriggerMask) << "\n"
 533             << "  IpZ range:           " << fMinIpZ <<" - "<< fMaxIpZ << "\n"
 534             << "  Max Nch:             " << fMaxN << "\n"
 535             << "  Use phi acceptance:  " << fUsePhiAcc << "\n"
 536             << "  Bins:" << std::endl;
 537   TIter   next(&fBins);
 538   EtaBin* bin = 0;
 539   while ((bin = static_cast<EtaBin*>(next()))) {
 540     std::cout << "   " << bin->GetName() << std::endl;
 541   }
 542 }
 543
 544 //====================================================================
 545 AliForwardMultDists::EtaBin::EtaBin()
 546   : fName(""),
 547     fMinEta(0),
 548     fMaxEta(0),
 549     fMinBin(0),
 550     fMaxBin(0),
 551     fSum(0),
 552     fCorr(0),
 553     fResponse(0),
 554     fTruth(0),
 555     fCoverage(0)
 556 {
 557   /**
 558    * I/O constructor
 559    */
 560 }
 561
 562 //____________________________________________________________________
 563 AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta)
 564   : fName(""),
 565     fMinEta(minEta),
 566     fMaxEta(maxEta),
 567     fMinBin(0),
 568     fMaxBin(0),
 569     fSum(0),
 570     fCorr(0),
 571     fResponse(0),
 572     fTruth(0),
 573     fCoverage(0)
 574 {
 575   /**
 576    * User constructor
 577    *
 578    * @param minEta Least @f$\eta@f$ to consider
 579    * @param maxEta Largest @f$\eta@f$ to consider
 580    */
 581   fName = TString::Format("%+05.2f_%+05.2f", fMinEta, fMaxEta);
 582   fName.ReplaceAll("-", "m");
 583   fName.ReplaceAll("+", "p");
 584   fName.ReplaceAll(".", "d");
 585 }
 586 //____________________________________________________________________
 587 AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
 588   : TObject(o),
 589     fName(o.fName),
 590     fMinEta(o.fMinEta),
 591     fMaxEta(o.fMaxEta),
 592     fMinBin(o.fMinBin),
 593     fMaxBin(o.fMaxBin),
 594     fSum(o.fSum),
 595     fCorr(o.fCorr),
 596     fResponse(o.fResponse),
 597     fTruth(o.fTruth),
 598     fCoverage(o.fCoverage)
 599 {}
 600 //____________________________________________________________________
 601 AliForwardMultDists::EtaBin&
 602 AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
 603 {
 604   if (&o == this) return *this;
 605
 606   fName         = o.fName;
 607   fMinEta       = o.fMinEta;
 608   fMaxEta       = o.fMaxEta;
 609   fMinBin       = o.fMinBin;
 610   fMaxBin       = o.fMaxBin;
 611   fSum          = o.fSum;
 612   fCorr         = o.fCorr;
 613   fResponse     = o.fResponse;
 614   fTruth        = o.fTruth;
 615   fCoverage     = o.fCoverage;
 616
 617   return *this;
 618 }
 619 //____________________________________________________________________
 620 Bool_t
 621 AliForwardMultDists::EtaBin::IsSymmetric() const
 622 {
 623   return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
 624 }
 625 //____________________________________________________________________
 626 Bool_t
 627 AliForwardMultDists::EtaBin::IsNegative() const
 628 {
 629   return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
 630 }
 631 //____________________________________________________________________
 632 Bool_t
 633 AliForwardMultDists::EtaBin::IsPositive() const
 634 {
 635   return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
 636 }
 637 //____________________________________________________________________
 638 const char*
 639 AliForwardMultDists::EtaBin::ParentName() const
 640 {
 641   if      (IsSymmetric()) return "symmetric";
 642   else if (IsNegative())  return "negative";
 643   else if (IsPositive())  return "positive";
 644   return "other";
 645 }
 646 //____________________________________________________________________
 647 TList*
 648 AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
 649 {
 650   const char* parent = ParentName();
 651   TObject*    op     = l->FindObject(parent);
 652
 653   if (op) return static_cast<TList*>(op);
 654   if (!create) return 0;
 655
 656   // Info("FindParent", "Parent %s not found in %s, creating and adding",
 657   //      parent, l->GetName());
 658   TList* p = new TList;
 659   p->SetName(parent);
 660   p->SetOwner();
 661   l->Add(p);
 662
 663   TList* lowEdges = new TList;
 664   lowEdges->SetName("lowEdges");
 665   lowEdges->SetOwner();
 666   p->Add(lowEdges);
 667
 668   TList* highEdges = new TList;
 669   highEdges->SetName("highEdges");
 670   highEdges->SetOwner();
 671   p->Add(highEdges);
 672
 673   return p;
 674 }
 675
 676 //____________________________________________________________________
 677 void
 678 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
 679                                           UShort_t max, Bool_t useMC)
 680 {
 681   /**
 682    * Set-up internal structures on first event.
 683    *
 684    * @param list  List to add information to
 685    * @param hist  Template histogram
 686    * @param max   Maximum number of particles
 687    * @param useMC Whether to set-up for MC input
 688    */
 689   TList* p = FindParent(list, true);
 690   TList* l = new TList;
 691   l->SetName(GetName());
 692   l->SetOwner();
 693   p->Add(l);
 694   TList* le = static_cast<TList*>(p->FindObject("lowEdges"));
 695   TList* he = static_cast<TList*>(p->FindObject("highEdges"));
 696   if (!le || !he) {
 697     AliError("Failed to get bin edge lists from parent");
 698     return;
 699   }
 700   else {
 701     Int_t n = le->GetEntries();
 702     TParameter<double>* lp =
 703       new TParameter<double>(Form("minEta%02d", n), fMinEta);
 704     TParameter<double>* hp =
 705       new TParameter<double>(Form("maxEta%02d", n), fMaxEta);
 706     lp->SetMergeMode('f');
 707     hp->SetMergeMode('f');
 708     le->Add(lp);
 709     he->Add(hp);
 710   }
 711
 712   fMinBin = hist.GetXaxis()->FindBin(fMinEta);
 713   fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
 714
 715   fSum = new TH1D("rawDist",
 716                   Form("Raw P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f",
 717                        fMinEta, fMaxEta),
 718                   max+1, -.5, max+.5);
 719   fSum->SetDirectory(0);
 720   fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
 721   fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
 722   fSum->SetFillColor(kRed+1);
 723   fSum->SetFillStyle(0);
 724   // fSum->SetOption("hist e");
 725   fSum->SetMarkerStyle(20);
 726   fSum->SetMarkerColor(kRed+1);
 727   fSum->SetLineColor(kBlack);
 728   fSum->Sumw2();
 729   l->Add(fSum);
 730
 731   fCorr = new TH2D("corr", "Correlation of SPD and FMD signals",
 732                    max+2, -1.5, max+.5, max+2, -1.5, max+.5);
 733   fCorr->SetDirectory(0);
 734   fCorr->GetXaxis()->SetTitle("Forward");
 735   fCorr->GetYaxis()->SetTitle("Central");
 736   fCorr->SetOption("colz");
 737   l->Add(fCorr);
 738
 739   fCoverage = new TH1D("coverage", "Fraction of bins with coverage",
 740                        101, -.5, 100.5);
 741   fCoverage->SetDirectory(0);
 742   fCoverage->SetXTitle("Fraction of bins [%]");
 743   fCoverage->SetFillStyle(3001);
 744   fCoverage->SetFillColor(kGreen+1);
 745   l->Add(fCoverage);
 746
 747   if (!useMC) return;
 748   fResponse = new TH2D("response", "Reponse matrix",
 749                        max+1, -.5, max+.5, max+1, -.5, max+.5);
 750   fResponse->SetDirectory(0);
 751   fResponse->SetXTitle("MC");
 752   fResponse->SetYTitle("Signal");
 753   fResponse->SetOption("colz");
 754   l->Add(fResponse);
 755
 756   fTruth = static_cast<TH1D*>(fSum->Clone("truth"));
 757   fTruth->SetTitle(Form("True P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f",
 758                         fMinEta, fMaxEta));
 759   fTruth->SetLineColor(kBlack);
 760   fTruth->SetFillColor(kBlue+1);
 761   fTruth->SetFillStyle(0);
 762   fTruth->SetDirectory(0);
 763   /// fTruth->SetOption("");
 764   fTruth->SetMarkerColor(kBlue+1);
 765   fTruth->SetMarkerStyle(24);
 766   // fTruth->Sumw2();
 767   l->Add(fTruth);
 768 }
 769 //____________________________________________________________________
 770 void
 771 AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
 772                                      const TH1& sumCentral,
 773                                      const TH2& forward,
 774                                      const TH2& central,
 775                                      const TH1* mc)
 776 {
 777   /**
 778    * Process a single event
 779    *
 780    * @param sumForward  Projection of forward data
 781    * @param sumCentral  Projection of the central data
 782    * @param forward     The original forward data
 783    * @param central     The original central data
 784    */
 785   Double_t sum        = 0;
 786   Double_t e2sum      = 0;
 787   Int_t    covered    = 0;
 788   Double_t fsum       = -1;
 789   Double_t csum       = -1;
 790   Double_t mcsum      = 0;
 791   Double_t mce2sum    = 0;
 792   for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) {
 793     Double_t cF = sumForward.GetBinContent(iEta);
 794     Double_t eF = sumForward.GetBinError(iEta);
 795     Bool_t   bF = forward.GetBinContent(iEta, 0) > 0;
 796     Double_t cC = sumCentral.GetBinContent(iEta);
 797     Double_t eC = sumCentral.GetBinError(iEta);
 798     Bool_t   bC = central.GetBinContent(iEta, 0) > 0;
 799     Double_t cM = (mc ? mc->GetBinContent(iEta) : 0);
 800     Double_t eM = (mc ? mc->GetBinError(iEta)   : 0);
 801
 802     /*
 803       Info("",
 804       "bF=%d cF=%7.4f+/-%8.5f "
 805       "bC=%d cC=%7.4f+/-%8.5f "
 806       "cM=%7.4f+/-%48.5f",
 807       bF, cF, eF, bC, cC, eC, cM, eM);
 808     */
 809     if (mc) {
 810       mcsum   += cM;
 811       mce2sum += eM * eM;
 812     }
 813
 814     Double_t c  = 0;
 815     Double_t e  = 0;
 816     // If we have an overlap - as given by the eta-coverage,
 817     // calculate the mean
 818     if (bF & bC) {
 819       c = (cF + cC) / 2;
 820       e = TMath::Sqrt(eF*eF + eC*eC);
 821       covered++;
 822     }
 823     // Else, if we have coverage from forward, use that
 824     else if (bF) {
 825       c = cF;
 826       e = eF;
 827       covered++;
 828     }
 829     // Else, if we have covrage from central, use that
 830     else if (bC) {
 831       c = cC;
 832       e = eC;
 833       covered++;
 834     }
 835     // Else, we have incomplete coverage
 836
 837     if (bF) {
 838       if (fsum < 0) fsum = 0;
 839       fsum += cF;
 840     }
 841     if (bC) {
 842       if (csum < 0) csum = 0;
 843       csum += cC;
 844     }
 845
 846     sum   += c;
 847     e2sum += e*e;
 848   }
 849
 850   // Fill with weight
 851   Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
 852   fSum->Fill(sum, w);
 853   fCorr->Fill(fsum, csum);
 854
 855   Int_t nTotal = fMaxBin - fMinBin + 1;
 856   fCoverage->Fill(100*float(covered) / nTotal);
 857   // Info(GetName(), "covered: %3d, nTotal: %3d -> %3d%% mc: %3d",
 858   //      covered, nTotal, int(100*float(covered)/nTotal), int(mcsum));
 859
 860   if (mc) {
 861     if (fTruth) {
 862       w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
 863       fTruth->Fill(mcsum, w);
 864     }
 865     if (fResponse)
 866       fResponse->Fill(mcsum, sum);
 867   }
 868 }
 869 //____________________________________________________________________
 870 void
 871 AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
 872 {
 873   /**
 874    * Called at the end of the final processing of the job on the
 875    * full data set (merged data)
 876    *
 877    * @param in    Input list
 878    * @param out   Output list
 879    * @param maxN  Maximum number of @f$N_{ch}@f$ to consider
 880    */
 881   TList* ip = FindParent(in, false);
 882   if (!ip) {
 883     AliErrorF("Parent folder %s not found in input", ParentName());
 884     return;
 885   }
 886
 887   TList* i = dynamic_cast<TList*>(ip->FindObject(GetName()));
 888   if (!i) {
 889     AliErrorF("List %s not found in input", GetName());
 890     return;
 891   }
 892
 893   TList* op = FindParent(out, true);
 894   TList* o  = static_cast<TList*>(i->Clone());
 895   o->SetOwner();
 896   op->Add(o);
 897
 898   fSum   = static_cast<TH1*>(o->FindObject("rawDist"));
 899   fTruth = static_cast<TH1*>(o->FindObject("truth"));
 900   TH1*  hists[] = { fSum, fTruth, 0 };
 901   TH1** phist   = hists;
 902   while (*phist) {
 903     TH1* h = *phist;
 904     if (h) {
 905       Double_t intg = h->Integral(1, maxN);
 906       h->Scale(1. / intg);
 907     }
 908     phist++;
 909   }
 910 }
 911 //====================================================================
 912 //
 913 // EOF
 914 //