[u/mrichter/AliRoot.git] / 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"
#include "AliFMDMultCuts.h"
class AliESDFMD;
class TH2D;
class TH1D;
class TProfile;
class AliFMDCorrELossFit;

/** 
 * This class calculates the inclusive charged particle density
 * in each for the 5 FMD rings. 
 *
 * @par Input:
 *   - AliESDFMD object possibly corrected for sharing
 *
 * @par Output:
 *   - 5 RingHistos objects - each with a number of vertex dependent 
 *     2D histograms of the inclusive charge particle density 
 * 
 * @par Corrections used: 
 *   - AliFMDAnaCalibEnergyDistribution 
 *   - AliFMDDoubleHitCorrection 
 *   - AliFMDDeadCorrection 
 *
 * @ingroup pwg2_forward_algo
 * @ingroup pwg2_forward_aod
 */
class AliFMDDensityCalculator : public TNamed
{
public:
  /** 
   * Constructor 
   */
  AliFMDDensityCalculator();
  /** 
   * Constructor 
   * 
   * @param name Name of object
   */
  AliFMDDensityCalculator(const char* name);
  /** 
   * Copy constructor 
   * 
   * @param o Object to copy from 
   */
  AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
  /** 
   * Destructor 
   */
  virtual ~AliFMDDensityCalculator();
  /** 
   * Assignement operator
   * 
   * @param o Object to assign from 
   * 
   * @return Reference to this object
   */
  AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
  /** 
   * Initialize this sub-algorithm
   * 
   * @param etaAxis Not used 
   */
  virtual void Init(const TAxis& etaAxis);
  /** 
   * Do the calculations 
   * 
   * @param fmd      AliESDFMD object (possibly) corrected for sharing
   * @param hists    Histogram cache
   * @param vtxBin   Vertex bin 
   * @param lowFlux  Low flux flag. 
   * 
   * @return true on successs 
   */
  virtual Bool_t Calculate(const AliESDFMD& fmd, 
			   AliForwardUtil::Histos& hists, 
			   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 
   */
  virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
  /** 
   * Output diagnostic histograms to directory 
   * 
   * @param dir List to write in
   */  
  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 whether to use the phi acceptance correction. 
   * 
   * How the phi acceptance is used depends on the value passed.  
   * - 0:  No phi acceptance 
   * - 1:  Phi acceptance correction done to estimate of particles 
   * - 2:  Phi acceptance correction done to energy deposited 
   *
   * @param u If >0, use the phi acceptance (default is false)
   */
  void SetUsePhiAcceptance(UShort_t u) { fUsePhiAcceptance = 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) { fCuts.SetMultCuts(cut,cut,cut,cut,cut); }
  /** 
   * Set the lower multiplicity cuts 
   * 
   * @param fmd1i Lower mulitplicyt cut for FMD1i
   * @param fmd2i Lower mulitplicyt cut for FMD2i 
   * @param fmd2o Lower mulitplicyt cut for FMD2o 
   * @param fmd3i Lower mulitplicyt cut for FMD3i 
   * @param fmd3o Lower mulitplicyt cut for FMD3o 
   */
  void SetMultCuts(Double_t fmd1i, 
		   Double_t fmd2i, 
		   Double_t fmd2o, 
		   Double_t fmd3i, 
		   Double_t fmd3o); 
  /** 
   * Set the luming factors used in the Poisson method
   * 
   * @param eta Must be 1 or larger 
   * @param phi Must be 1 or larger 
   */
  void SetLumping(Int_t eta, Int_t phi) { 
    fEtaLumping = (eta < 1 ? 1 : eta); 
    fPhiLumping = (phi < 1 ? 1 : phi); 
  }
  /** 
   * Set the number of landau width to subtract from the most probably
   * value to get the low cut.
   * 
   * @param n Number of @f$ \xi@f$ 
   */
  void SetNXi(Double_t nXi) { fCuts.SetNXi(nXi); /* fNXi = nXi;*/ } 
  /** 
   * Whether to include sigma in the number subtracted from the MPV to
   * get the low cut
   * 
   * @param u If true, then low cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$ 
   */
  void SetIncludeSigma(Bool_t u) { fCuts.SetIncludeSigma(u); /*fIncludeSigma = u;*/ }
  /** 
   * 
   * Set the fraction of MPV
   * 
   * @param cut if true cut at fraction of MPV 
   */
  void SetFractionOfMPV(Double_t cut) { fCuts.SetMPVFraction(cut); /*fFractionOfMPV = 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(UShort_t d, Char_t r, Double_t eta, 
		      Bool_t errors=true) const;
  /** 
   * 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(UShort_t d, Char_t r, Int_t ieta, 
		      Bool_t errors=true) const;
  /** 
   * Print information 
   * 
   * @param option Print options 
   *   - max  Print max weights 
   */
  void Print(Option_t* option="") const;
  AliFMDMultCuts& GetCuts() { return fCuts; }
  void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
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
   * 
   * @param mult     Signal
   * @param d        Detector
   * @param r        Ring 
   * @param s        Sector 
   * @param t        Strip (not used)
   * @param v        Vertex bin 
   * @param eta      Pseudo-rapidity 
   * @param lowFlux  Low-flux flag 
   * 
   * @return The number of particles 
   */
  virtual Float_t NParticles(Float_t mult, 
			     UShort_t d, Char_t r, UShort_t s, UShort_t t, 
			     UShort_t v, Float_t eta, Bool_t lowFlux) const;
  /** 
   * Get the inverse correction factor.  This consist of
   * 
   * - acceptance correction (corners of sensors) 
   * - double hit correction (for low-flux events) 
   * - dead strip correction 
   * 
   * @param d        Detector
   * @param r        Ring 
   * @param s        Sector 
   * @param t        Strip (not used)
   * @param v        Vertex bin 
   * @param eta      Pseudo-rapidity 
   * @param lowFlux  Low-flux flag 
   * 
   * @return 
   */
  virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
			     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
   * 
   * @param r  Ring type ('I' or 'O')
   * @param t  Strip number 
   * 
   * @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
   */
  struct RingHistos : public AliForwardUtil::RingHistos
  { 
    /** 
     * Default CTOR
     */
    RingHistos();
    /** 
     * Constructor
     * 
     * @param d detector
     * @param r ring 
     */
    RingHistos(UShort_t d, Char_t r);
    /** 
     * Copy constructor 
     * 
     * @param o Object to copy from 
     */
    RingHistos(const RingHistos& o);
    /** 
     * Assignment operator 
     * 
     * @param o Object to assign from 
     * 
     * @return Reference to this 
     */
    RingHistos& operator=(const RingHistos& o);
    /** 
     * 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(const TH2D* h, Int_t etaLumping, Int_t phiLumping);
    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
    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 
    TH1D*     fELoss;          // Energy loss as seen by this 
    TH1D*     fELossUsed;      // Energy loss in strips with signal 
    Double_t  fMultCut;        // If set, use this
    
    ClassDef(RingHistos,5);
  };
  /** 
   * Get the ring histogram container 
   * 
   * @param d Detector
   * @param r Ring 
   * 
   * @return Ring histogram container 
   */
  RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
  TList    fRingHistos;    //  List of histogram containers
  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 
  UShort_t fUsePhiAcceptance; // Whether to correct for corners 
  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 
  TH2D*    fMaxWeights;    //  Histogram of max weights
  TH2D*    fLowCuts;       //  Histogram of low cuts
  Int_t    fEtaLumping;    //  How to lump eta bins for Poisson 
  Int_t    fPhiLumping;    //  How to lump phi bins for Poisson 
  Int_t    fDebug;         //  Debug level 
  AliFMDMultCuts fCuts; // Cuts

  ClassDef(AliFMDDensityCalculator,6); // Calculate Nch density 
};

#endif
// Local Variables:
//   mode: C++
// End:
Commit	Line	Data
7984e5f7	1	// This class calculates the inclusive charged particle density
	2	// in each for the 5 FMD rings.
	3	//
	4	#ifndef ALIFMDDENSITYCALCULATOR_H
	5	#define ALIFMDDENSITYCALCULATOR_H
ffca499d	6	/**
	7	* @file AliFMDDensityCalculator.h
	8	* @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
	9	* @date Wed Mar 23 14:02:09 2011
	10	*
	11	* @brief
	12	*
	13	*
	14	* @ingroup pwg2_forward_aod
	15	*/
7e4038b5	16	#include <TNamed.h>
7e4038b5	17	#include <TList.h>
1174780f	18	#include <TArrayI.h>
7e4038b5	19	#include "AliForwardUtil.h"
d2638bb7	20	#include "AliFMDMultCuts.h"
7e4038b5	21	class AliESDFMD;
7e4038b5	22	class TH2D;
dd497217	23	class TH1D;
0bd4b00f	24	class TProfile;
1174780f	25	class AliFMDCorrELossFit;
7e4038b5	26
	27	/**
	28	* This class calculates the inclusive charged particle density
	29	* in each for the 5 FMD rings.
	30	*
	31	* @par Input:
	32	* - AliESDFMD object possibly corrected for sharing
	33	*
	34	* @par Output:
	35	* - 5 RingHistos objects - each with a number of vertex dependent
	36	* 2D histograms of the inclusive charge particle density
	37	*
	38	* @par Corrections used:
	39	* - AliFMDAnaCalibEnergyDistribution
	40	* - AliFMDDoubleHitCorrection
	41	* - AliFMDDeadCorrection
	42	*
8eb443e1	43	* @ingroup pwg2_forward_algo
ffca499d	44	* @ingroup pwg2_forward_aod
7e4038b5	45	*/
	46	class AliFMDDensityCalculator : public TNamed
	47	{
	48	public:
	49	/**
	50	* Constructor
	51	*/
	52	AliFMDDensityCalculator();
	53	/**
	54	* Constructor
	55	*
	56	* @param name Name of object
	57	*/
	58	AliFMDDensityCalculator(const char* name);
	59	/**
	60	* Copy constructor
	61	*
	62	* @param o Object to copy from
	63	*/
	64	AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
	65	/**
	66	* Destructor
	67	*/
	68	virtual ~AliFMDDensityCalculator();
	69	/**
	70	* Assignement operator
	71	*
	72	* @param o Object to assign from
	73	*
	74	* @return Reference to this object
	75	*/
7c1a1f1d	76	AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
1174780f	77	/**
	78	* Initialize this sub-algorithm
	79	*
	80	* @param etaAxis Not used
	81	*/
	82	virtual void Init(const TAxis& etaAxis);
7e4038b5	83	/**
	84	* Do the calculations
	85	*
	86	* @param fmd AliESDFMD object (possibly) corrected for sharing
	87	* @param hists Histogram cache
	88	* @param vtxBin Vertex bin
	89	* @param lowFlux Low flux flag.
	90	*
	91	* @return true on successs
	92	*/
	93	virtual Bool_t Calculate(const AliESDFMD& fmd,
	94	AliForwardUtil::Histos& hists,
0bd4b00f	95	UShort_t vtxBin, Bool_t lowFlux);
7e4038b5	96	/**
	97	* Scale the histograms to the total number of events
	98	*
c389303e	99	* @param dir where to put the output
7e4038b5	100	* @param nEvents Number of events
7e4038b5	101	*/
fb3430ac	102	virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
7e4038b5	103	/**
	104	* Output diagnostic histograms to directory
	105	*
	106	* @param dir List to write in
	107	*/
8eb443e1	108	virtual void DefineOutput(TList* dir);
ea3e5d95	109	/**
	110	* Set the debug level. The higher the value the more output
	111	*
	112	* @param dbg Debug level
	113	*/
	114	void SetDebug(Int_t dbg=1) { fDebug = dbg; }
0bd4b00f	115	/**
	116	* Maximum particle weight to use
	117	*
	118	* @param m
	119	*/
	120	void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
9d05ffeb	121	/**
	122	* Set whether to use poisson statistics to estimate the
	123	* number of particles that has hit within a region. If this is true,
	124	* then the average charge particle density is given by
	125	* @f[
ffca499d	126	* \lambda = -\log\left(\frac{N_e}{N_t}\right)
9d05ffeb	127	* @f]
	128	* where $N_e$ is the number of strips within the region that has no
	129	* hits over threshold, and $N_t$ is the total number of strips in the
	130	* region/
	131	*
	132	* @param u Whether to use poisson statistics to estimate the
	133	* number of particles that has hit within a region.
	134	*/
	135	void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
5bb5d1f6	136	/**
	137	* Set whether to use the phi acceptance correction.
	138	*
0082a8fc	139	* How the phi acceptance is used depends on the value passed.
	140	* - 0: No phi acceptance
	141	* - 1: Phi acceptance correction done to estimate of particles
	142	* - 2: Phi acceptance correction done to energy deposited
	143	*
	144	* @param u If >0, use the phi acceptance (default is false)
5bb5d1f6	145	*/
0082a8fc	146	void SetUsePhiAcceptance(UShort_t u) { fUsePhiAcceptance = u; }
0bd4b00f	147	/**
	148	* Set the lower multiplicity cut. This overrides the setting in
	149	* the energy loss fits.
	150	*
	151	* @param cut Cut to use
	152	*/
d2638bb7	153	void SetMultCut(Double_t cut) { fCuts.SetMultCuts(cut,cut,cut,cut,cut); }
5bb5d1f6	154	/**
	155	* Set the lower multiplicity cuts
	156	*
	157	* @param fmd1i Lower mulitplicyt cut for FMD1i
	158	* @param fmd2i Lower mulitplicyt cut for FMD2i
	159	* @param fmd2o Lower mulitplicyt cut for FMD2o
	160	* @param fmd3i Lower mulitplicyt cut for FMD3i
	161	* @param fmd3o Lower mulitplicyt cut for FMD3o
	162	*/
	163	void SetMultCuts(Double_t fmd1i,
	164	Double_t fmd2i,
	165	Double_t fmd2o,
	166	Double_t fmd3i,
	167	Double_t fmd3o);
b6b35c77	168	/**
	169	* Set the luming factors used in the Poisson method
	170	*
	171	* @param eta Must be 1 or larger
	172	* @param phi Must be 1 or larger
	173	*/
	174	void SetLumping(Int_t eta, Int_t phi) {
	175	fEtaLumping = (eta < 1 ? 1 : eta);
	176	fPhiLumping = (phi < 1 ? 1 : phi);
	177	}
5bb5d1f6	178	/**
	179	* Set the number of landau width to subtract from the most probably
	180	* value to get the low cut.
	181	*
	182	* @param n Number of @f$ \xi@f$
	183	*/
d2638bb7	184	void SetNXi(Double_t nXi) { fCuts.SetNXi(nXi); /* fNXi = nXi;*/ }
5bb5d1f6	185	/**
	186	* Whether to include sigma in the number subtracted from the MPV to
	187	* get the low cut
	188	*
	189	* @param u If true, then low cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$
	190	*/
d2638bb7	191	void SetIncludeSigma(Bool_t u) { fCuts.SetIncludeSigma(u); /fIncludeSigma = u;/ }
	192	/**
	193	*
	194	* Set the fraction of MPV
	195	*
	196	* @param cut if true cut at fraction of MPV
	197	*/
	198	void SetFractionOfMPV(Double_t cut) { fCuts.SetMPVFraction(cut); /fFractionOfMPV = cut;/ }
5bb5d1f6	199	/**
	200	* Get the multiplicity cut. If the user has set fMultCut (via
	201	* SetMultCut) then that value is used. If not, then the lower
	202	* value of the fit range for the enery loss fits is returned.
	203	*
	204	* @return Lower cut on multiplicity
	205	*/
	206	Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
	207	Bool_t errors=true) const;
0bd4b00f	208	/**
	209	* Get the multiplicity cut. If the user has set fMultCut (via
	210	* SetMultCut) then that value is used. If not, then the lower
	211	* value of the fit range for the enery loss fits is returned.
	212	*
	213	* @return Lower cut on multiplicity
	214	*/
5bb5d1f6	215	Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
5bb5d1f6	216	Bool_t errors=true) const;
0bd4b00f	217	/**
	218	* Print information
	219	*
1174780f	220	* @param option Print options
1174780f	221	* - max Print max weights
0bd4b00f	222	*/
0bd4b00f	223	void Print(Option_t* option="") const;
d2638bb7	224	AliFMDMultCuts& GetCuts() { return fCuts; }
d2638bb7	225	void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
7e4038b5	226	protected:
1174780f	227	/**
	228	* Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
	229	*
	230	* @param cor Correction
	231	* @param d Detector
	232	* @param r Ring
	233	* @param iEta Eta bin
	234	*/
fb3430ac	235	Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
1174780f	236	UShort_t d, Char_t r, Int_t iEta) const;
	237
	238	/**
	239	* Find the max weights and cache them
	240	*
	241	*/
	242	void CacheMaxWeights();
	243	/**
	244	* Find the (cached) maximum weight for FMD<i>dr</i> in
	245	* @f$\eta@f$ bin @a iEta
	246	*
	247	* @param d Detector
	248	* @param r Ring
	249	* @param iEta Eta bin
	250	*
	251	* @return max weight or <= 0 in case of problems
	252	*/
	253	Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
	254	/**
	255	* Find the (cached) maximum weight for FMD<i>dr</i> iat
	256	* @f$\eta@f$
	257	*
	258	* @param d Detector
	259	* @param r Ring
	260	* @param eta Eta bin
	261	*
	262	* @return max weight or <= 0 in case of problems
	263	*/
	264	Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
	265
7e4038b5	266	/**
	267	* Get the number of particles corresponding to the signal mult
	268	*
	269	* @param mult Signal
	270	* @param d Detector
	271	* @param r Ring
	272	* @param s Sector
	273	* @param t Strip (not used)
	274	* @param v Vertex bin
	275	* @param eta Pseudo-rapidity
	276	* @param lowFlux Low-flux flag
	277	*
	278	* @return The number of particles
	279	*/
	280	virtual Float_t NParticles(Float_t mult,
	281	UShort_t d, Char_t r, UShort_t s, UShort_t t,
0bd4b00f	282	UShort_t v, Float_t eta, Bool_t lowFlux) const;
7e4038b5	283	/**
	284	* Get the inverse correction factor. This consist of
	285	*
	286	* - acceptance correction (corners of sensors)
	287	* - double hit correction (for low-flux events)
	288	* - dead strip correction
	289	*
	290	* @param d Detector
	291	* @param r Ring
	292	* @param s Sector
	293	* @param t Strip (not used)
	294	* @param v Vertex bin
	295	* @param eta Pseudo-rapidity
	296	* @param lowFlux Low-flux flag
	297	*
	298	* @return
	299	*/
	300	virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
0bd4b00f	301	UShort_t v, Float_t eta, Bool_t lowFlux) const;
7e4038b5	302	/**
	303	* Get the acceptance correction for strip @a t in an ring of type @a r
	304	*
	305	* @param r Ring type ('I' or 'O')
	306	* @param t Strip number
	307	*
	308	* @return Inverse acceptance correction
	309	*/
	310	virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
