#include <TH2D.h>
#include <TBrowser.h>
#include <TROOT.h>
+#include <TMath.h>
+#include <TList.h>
#include <iostream>
+//
+// A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
+// using Poisson statistics.
+//
+// The input is assumed to be binned in @f$(\eta,\varphi)@f$ as
+// described by the 2D histogram passwd to the Reset member function.
+//
+// The data is grouped in to regions as defined by the parameters
+// fEtaLumping and fPhiLumping. The total number of cells and number
+// of empty cells is then calculate in each region. The mean
+// multiplicity over the region is then determined as
+//
+// @f[
+// \lange m\rangle = -\log\left(\frac{e}{t}\right)
+// @f]
+// where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
+// total number of cells in the region. A correction for counting
+// statistics, is then applied
+// @f{eqnarray*}
+// c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
+// &=& \frac{1}{1 - \frac{e}{t}}
+// @f{eqnarray*}
+// and the final number in each cell is then
+// @f[
+// h_i c \lange m\rangle
+// @f]
+// where @f$h_i@f$ is the number of hits in the cell @f$i@f$
+//
+//
+
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator()
: TNamed(),
fPhiLumping(5),
fTotal(0),
fEmpty(0),
- fBasic(0)
-{}
+ fBasic(0),
+ fEmptyVsTotal(0),
+ fMean(0),
+ fOcc(0),
+ fCorr(0)
+{
+ //
+ // CTOR
+ //
+}
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator(const char*)
fPhiLumping(5),
fTotal(0),
fEmpty(0),
- fBasic(0)
-{}
+ fBasic(0),
+ fEmptyVsTotal(0),
+ fMean(0),
+ fOcc(0),
+ fCorr(0)
+{
+ //
+ // CTOR
+ //
+}
+//____________________________________________________________________
+AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
+ : TNamed(o),
+ fEtaLumping(o.fEtaLumping),
+ fPhiLumping(o.fPhiLumping),
+ fTotal(0),
+ fEmpty(0),
+ fBasic(0),
+ fEmptyVsTotal(0),
+ fMean(0),
+ fOcc(0),
+ fCorr(0)
+{
+ Init();
+ Reset(o.fBasic);
+}
+
+//____________________________________________________________________
+AliPoissonCalculator::~AliPoissonCalculator()
+{
+ CleanUp();
+}
+
+//____________________________________________________________________
+void
+AliPoissonCalculator::CleanUp()
+{
+ if (fTotal) delete fTotal; fTotal = 0;
+ if (fEmpty) delete fEmpty; fEmpty = 0;
+ if (fBasic) delete fBasic; fBasic = 0;
+ if (fEmptyVsTotal) delete fEmptyVsTotal; fEmptyVsTotal = 0;
+ if (fMean) delete fMean; fMean = 0;
+ if (fOcc) delete fOcc; fOcc = 0;
+ if (fCorr) delete fCorr; fCorr = 0;
+}
+//____________________________________________________________________
+AliPoissonCalculator&
+AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
+{
+ TNamed::operator=(o);
+ fEtaLumping = o.fEtaLumping;
+ fPhiLumping = o.fPhiLumping;
+ CleanUp();
+ Init(-1,-1);
+ Reset(o.fBasic);
+ return *this;
+}
+
+//____________________________________________________________________
+void
+AliPoissonCalculator::Init(Int_t etaLumping, Int_t phiLumping)
+{
+ //
+ // Initialize
+ //
+ if (etaLumping > 0) SetEtaLumping(etaLumping);
+ if (phiLumping > 0) SetEtaLumping(phiLumping);
+ if (fEmptyVsTotal) return;
+
+ Int_t n = fEtaLumping * fPhiLumping + 1;
+ fEmptyVsTotal = new TH2D("emptyVsTotal",
+ "# of empty # bins vs total # bins",
+ n, -.5, n-.5, n, -.5, n-.5);
+ fEmptyVsTotal->SetXTitle("Total # bins");
+ fEmptyVsTotal->SetYTitle("# empty bins");
+ fEmptyVsTotal->SetZTitle("Correlation");
+ fEmptyVsTotal->SetOption("colz");
+ fEmptyVsTotal->SetDirectory(0);
+
+ n = (fEtaLumping + fPhiLumping);
+ fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
+ 10*n+1, -.1, n+.1);
+ fMean->SetXTitle("#bar{N_{ch}}=#log(empty/total)");
+ fMean->SetYTitle("Events");
+ fMean->SetFillColor(kRed+1);
+ fMean->SetFillStyle(3001);
+ fMean->SetDirectory(0);
+
+ fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
+ 100, 0, 100);
+ fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
+ fOcc->SetYTitle("Events");
+ fOcc->SetFillColor(kBlue+1);
+ fOcc->SetFillStyle(3001);
+ fOcc->SetDirectory(0);
+
+ fCorr = new TH2D("correction", "Correction as function of mean N_{ch}",
+ 10*n+1, -.1, n+.1, 100, 0, 10);
+ fCorr->SetXTitle("#bar{N_{ch}}");
+ fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
+ fCorr->SetZTitle("Events");
+ fCorr->SetOption("colz");
+ fCorr->SetDirectory(0);
+}
+
+//____________________________________________________________________
+void
+AliPoissonCalculator::Output(TList* d)
+{
+ if (!d) return;
+ if (!fEmptyVsTotal) Init();
+ d->Add(fEmptyVsTotal);
+ d->Add(fMean);
+ d->Add(fOcc);
+ d->Add(fCorr);
+}
//____________________________________________________________________
void
AliPoissonCalculator::Reset(const TH2D* base)
{
+ //
+ // Reset histogram
+ //
+ if (!base) return;
if (fBasic && fTotal && fEmpty) {
fBasic->Reset();
fTotal->Reset();
Double_t etaMin = base->GetXaxis()->GetXmin();
Double_t etaMax = base->GetXaxis()->GetXmax();
Int_t nPhiF = base->GetNbinsY();
- Int_t nPhi = nPhiF / phiLumping;
+ Int_t nPhi = nPhiF / fPhiLumping;
Double_t phiMin = base->GetYaxis()->GetXmin();
Double_t phiMax = base->GetYaxis()->GetXmax();
AliPoissonCalculator::Fill(Double_t eta, Double_t phi, Bool_t hit,
Double_t weight)
{
+ //
+ // Fill in an observation
+ //
+ // Parameters:
+ // eta Eta value
+ // phi Phi value
+ // hit True if hit
+ // weight Weight if this
+ //
fTotal->Fill(eta, phi);
if (hit) fBasic->Fill(eta, phi, weight);
else fEmpty->Fill(eta, phi);
}
//____________________________________________________________________
-void
-AliPoissonCalculator::Result(TH2D* output)
+Double_t
+AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
+{
+ if (total <= 0) return 0;
+ if (empty < .001) empty = .001;
+ return -TMath::Log(empty/total);
+}
+//____________________________________________________________________
+Double_t
+AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
+{
+ if (total <= 0) return 0;
+ if (TMath::Abs(empty-total) < .001) empty = total - .001;
+ return 1 / (1 - empty / total);
+}
+
+//____________________________________________________________________
+TH2D*
+AliPoissonCalculator::Result()
{
+ //
+ // Calculate result and store in @a output
+ //
+ // Return:
+ // The result histogram (fBase overwritten)
+ //
for (Int_t ieta = 1; ieta <= fBasic->GetNbinsX(); ieta++) {
for (Int_t iphi = 1; iphi <= fBasic->GetNbinsY(); iphi++) {
Double_t eta = fBasic->GetXaxis()->GetBinCenter(ieta);
Double_t hits = fBasic->GetBinContent(ieta,iphi);
// Mean in region of interest
- Double_t poissonM = (total <= 0 || empty <= 0 ? 0 :
- -TMath::Log(empty / total));
-
- // Note, that given filled=total-empty, and
- //
- // m = -log(empty/total)
- // = -log(1 - filled/total)
- //
- // v = m / (1 - exp(-m))
- // = -total/filled * (log(total-filled)-log(total))
- // = -total / (total-empty) * log(empty/total)
- // = total (log(total)-log(empty)) / (total-empty)
- //
- Double_t poissonV = hits;
- if(poissonM > 0)
- // Correct for counting statistics and weight by counts
- poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
- Double_t poissonE = TMath::Sqrt(hits);
+ Double_t poissonM = CalculateMean(empty, total);
+ Double_t poissonC = CalculateCorrection(empty, total);
+ Double_t poissonV = hits * poissonM * poissonC;
+ Double_t poissonE = TMath::Sqrt(poissonV);
if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
- output->SetBinContent(ieta,iphi,poissonV);
- output->SetBinError(ieta,iphi,poissonE);
+ fBasic->SetBinContent(ieta,iphi,poissonV);
+ fBasic->SetBinError(ieta,iphi,poissonE);
}
}
+ for (Int_t ieta = 1; ieta <= fEmpty->GetNbinsX(); ieta++) {
+ for (Int_t iphi = 1; iphi <= fEmpty->GetNbinsY(); iphi++) {
+ Double_t empty = fEmpty->GetBinContent(ieta, iphi);
+ Double_t total = fTotal->GetBinContent(ieta, iphi);
+ Double_t mean = CalculateMean(empty, total);
+ Double_t corr = CalculateCorrection(empty, total);
+ fEmptyVsTotal->Fill(total, empty);
+ fMean->Fill(mean);
+ fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
+ fCorr->Fill(mean, corr);
+ }
+ }
+ return fBasic;
}
//____________________________________________________________________
void
AliPoissonCalculator::Print(const Option_t*) const
{
+ //
+ // Print information
+ //
+ // Parameters:
+ // option Not used
+ //
char ind[gROOT->GetDirLevel()+3];
for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
ind[gROOT->GetDirLevel()] = '\0';
void
AliPoissonCalculator::Browse(TBrowser* b)
{
- if (fBasic) b->Add(fBasic);
- if (fTotal) b->Add(fTotal);
- if (fEmpty) b->Add(fEmpty);
-
+ //
+ // Browse this object
+ //
+ // Parameters:
+ // b Object to browse
+ //
+ if (fTotal) b->Add(fTotal);
+ if (fEmpty) b->Add(fEmpty);
+ if (fBasic) b->Add(fBasic);
+ if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
+ if (fMean) b->Add(fMean);
+ if (fOcc) b->Add(fOcc);
+ if (fCorr) b->Add(fCorr);
}
//
// EOF
#define ALIPOISSONCALCULATOR_H
#include <TNamed.h>
class TH2D;
+class TH1D;
class TBrowser;
+class TList;
-
+/**
+ * A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
+ * using Poisson statistics.
+ *
+ * The input is assumed to be binned in @f$(\eta,\varphi)@f$ as
+ * described by the 2D histogram passwd to the Reset member function.
+ *
+ * The data is grouped in to regions as defined by the parameters
+ * fEtaLumping and fPhiLumping. The total number of cells and number
+ * of empty cells is then calculate in each region. The mean
+ * multiplicity over the region is then determined as
+ *
+ * @f[
+ * \langle m\rangle = -\log\left(\frac{e}{t}\right)
+ * @f]
+ * where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
+ * total number of cells in the region. A correction for counting
+ * statistics, is then applied
+ * @f{eqnarray*}{
+ * c &=& \frac{1}{1 - \exp{-\langle m\rangle}}\\ &=&
+ * \frac{1}{1 - \frac{e}{t}}
+ * @f}
+ * and the final number in each cell is then
+ * @f$h_i c \langle m\rangle@f$
+ * where @f$h_i@f$ is the number of hits in the cell @f$i@f$
+ *
+ */
class AliPoissonCalculator : public TNamed
{
public:
+ /**
+ * Constructor
+ */
AliPoissonCalculator();
+ /**
+ * Constructor
+ *
+ */
AliPoissonCalculator(const char*);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliPoissonCalculator(const AliPoissonCalculator& o);
+
+ /**
+ * Destructor
+ */
virtual ~AliPoissonCalculator();
-
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliPoissonCalculator& operator=(const AliPoissonCalculator& o);
+ /**
+ * Set the number of eta bins to group into a region
+ *
+ * @param n Number of eta bins per region
+ */
void SetEtaLumping(UShort_t n) { fEtaLumping = n; } //*MENU*
+ /**
+ * Set the number of phi bins to group into a region
+ *
+ * @param n Number of phi bins per region
+ */
void SetPhiLumping(UShort_t n) { fPhiLumping = n; } //*MENU*
-
+ /**
+ * Intialize this object
+ *
+ * @param etaLumping If larger than 0, set the eta lumping to this
+ * @param phiLumping If larger than 0, set the phi lumping to this
+ */
+ void Init(Int_t etaLumping=-1, Int_t phiLumping=-1);
+ /**
+ * Output stuff to the passed list
+ *
+ * @param d List to add output histograms to
+ */
+ void Output(TList* d);
+ /**
+ * Reset the cache histogram
+ *
+ * @param base Basic histogram
+ */
void Reset(const TH2D* base);
+ /**
+ * Fill in an observation
+ *
+ * @param eta Eta value
+ * @param phi Phi value
+ * @param hit True if hit
+ * @param weight Weight if this
+ */
void Fill(Double_t eta, Double_t phi, Bool_t hit, Double_t weight=1);
- void Result(TH2D* output);
-
- void IsFolder() const { return kTRUE; }
+ /**
+ * Calculate result and store in @a output
+ *
+ * @return The result histogram (fBase overwritten)
+ */
+ TH2D* Result();
+ /**
+ * @return Always true
+ */
+ Bool_t IsFolder() const { return kTRUE; }
+ /**
+ * Print information
+ *
+ * @param option Not used
+ */
void Print(const Option_t* option="") const;
+ /**
+ * Browse this object
+ *
+ * @param b Object to browse
+ */
void Browse(TBrowser* b);
-protected:
- AliPoissonCalculator(const AliPoissonCalculator& o);
- AliPoissonCalculator& operator=(const AliPoissonCalculator& o);
-
- UShort_t fEtaLumping;
- UShort_t fPhiLumping;
- TH2D* fTotal; // Total number of strips in a region
- TH2D* fEmpty; // Total number of strips in a region
- TH2D* fBasic; // Total number basic hits in a region
+ TH2D* GetEmptyVsTotal() const { return fEmptyVsTotal; }
+ TH1D* GetMean() const { return fMean; }
+ TH1D* GetOccupancy() const { return fOcc; }
+ TH2D* GetCorrection() const { return fCorr; }
+protected:
+ /**
+ * Clean up allocated space
+ *
+ */
+ void CleanUp();
+ /**
+ * Calculate the mean
+ *
+ * This is based on the fact that for a Poisson
+ * @f[
+ * P(n;\lambda) = \frac{-\lambda^n e^{-\lambda}}{n!}
+ * @f]
+ * we have the probability for 0 observation
+ * @f[
+ * P(0;\lambda) = e^{-\lambda} = \frac{N_{empty}}{N_{total}}
+ * @f]
+ * and so we get that the mean is the defined region is
+ * @f[
+ * \lambda = -\log\left(\frac{N_{empty}}{N_{total}}\right)
+ * @f]
+ *
+ * Note the boundary conditions
+ * - @f$N_{total}=0 \rightarrow\lambda=0@f$
+ * - @f$N_{empty}<\epsilon\rightarrow N_{empty} = \epsilon@f$
+ *
+ * @param empty Number of empty bins
+ * @param total Total number of bins
+ *
+ * @return The mean in the defined region
+ */
+ Double_t CalculateMean(Double_t empty, Double_t total) const;
+ /**
+ * The mean @f$\lambda@f$ calculated above is not the full story.
+ * In addition it needs to be corrected using the expression
+ * @f[
+ * \frac{1}{1-e^{\lambda}} =
+ * \frac{1}{1-\frac{N_{empty}}{N_{total}}}
+ * @f]
+ *
+ * Note the boundary conditions
+ * - @f$N_{total}=0 \rightarrow\lambda=0@f$
+ * - @f$|N_{total}-N_{empty}|<\epsilon\rightarrow N_{empty} =
+ * N_{total}-\epsilon@f$
+ *
+ * @param empty Number of empty bins
+ * @param total Total number of bins
+ *
+ * @return The correction to the mean.
+ */
+ Double_t CalculateCorrection(Double_t empty, Double_t total) const;
+ UShort_t fEtaLumping; // Grouping of eta bins
+ UShort_t fPhiLumping; // Grouping of phi bins
+ TH2D* fTotal; // Total number of strips in a region
+ TH2D* fEmpty; // Total number of strips in a region
+ TH2D* fBasic; // Total number basic hits in a region
+ TH2D* fEmptyVsTotal; // Empty versus total cells
+ TH1D* fMean; // Mean calculated by poisson method
+ TH1D* fOcc; // Histogram of occupancies
+ TH2D* fCorr; // Correction as a function of mean
ClassDef(AliPoissonCalculator,1) // Calculate N_ch using Poisson
};
+#endif
// Local Variables:
// mode: C++
// End:
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff