-#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDDENSITYCALCULATOR_H
-#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDDENSITYCALCULATOR_H
+// This class calculates the inclusive charged particle density
+// in each for the 5 FMD rings.
+//
+#ifndef ALIFMDDENSITYCALCULATOR_H
+#define ALIFMDDENSITYCALCULATOR_H
+/**
+ * @file AliFMDDensityCalculator.h
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date Wed Mar 23 14:02:09 2011
+ *
+ * @brief
+ *
+ *
+ * @ingroup pwg2_forward_aod
+ */
#include <TNamed.h>
#include <TList.h>
+#include <TArrayI.h>
#include "AliForwardUtil.h"
class AliESDFMD;
class TH2D;
+class TH1D;
+class TProfile;
+class AliFMDCorrELossFit;
/**
* This class calculates the inclusive charged particle density
* - AliFMDDoubleHitCorrection
* - AliFMDDeadCorrection
*
- * @ingroup pwg2_forward_analysis
+ * @ingroup pwg2_forward_algo
+ * @ingroup pwg2_forward_aod
*/
class AliFMDDensityCalculator : public TNamed
{
*
* @return Reference to this object
*/
- AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator&);
+ AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
+ /**
+ * Initialize this sub-algorithm
+ *
+ * @param etaAxis Not used
+ */
+ virtual void Init(const TAxis& etaAxis);
/**
* Do the calculations
*
*/
virtual Bool_t Calculate(const AliESDFMD& fmd,
AliForwardUtil::Histos& hists,
- Int_t vtxBin, Bool_t lowFlux);
+ UShort_t vtxBin, Bool_t lowFlux);
/**
* Scale the histograms to the total number of events
*
+ * @param dir where to put the output
* @param nEvents Number of events
*/
- void ScaleHistograms(TList* dir, Int_t nEvents);
+ virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
/**
* Output diagnostic histograms to directory
*
* @param dir List to write in
*/
- void DefineOutput(TList* dir);
+ virtual void DefineOutput(TList* dir);
+ /**
+ * Set the debug level. The higher the value the more output
+ *
+ * @param dbg Debug level
+ */
+ void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /**
+ * Maximum particle weight to use
+ *
+ * @param m
+ */
+ void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
+ /**
+ * Set whether to use poisson statistics to estimate the
+ * number of particles that has hit within a region. If this is true,
+ * then the average charge particle density is given by
+ * @f[
+ * \lambda = -\log\left(\frac{N_e}{N_t}\right)
+ * @f]
+ * where $N_e$ is the number of strips within the region that has no
+ * hits over threshold, and $N_t$ is the total number of strips in the
+ * region/
+ *
+ * @param u Whether to use poisson statistics to estimate the
+ * number of particles that has hit within a region.
+ */
+ void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
+ /**
+ * Set the lower multiplicity cut. This overrides the setting in
+ * the energy loss fits.
+ *
+ * @param cut Cut to use
+ */
+ void SetMultCut(Double_t cut) { fMultCut = cut; }
+ /**
+ * Get the multiplicity cut. If the user has set fMultCut (via
+ * SetMultCut) then that value is used. If not, then the lower
+ * value of the fit range for the enery loss fits is returned.
+ *
+ * @return Lower cut on multiplicity
+ */
+ Double_t GetMultCut() const;
+ /**
+ * Print information
+ *
+ * @param option Print options
+ * - max Print max weights
+ */
+ void Print(Option_t* option="") const;
protected:
+ /**
+ * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
+ *
+ * @param cor Correction
+ * @param d Detector
+ * @param r Ring
+ * @param iEta Eta bin
+ */
+ Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
+ UShort_t d, Char_t r, Int_t iEta) const;
+
+ /**
+ * Find the max weights and cache them
+ *
+ */
+ void CacheMaxWeights();
+ /**
+ * Find the (cached) maximum weight for FMD<i>dr</i> in
+ * @f$\eta@f$ bin @a iEta
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param iEta Eta bin
+ *
+ * @return max weight or <= 0 in case of problems
+ */
+ Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
+ /**
+ * Find the (cached) maximum weight for FMD<i>dr</i> iat
+ * @f$\eta@f$
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param eta Eta bin
+ *
+ * @return max weight or <= 0 in case of problems
+ */
+ Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
+
/**
* Get the number of particles corresponding to the signal mult
*
*/
virtual Float_t NParticles(Float_t mult,
UShort_t d, Char_t r, UShort_t s, UShort_t t,
- Int_t v, Float_t eta, Bool_t lowFlux) const;
+ UShort_t v, Float_t eta, Bool_t lowFlux) const;
/**
* Get the inverse correction factor. This consist of
*
* @return
*/
virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
- Int_t v, Float_t eta, Bool_t lowFlux) const;
+ UShort_t v, Float_t eta, Bool_t lowFlux) const;
/**
* Get the acceptance correction for strip @a t in an ring of type @a r
*
* @return Inverse acceptance correction
*/
virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
-
+ /**
+ * Generate the acceptance corrections
+ *
+ * @param r Ring to generate for
+ *
+ * @return Newly allocated histogram of acceptance corrections
+ */
+ virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const;
/**
* Internal data structure to keep track of the histograms
*/
* Destructor
*/
~RingHistos();
+ /**
+ * Initialize the object
+ *
+ * @param eAxis
+ */
+ void Init(const TAxis& eAxis);
+ /**
+ * Make output
+ *
+ * @param dir Where to put it
+ */
void Output(TList* dir);
/**
* Scale the histograms to the total number of events
*
+ * @param dir Where the output is
* @param nEvents Number of events
*/
void ScaleHistograms(TList* dir, Int_t nEvents);
+ /**
+ * Create Poisson histograms
+ */
+ void ResetPoissonHistos();
TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
+ TProfile* fEtaVsN; // Average uncorrected Nch vs eta
+ TProfile* fEtaVsM; // Average corrected Nch vs eta
+ TProfile* fCorr; // Average correction vs eta
TH2D* fDensity; // Distribution inclusive Nch
- ClassDef(RingHistos,1);
+ TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
+ TH2D* fTotalStrips; // Total number of strips in a region
+ TH2D* fEmptyStrips; // Total number of strips in a region
+ TH2D* fBasicHits ; // Total number basic hits in a region
+ TH2D* fEmptyVsTotal; // # of empty strips vs total number of strips
+
+
+ ClassDef(RingHistos,3);
};
/**
* Get the ring histogram container
*/
RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
- TList fRingHistos; // List of histogram containers
- Double_t fMultCut; // Low cut on scaled energy loss
+ TList fRingHistos; // List of histogram containers
+ Double_t fMultCut; // Low cut on scaled energy loss
+ TH1D* fSumOfWeights; // Histogram
+ TH1D* fWeightedSum; // Histogram
+ TH1D* fCorrections; // Histogram
+ UShort_t fMaxParticles; // Maximum particle weight to use
+ Bool_t fUsePoisson; // If true, then use poisson statistics
+ TH1D* fAccI; // Acceptance correction for inner rings
+ TH1D* fAccO; // Acceptance correction for outer rings
+ TArrayI fFMD1iMax; // Array of max weights
+ TArrayI fFMD2iMax; // Array of max weights
+ TArrayI fFMD2oMax; // Array of max weights
+ TArrayI fFMD3iMax; // Array of max weights
+ TArrayI fFMD3oMax; // Array of max weights
+ Int_t fDebug; // Debug level
- ClassDef(AliFMDDensityCalculator,1); // Calculate Nch density
+ ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density
};
#endif