0bd4b00f	311	/**
	312	* Generate the acceptance corrections
	313	*
	314	* @param r Ring to generate for
	315	*
	316	* @return Newly allocated histogram of acceptance corrections
	317	*/
	318	virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const;
7e4038b5	319	/**
	320	* Internal data structure to keep track of the histograms
	321	*/
9d99b0dd	322	struct RingHistos : public AliForwardUtil::RingHistos
7e4038b5	323	{
	324	/**
	325	* Default CTOR
	326	*/
	327	RingHistos();
	328	/**
	329	* Constructor
	330	*
	331	* @param d detector
	332	* @param r ring
	333	*/
	334	RingHistos(UShort_t d, Char_t r);
	335	/**
	336	* Copy constructor
	337	*
	338	* @param o Object to copy from
	339	*/
	340	RingHistos(const RingHistos& o);
	341	/**
	342	* Assignment operator
	343	*
	344	* @param o Object to assign from
	345	*
	346	* @return Reference to this
	347	*/
	348	RingHistos& operator=(const RingHistos& o);
	349	/**
	350	* Destructor
	351	*/
	352	~RingHistos();
9d05ffeb	353	/**
	354	* Initialize the object
	355	*
	356	* @param eAxis
	357	*/
	358	void Init(const TAxis& eAxis);
c389303e	359	/**
	360	* Make output
	361	*
	362	* @param dir Where to put it
	363	*/
7e4038b5	364	void Output(TList* dir);
9d99b0dd	365	/**
	366	* Scale the histograms to the total number of events
	367	*
c389303e	368	* @param dir Where the output is
9d99b0dd	369	* @param nEvents Number of events
	370	*/
	371	void ScaleHistograms(TList* dir, Int_t nEvents);
e308a636	372	/**
	373	* Create Poisson histograms
	374	*/
b6b35c77	375	void ResetPoissonHistos(const TH2D* h, Int_t etaLumping, Int_t phiLumping);
	376	TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
	377	TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
	378	TProfile* fEtaVsN; // Average uncorrected Nch vs eta
	379	TProfile* fEtaVsM; // Average corrected Nch vs eta
	380	TProfile* fCorr; // Average correction vs eta
	381	TH2D* fDensity; // Distribution inclusive Nch
9d05ffeb	382	TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
b6b35c77	383	TH2D* fTotalStrips; // Total number of strips in a region
	384	TH2D* fEmptyStrips; // Total number of strips in a region
	385	TH2D* fBasicHits ; // Total number basic hits in a region
f7cfc454	386	TH2D* fEmptyVsTotal; // # of empty strips vs total number of # strips
	387	TH1D* fELoss; // Energy loss as seen by this
	388	TH1D* fELossUsed; // Energy loss in strips with signal
	389	Double_t fMultCut; // If set, use this
9fde7142	390
f7cfc454	391	ClassDef(RingHistos,5);
7e4038b5	392	};
	393	/**
	394	* Get the ring histogram container
	395	*
	396	* @param d Detector
	397	* @param r Ring
	398	*
	399	* @return Ring histogram container
	400	*/
	401	RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
