FastCentEstimators.C

Go to the documentation of this file.

#ifndef FASTCENTESTIMATOR_C
#define FASTCENTESTIMATOR_C
#ifndef __CINT__
# include <TObject.h>
# include <TString.h>
# include <TH1.h>
# include <TMath.h>
# include <TCollection.h>
# include <TTree.h>
# include <TParticle.h>
# include <TError.h>
# include "FastShortHeader.C"
# include <TH2.h>
#else
class TH1;
// class TH2;
class TCollection;
class TTree;
class TParticle;
class FastShortHeader;
#endif 

enum {
  kNeutron = 2112,
  kProton = 2212
};

//====================================================================
struct FastCentEstimator : public TObject
{
  TString fName;
  Bool_t fVerbose;
  FastCentEstimator(const char* name="")
    : TObject(), fName(name), fVerbose(false)
  {}
  virtual ~FastCentEstimator() {}
  const char* GetName() const { return fName.Data(); }
  void SetVerbose(Bool_t verb) { fVerbose = verb; }
  virtual void Setup(TCollection* out, TTree* tree, UShort_t sNN,
         Bool_t tgtA, Bool_t projA) = 0;
  virtual void PreEvent() {}
  virtual void ProcessHeader(FastShortHeader&) {}
  virtual void Process(const TParticle* p) = 0;
  virtual void PostEvent() {}
  virtual void Terminate(TCollection* out) = 0;
  virtual void Print(Option_t* option="") const
  {
    Printf("%s: %s", ClassName(), GetName());
  }
  static Double_t Theta(const TParticle* p)
  {
    Double_t pT    = p->Pt();
    Double_t pZ    = p->Pz();
    Double_t theta = TMath::ATan2(pT, pZ);
    return theta;
  }
  static Double_t Eta(const TParticle* p)
  {
    Double_t theta = Theta(p);
    Double_t tanth = TMath::Tan(theta/2);
    if (tanth < 1e-100) return +1e10;
    if (tanth > 1e100)  return -1e10;
    Double_t eta   = -TMath::Log(tanth);
    return eta;
  }
  static Double_t Phi(const TParticle* p)
  {
    Double_t px  = p->Px();
    Double_t py  = p->Py();
    Double_t phi = TMath::ATan2(px, py);
    return phi;
  }
  static Bool_t IsPrimary(const TParticle* p)
  {
    return p->TestBit(BIT(14));
  }
  static Bool_t IsWeakDecay(const TParticle* p)
  {
    return p->TestBit(BIT(15));
  }
  static Bool_t IsCharged(const TParticle* p)
  {
    return p->TestBit(BIT(16));
  }
  ClassDef(FastCentEstimator,1);
};




