FlowMC/Centrality/CentralityCalculator.h

#ifndef CENTRALITYCALCULATOR_H
#define CENTRALITYCALCULATOR_H

#include <Rtypes.h>
#include <TH1F.h>
#include <TRandom3.h>

#include <iostream>

//#include "../Utils/Funcs.cc"

struct Lim{
  unsigned short L[10];
  unsigned short H[10];
};

// Calculate centrality bins borders
//S----------------------------------------------------------------
Lim CentralityBorders(TH1F* nCh){
  // printf(">>> Centrality file: %s\n", refMultFileName.Data());
  // TFile *cFile = new TFile(refMultFileName.Data(), "READ");
  // TH1F *nCh = (TH1F*) cFile -> Get("hMultPHENIX");

  int    HistLim  = nCh -> GetNbinsX();
  double Integral = nCh -> Integral(1, HistLim);
  double tmpIntegral;
  int    RightLim;
  Lim    lim;
  double efficiency = 0.93;
  double BinWidth   = 10.;
  double Nbins      = 100. / BinWidth;
  double fraction   = (Integral / 100.) * efficiency * BinWidth;

  printf("Integral: %8.0f\n", Integral);
  printf("Centrality bins width: %2.1f, number of bins: %2.0f\n", BinWidth, Nbins);
  printf("Trigger efficiency: %1.2f\n", efficiency);
  printf("Fraction of events in each centrality bin: %8.2f\n", fraction);

  RightLim = HistLim;
  int j = 0;
  for(int i = HistLim; i > 0; --i){
    tmpIntegral = nCh -> Integral(i, RightLim);
    if(tmpIntegral >= fraction){
      printf("[%d,%d]: %8.2f\n", i, RightLim, tmpIntegral);
      lim.H[j] = RightLim;
      lim.L[j] = i;
      ++j;
      RightLim = i - 1;
    }
  }
  // cFile -> Close();
  return lim;
}
//E----------------------------------------------------------------

class CentralityCalculator
{
private:
  TH1F *fHmult;
  static const Int_t NmultiplicityBins = 100;
  Double_t multiplicity_bins[NmultiplicityBins+1];
  TRandom3 *RNG;
  TString  fModel;
  Double_t fEnergy;
  Int_t GetCent(Int_t multiplicity);
  Int_t Integrate(Int_t max_bin, Float_t sum);
  Double_t fEfficiency;
  Lim      lim;

public:
  CentralityCalculator();
  virtual ~CentralityCalculator();

  void  SetHistogram(TH1F *const& _h) { fHmult = _h; }
  void  SetModel(TString _s) { fModel = _s; }
  void  SetEfficiency(Double_t _a) { fEfficiency = _a; }
  void  SetEnergy(Double_t _a) { fEnergy = _a; }
  Int_t GetCentrality(Int_t multiplicity);
  void  Calc();
  void  DrawHistogram(TString outName);
  void  GenerateMacro(TString outName = "");
};

#endif