ea3e5d95	402	TList fRingHistos; // List of histogram containers
1174780f	403	TH1D* fSumOfWeights; // Histogram
	404	TH1D* fWeightedSum; // Histogram
	405	TH1D* fCorrections; // Histogram
0bd4b00f	406	UShort_t fMaxParticles; // Maximum particle weight to use
9d05ffeb	407	Bool_t fUsePoisson; // If true, then use poisson statistics
0082a8fc	408	UShort_t fUsePhiAcceptance; // Whether to correct for corners
0bd4b00f	409	TH1D* fAccI; // Acceptance correction for inner rings
0bd4b00f	410	TH1D* fAccO; // Acceptance correction for outer rings
1174780f	411	TArrayI fFMD1iMax; // Array of max weights
	412	TArrayI fFMD2iMax; // Array of max weights
	413	TArrayI fFMD2oMax; // Array of max weights
	414	TArrayI fFMD3iMax; // Array of max weights
	415	TArrayI fFMD3oMax; // Array of max weights
5bb5d1f6	416	TH2D* fMaxWeights; // Histogram of max weights
5bb5d1f6	417	TH2D* fLowCuts; // Histogram of low cuts
b6b35c77	418	Int_t fEtaLumping; // How to lump eta bins for Poisson
b6b35c77	419	Int_t fPhiLumping; // How to lump phi bins for Poisson
ea3e5d95	420	Int_t fDebug; // Debug level
d2638bb7	421	AliFMDMultCuts fCuts; // Cuts
7e4038b5	422
d2638bb7	423	ClassDef(AliFMDDensityCalculator,6); // Calculate Nch density
7e4038b5	424	};
	425
	426	#endif
	427	// Local Variables:
	428	// mode: C++
	429	// End:
	430