//____________________________________________________________________
struct Fast1DCentEstimator : public FastCentEstimator
{
  ULong64_t fCache;
  TH1* fHistogram; 
  Bool_t fFromTop; 
  Fast1DCentEstimator(const char* name="")
    : FastCentEstimator(name), fHistogram(0), fFromTop(true)
  {}
  virtual ~Fast1DCentEstimator() {}
  virtual const char* MultSpec() const { return "l"; }
  void Setup(TCollection* l, TTree* tree, UShort_t,
       Bool_t, Bool_t)
  {
    if (fHistogram && l) l->Add(fHistogram);
    TString leaves; leaves.Form("value/%s", MultSpec());
    if (tree) tree->Branch(GetName(), &fCache, leaves.Data());
  }
  virtual void PreEvent()
  {
    fCache = 0;
  }
  virtual void PostEvent()
  {
    if (fVerbose) Info("PostEvent", " Got %lld %s particles",
           fCache, GetName());
    fHistogram->Fill(fCache);
  }
  virtual TH1* GetHistogram(TCollection* l) = 0;
  virtual void Terminate(TCollection* out)
  {
    TH1* h    = GetHistogram(out);
    TH1* cent = static_cast<TH1*>(h->Clone(GetName()));
    cent->SetDirectory(0);
    cent->SetYTitle("Centrality [%]");
    cent->SetTitle(Form("%s mapping", GetName()));
    cent->Reset();
    out->Add(cent);
    
    Int_t    nX         = h->GetNbinsX();
    Double_t total      = h->Integral(1,nX);
    Int_t    dBin       = (fFromTop ? -1 : 1);
    Int_t    end        = (fFromTop ? 0  : nX+1);
    Int_t    start      = (fFromTop ? nX : 1);
    if (fVerbose)
      Info("Teminate", "Integrating %s from bin %d to 1",
     h->GetName(), nX);
    for (Int_t i = start; i != end; i += dBin) {
      Double_t curInt  = (fFromTop ?
        h->Integral(i, start) :
        h->Integral(start,i));
      if (curInt < 0) continue;
      Double_t curCent = curInt / total * 100;
      cent->SetBinContent(i, curCent);
      if (fVerbose)
  Info("Terminate", "Bin %3d -> %9f/%9f -> %5.1f%%",
       i, curInt, total, curCent);
    }    
  }
  virtual void Print(Option_t* option="nh") const
  {
    TString opt(option); opt.ToLower();
    if (opt.Contains("n")) 
      Printf("1D Estimator: %s (%screasing)",
       GetName(),(fFromTop ? "de" : "in"));
    if (opt.Contains("h") && fHistogram) {
      Int_t nBin = fHistogram->GetNbinsX();
      Printf("  %d bins between %f and %f",nBin, 
       fHistogram->GetXaxis()->GetXmin(),
       fHistogram->GetXaxis()->GetXmax());
    }
  }
  ClassDef(Fast1DCentEstimator,1);
};

//____________________________________________________________________
struct BCentEstimator : public Fast1DCentEstimator
{
  // TH2* fBvsC;
  const Double_t fkFactor;
  BCentEstimator()
    : Fast1DCentEstimator("B"),
      // fBvsC(0),
      fkFactor(1000)
  {}
  virtual void Setup(TCollection* out, TTree* tree, UShort_t sNN,
         Bool_t tgtA, Bool_t projA)
  {
    fHistogram = MakeHistogram(sNN, tgtA, projA);
#if 0
    if (fHistogram &&
  fHistogram->GetXaxis()) {
      if (fHistogram->GetXaxis()->GetXbins() &&
    fHistogram->GetXaxis()->GetXbins()->GetArray()) {
  fBvsC = new TH2D("bVsC", "Impact parameter vs Centrality",
       fHistogram->GetXaxis()->GetNbins(),
       fHistogram->GetXaxis()->GetXbins()->GetArray(),
       20, 0, 100);
      }
      else {
  fBvsC = new TH2D("bVsC", "Impact parameter vs Centrality",
       fHistogram->GetXaxis()->GetNbins(),
       fHistogram->GetXaxis()->GetXmin(),
       fHistogram->GetXaxis()->GetXmax(),
       20, 0, 100);
      }
    }
    if (fBvsC) {
      fBvsC->SetDirectory(0);
      fBvsC->SetXTitle("b [fm]");
      fBvsC->SetYTitle("Centrality [%]");
      out->Add(fBvsC);
    }
    else
      Warning("BCentEstimator", "Couldn't make bVsC histogram");
#endif
    Fast1DCentEstimator::Setup(out, tree, sNN, tgtA, projA);
    // if (tree) tree->Branch(GetName(), &fB, "value/D");    
  }
  virtual TH1* GetHistogram(TCollection* l)
  {
    return static_cast<TH1*>(l->FindObject(Form("raw%s",GetName())));
  }
  TH1* MakeHistogram(UShort_t sNN, Bool_t tgtA, Bool_t projA)
  {
    TArrayD cents;
    TArrayD bins; // In B
    if (tgtA && projA) { // Pb-Pb
      if (sNN == 2760) {
  // PbPb @ 2.76TeV only 
  // Updated 4th of November 2014 from 
  // cern.ch/twiki/bin/view/ALICE/CentStudies
  //        #Tables_with_centrality_bins_AN1
  Double_t bs[] = { 0,      1.57,  2.22,  2.71,  3.13,
        3.50,   4.94,  6.05,  6.98,  7.81,
        8.55,   9.23,  9.88, 10.47, 11.04,
        11.58, 12.09, 12.58, 13.05, 13.52,
        13.97, 14.43, 14.96, 15.67, 20.00 };
  Double_t cs[] = { 0.5,   1.5,   2.5,   3.5,   4.5,
        7.5,   12.5,  17.5,  22.5,  27.5,
        32.5,  37.5,  42.5,  47.5,  52.5,
        57.5,  62.5,  67.5,  72.5,  77.5,
        82.5,  87.5,  92.5,  97.5 };
  cents.Set(24,cs);
  bins.Set(25,bs);
      }
      else if (sNN == 5023) {
  // PbPb @ 5.02TeV only 
  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentralityCodeSnippets
  Double_t bs[] = { 0.00, 1.56, 2.22, 2.71, 3.13,
        3.51, 3.84, 4.15, 4.43, 4.71,
        4.96, 6.08, 7.01, 7.84, 8.59,
        9.27, 9.92, 10.5, 11.1, 11.6,
        12.1, 12.6, 13.1, 13.6, 14.0,
        14.5, 15.0, 15.7, 19.6  }; // 29
  Double_t cs[] = { 0.5,   1.5,  2.5,  3.5,  4.5,
        5.5,   6.5,  7.5,  8.5,  9.5,
        12.5, 17.5, 22.5, 27.5, 32.5,
        37.5, 42.5, 47.5, 52.5, 57.5,
        62.5, 67.5, 72.5, 77.5, 82.5,
        87.5, 92.5, 97.5  };    
  cents.Set(28,cs);
  bins.Set(29,bs);
      }
    }
    else if (tgtA || projA) { // p-Pb or Pb-p
      if (sNN == 5023) {
  Double_t cs[] = { 2.5,     7.5,     15.,      30.,
        50.,     70.,     90. };
  Double_t bs[] = { 0,       1.83675, 2.59375,  3.66875,
        5.18625, 6.35475, 7.40225, 13.8577 };
  cents.Set(7, cs);
  bins.Set(8, bs);
      }
    }
    if (bins.GetSize() <= 0 || cents.GetSize() <= 0 )  // Nothing defined
      return 0;
    printf("b bins: ");
    for (Int_t i = 0; i < bins.GetSize(); i++) {
      bins[i] *= fkFactor; // Scale to 1/1000 fm
      printf("%s%7.1f", i!=0 ? "-" : "", bins[i]);
    }
    Printf("");
    TH1* h = new TH1D(Form("raw%s",GetName()), "B to Centrality",
          bins.GetSize()-1, bins.GetArray());
    h->SetDirectory(0);
    h->SetXTitle("b\\hbox{ [10^{3}fm]}");
    h->SetYTitle("c\\hbox{ [\\%]}");
    h->SetBinContent(0,1);
    for (Int_t i = 1; i <= cents.GetSize(); i++) {
      h->SetBinContent(i, cents[i-1]);
    }
    return h;
  }
  void PreEvent() { fCache = 100*fkFactor; };
  void ProcessHeader(FastShortHeader& h)
  {
    // In Ap (EPOS) we have many spec in the target, meaning they will
    // be detected on the C side (p is the projectile, A is the target)
    // if (!fSpectators) return;
    fCache = h.fB * fkFactor;
    Info("", "Cache=%ld", fCache);
    // if (fBvsC) fBvsC->Fill(fCache, h.fC);
  }
  virtual void Process(const TParticle*) {};
  void PostEvent() {};
  virtual void Terminate(TCollection* out)
  {
    TH1* h    = GetHistogram(out);
    if (!h) return;
    TH1* cent = static_cast<TH1*>(h->Clone(GetName()));
    cent->SetDirectory(0);
    cent->SetYTitle("Centrality [%]");
    cent->SetTitle(Form("%s mapping", GetName()));
    out->Add(cent);
    // Scale to number of workers
    Double_t scale = cent->GetBinContent(0);
    cent->Scale(1./scale);
  }
  Double_t GetCentrality(Double_t b) const
  {
    Double_t ret = (fHistogram ?
        fHistogram->GetBinContent(fHistogram->FindBin(b*fkFactor)) :
        200);
    // Info("", "Look-up of %f (%f) -> %5.1f%%", b*fkFactor, b, ret);
    return ret;
  }
  ClassDef(BCentEstimator,2);
};
  
//____________________________________________________________________
struct FastNchCentEstimator : public Fast1DCentEstimator
{
  FastNchCentEstimator(const char* name="")
    : Fast1DCentEstimator(name)
  {}
  virtual ~FastNchCentEstimator() {}
  virtual void Process(const TParticle* p)
  {
    if (!IsCharged(p)) return;
    if (!Accept(p))    return;
    fCache++;
  }
  virtual Bool_t Accept(const TParticle* p) = 0;
  virtual void Print(Option_t* option="nh") const
  {
    TString opt(option); opt.ToLower();
    if (opt.Contains("n")) 
      Printf("1D Nch Estimator: %s (%screasing)",
       GetName(),(fFromTop ? "de" : "in"));
    opt.ReplaceAll("n", "");
    Fast1DCentEstimator::Print(opt);
  }
  ClassDef(FastNchCentEstimator,1);
};

//____________________________________________________________________
struct V0CentEstimator : public FastNchCentEstimator
{
  Short_t  fMode;
  Bool_t   fOnlyPrimary;
  Double_t fAEtaMin;
  Double_t fAEtaMax;
  Double_t fCEtaMin;
  Double_t fCEtaMax;
  V0CentEstimator(Short_t mode=0, Bool_t onlyPrimary=false) 
    : FastNchCentEstimator(Form("%s%s", 
        (mode < 0 ? "V0C" : 
         mode > 0 ? "V0A" : "V0M"), 
        (onlyPrimary ? "P" : ""))),
      fMode(mode),
      fOnlyPrimary(onlyPrimary)        
  {
    fAEtaMin = +2.8;
    fAEtaMax = +5.1;
    fCEtaMin = -3.7;
    fCEtaMax = -1.7;
  }
  void Flip(Bool_t onlySign=true)
  {
    if (!onlySign) {
      std::swap(fAEtaMin, fCEtaMax); // a=[-1.7,*] c=[*,+2.8]
      std::swap(fAEtaMax, fCEtaMin); // a=[*,-3.7] c=[+5.1,*]
    }
    else {
      std::swap(fAEtaMin, fAEtaMax);  // a=[+5.1,+2.8]
      std::swap(fCEtaMin, fCEtaMax);  // c=[-1.7,-3.7]
    }
    //                 onlySign     !onlySign
    fAEtaMin *= -1; // a=[-5.1,*]   a=[+1.7,*]
    fAEtaMax *= -1; // a=[*,-2.8]   a=[*,+3.7]
    fCEtaMin *= -1; // c=[+1.7,*]   c=[-5.1,*]
    fCEtaMax *= -1; // c=[*,+3.7]   c=[*,-2.8]
    Printf("Flipped %s", (onlySign ? "sign" : "acceptance"));
  }
  void Setup(TCollection* l, TTree* tree, UShort_t sNN,
       Bool_t tgtA, Bool_t projA)
  {
    Bool_t  isAA  = (tgtA && projA);
    Bool_t  isPA  = (tgtA ^ projA); // XOR
    UInt_t  max   = (isAA ? 13000 : isPA ? 800 : 300);
    UInt_t  dBin  = (isAA ? 10    : isPA ?   1 :   1);
    Color_t color = (fMode < 0 ? kRed : fMode > 0 ? kBlue : kGreen)+2;
    fHistogram = new TH1D(Form("raw%s",GetName()),
        Form("%s #it{N}_{ch} distribution", GetName()),
        max/dBin, 0, (fMode == 0 ? 2 : 1)*max);
    fHistogram->SetXTitle("#it{N}_{ch}");
    fHistogram->SetYTitle("Raw #it{P}(#it{N}_{ch})");
    fHistogram->SetDirectory(0);
    fHistogram->SetLineColor(color);
    fHistogram->SetFillColor(color);
    fHistogram->SetMarkerColor(color);
    fHistogram->SetMarkerStyle(20);
    fHistogram->SetFillStyle(3002);

    Fast1DCentEstimator::Setup(l, tree, sNN, tgtA, projA);
  }
  Bool_t Accept(const TParticle* p)
  {
    if (fOnlyPrimary && !IsPrimary(p)) return false;
    Double_t eta = Eta(p);
    Bool_t   v0A = ((eta >= fAEtaMin) && (eta <= fAEtaMax));
    Bool_t   v0C = ((eta >= fCEtaMin) && (eta <= fCEtaMax));
    if (fMode < 0) return v0C;
    if (fMode > 0) return v0A;
    return v0A || v0C;
  }
  virtual TH1* GetHistogram(TCollection* l)
  {
    return static_cast<TH1*>(l->FindObject(Form("raw%s",GetName())));
  }
  virtual void Print(Option_t* option="nah") const
  {
    TString opt(option); opt.ToLower();
    if (opt.Contains("n")) 
      Printf("V0 Estimator: %s (%screasing) mode=%d",
       GetName(),(fFromTop ? "de" : "in"), fMode);
    if (opt.Contains("a"))
  Printf("  A: eta=[%f,%f], C: eta=[%f,%f]",
         fAEtaMin, fAEtaMax, fCEtaMin, fCEtaMax);
    opt.ReplaceAll("n", "");
    FastNchCentEstimator::Print(opt);
  }
  ClassDef(V0CentEstimator,1);
};
//____________________________________________________________________
struct RefMultEstimator : public FastNchCentEstimator
{
  Double_t fEtaCut;
  RefMultEstimator(Double_t etaCut=0.8) 
    : FastNchCentEstimator(Form("RefMult%02dd%02d",
        Int_t(etaCut), Int_t(100*etaCut)%100)),
      fEtaCut(etaCut)
  {}
  void Setup(TCollection* l, TTree* tree, UShort_t sNN,
       Bool_t tgtA, Bool_t projA)
  {
    Bool_t  isAA  = (tgtA && projA);
    Bool_t  isPA  = (tgtA ^ projA); // XOR
    UInt_t  max   = (isAA ? 15000 : isPA ? 900 : 200);
    UInt_t  dBin  = (isAA ? 10    : isPA ?   1 :   1);
    Color_t color = kMagenta;
    fHistogram = new TH1D(Form("raw%s",GetName()),
        Form("#it{N}_{ch} |#it{#eta}|<%5.2f distribution",
             fEtaCut), max/dBin, 0, max);
    fHistogram->SetXTitle("#it{N}_{ch}");
    fHistogram->SetYTitle("Raw #it{P}(#it{N}_{ch})");
    fHistogram->SetDirectory(0);
    fHistogram->SetLineColor(color);
    fHistogram->SetFillColor(color);
    fHistogram->SetMarkerColor(color);
    fHistogram->SetMarkerStyle(20);
    fHistogram->SetFillStyle(3002);

    Fast1DCentEstimator::Setup(l, tree, sNN, tgtA, projA);
  }
  Bool_t Accept(const TParticle* p)
  {
    Double_t eta = Eta(p);
    if (!IsPrimary(p)) return false;
    if (TMath::Abs(eta) > fEtaCut) return false;
    return true;
  }
  virtual TH1* GetHistogram(TCollection* l)
  {
    return static_cast<TH1*>(l->FindObject(Form("raw%s",GetName())));
  }
  ClassDef(RefMultEstimator,1);
};
//____________________________________________________________________
struct ZNCentEstimator : public Fast1DCentEstimator
{
  // What to put in 
  Short_t   fMode;
  Bool_t    fNeutrons;
  Bool_t    fSpectators;
  Bool_t    fPrimary;
  Double_t  fMinEta;
  Double_t  fMaxEta;
  Double_t  fMaxPhi;
  ULong64_t fNSpec;
  ZNCentEstimator(Short_t mode=0,
      Bool_t neutrons=true,
      Bool_t spectators=false,
      Bool_t primary=false)
    : Fast1DCentEstimator(Form("Z%c%c%c%c",
             (neutrons   ? 'N' : 'P'),
             (mode < 0   ? 'C' : (mode > 0 ? 'A' : 'M')),
             (spectators ? 'S' : 'E'),
             (primary    ? 'P' : 'A'))),
      fMode(mode),
      fNeutrons(neutrons),
      fSpectators(spectators),
      fPrimary(primary),
      fMinEta(0),
      fMaxEta(0),
      fMaxPhi(-1)
  {
    fFromTop = (fSpectators ? false : true);
    if (fNeutrons) {
      const Double_t zN = 1161.3; // distance to ZN cm
      const Double_t dN = 7;      // Area of ZN
      const Double_t rN = TMath::Sqrt(2*dN*dN); // Radius of ZN
      const Double_t tN = TMath::ATan2(rN,zN); // Largest angle ZN
      fMinEta           = -TMath::Log(TMath::Tan(tN/2));
      fMaxEta           = TMath::Infinity();
      fMaxPhi           = -1;
    }
    else {
      const Double_t  zP  = 1156.3; // distance to ZC cm
      const Double_t  wP  = 20.8; // Width of cal
      const Double_t  hP  = 12; // Width of cal
      const Double_t  oP  = 19;
      const Double_t  rO  = TMath::Sqrt(TMath::Power(oP+wP/2,2)+hP*hP/4);
      const Double_t  rI  = TMath::Sqrt(TMath::Power(oP-wP/2,2)+hP*hP/4);
      const Double_t  tO  = TMath::ATan2(rO, zP);
      const Double_t  tI  = TMath::ATan2(rI, zP);
      fMinEta             = -TMath::Log(TMath::Tan(tO/2));
      fMaxEta             = -TMath::Log(TMath::Tan(tI/2));
      fMaxPhi             = TMath::ATan2(hP/2,oP);
    }
  }
  void Setup(TCollection* l, TTree* tree, UShort_t sNN,
       Bool_t tgtA, Bool_t projA)
  {
    Bool_t  isAA  = (tgtA && projA);
    Bool_t  isPA  = (tgtA ^ projA); // XOR
    UInt_t  max   = (isAA ? 2000 : isPA ? 300 : 30);
    UInt_t  dBin  = (isAA ? 10   : isPA ?   1 :  1);
    Color_t color = (fMode < 0 ? kRed : fMode > 0 ? kBlue : kGreen)+2;

    TString nTxt;
    nTxt.Form("#it{N}_{%s%c} (%s)",
        fSpectators ? "spec," : "", fNeutrons ? 'n' : 'p',
        fPrimary ? "primary" : "all");
    fHistogram = new TH1D(Form("raw%s",GetName()),
        Form("%s %s  distribution", GetName(), nTxt.Data()),
        max/dBin, 0, (fMode == 0 ? 2 : 1)*max);
    fHistogram->SetXTitle(nTxt);
    fHistogram->SetYTitle(Form("Raw #it{P}(%s)", nTxt.Data()));
    fHistogram->SetDirectory(0);
    fHistogram->SetLineColor(color);
    fHistogram->SetFillColor(color);
    fHistogram->SetMarkerColor(color);
    fHistogram->SetMarkerStyle(20);
    fHistogram->SetFillStyle(3002);

    Fast1DCentEstimator::Setup(l, tree, sNN, tgtA, projA);
  }
  virtual void Process(const TParticle* p)
  {
    // Check if we should count emitted 
    if (fSpectators) return;

    // Check if we should count only primaries
    if      (fPrimary && !IsPrimary(p)) return; 
    else if (!fPrimary && p->GetStatusCode()==4) return; 
    // if ((fFlags & kPrimary) != 0 && !IsPrimary(p)) return;

    // Check if this is a spectator 
    // Int_t status = p->GetStatusCode();
    // if ((fFlags & kSpectators) == 0 && (status==13 || status==14)) return;

    // Check particle type 
    Int_t aPdg = TMath::Abs(p->GetPdgCode());
    if (fNeutrons && aPdg != kNeutron)
      // Looking for neutrons 
      return;
    if (!fNeutrons && aPdg != kProton)
      // Looking for protons 
      return;

    // Check side
    Double_t eta = Eta(p);
    if (fMode < 0 && eta > 0) return; // Wrong side
    if (fMode > 0 && eta < 0) return; // Wrong side

    // Check acceptance 
    Double_t aeta = TMath::Abs(eta);
    if (aeta > fMaxEta || aeta < fMinEta) return; // Not acceptance

    if (fMaxPhi > 0) { 
      Double_t phi = TMath::Abs(Phi(p) - (eta < 0 ? 0 : TMath::Pi()));
      if (phi > fMaxPhi) return; // Not acceptance 
    }

    // Decrement
    fCache++;
    // Printf("%s: increment from %d (%f) -> %lld",
    //        GetName(), aPdg, eta, fCache);
  }
  void ProcessHeader(FastShortHeader& h)
  {
    // In Ap (EPOS) we have many spec in the target, meaning they will
    // be detected on the C side (p is the projectile, A is the target)
    // if (!fSpectators) return;
    fCache   = 0;
    fNSpec   = 0;
    Int_t nC = 0;
    Int_t nA = 0;
    if (fNeutrons) {
      nC = h.fNSpecNtgt;
      nA = h.fNSpecNproj;
    }
    else {
      nC = h.fNSpecPtgt;
      nA = h.fNSpecPproj;
    }
    if      (fMode < 0) fNSpec = nC;
    else if (fMode > 0) fNSpec = nA;
    else                fNSpec = nC + nA;
    // Printf("%s: Initial cache value: %d", GetName(), fNSpec);
  }
  virtual void PostEvent()
  {
    // Info(GetName(), "Nspec=%lld Nneu=%lld", fNSpec, fCache);
    if (fSpectators) fCache = fNSpec;
    else {
      if (fNSpec > fCache) {
  // Warning("PostEvent", "Nspec (%d) > Nneu (%d)", fNSpec, fCache);
  fNSpec = fCache;
      }
      if (!fPrimary) fCache -= fNSpec;
    }
    Fast1DCentEstimator::PostEvent();
  }
  virtual TH1* GetHistogram(TCollection* l)
  {
    return static_cast<TH1*>(l->FindObject(Form("raw%s",GetName())));
  }
  virtual void Print(Option_t* option="nah") const
  {
    TString opt(option); opt.ToLower();
    if (opt.Contains("n")) 
      Printf("ZN Estimator: %s (%screasing) mode=%d - %s %s",
       GetName(),(fFromTop ? "de" : "in"), fMode,
       (fSpectators ? "spectator" : "emitted"),
       (fPrimary    ? "primary"   : "all"),
       (fNeutrons   ? "neutrons"  : "protons"));
    if (opt.Contains("a"))
      Printf("  eta=[%f,%f], phi=%f", fMinEta, fMaxEta, fMaxPhi);
    opt.ReplaceAll("n", "");
    Fast1DCentEstimator::Print(opt);
  }
  ClassDef(ZNCentEstimator,1);
};
#endif
//
// EOF
//