1 // This class calculates the inclusive charged particle density
2 // in each for the 5 FMD rings.
4 #ifndef ALIFMDDENSITYCALCULATOR_H
5 #define ALIFMDDENSITYCALCULATOR_H
7 * @file AliFMDDensityCalculator.h
8 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
9 * @date Wed Mar 23 14:02:09 2011
14 * @ingroup pwg2_forward_aod
19 #include "AliForwardUtil.h"
20 #include "AliFMDMultCuts.h"
21 #include "AliPoissonCalculator.h"
26 class AliFMDCorrELossFit;
29 * This class calculates the inclusive charged particle density
30 * in each for the 5 FMD rings.
33 * - AliESDFMD object possibly corrected for sharing
36 * - 5 RingHistos objects - each with a number of vertex dependent
37 * 2D histograms of the inclusive charge particle density
39 * @par Corrections used:
40 * - AliFMDAnaCalibEnergyDistribution
41 * - AliFMDDoubleHitCorrection
42 * - AliFMDDeadCorrection
44 * @ingroup pwg2_forward_algo
45 * @ingroup pwg2_forward_aod
47 class AliFMDDensityCalculator : public TNamed
51 * How to correct for the missing phi coverage at the corners of the
58 /** Correct the calculated number charged particles */
60 /** Correct the energy loss */
66 AliFMDDensityCalculator();
70 * @param name Name of object
72 AliFMDDensityCalculator(const char* name);
76 * @param o Object to copy from
78 AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
82 virtual ~AliFMDDensityCalculator();
84 * Assignement operator
86 * @param o Object to assign from
88 * @return Reference to this object
90 AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
92 * Initialize this sub-algorithm
94 * @param etaAxis Not used
96 virtual void Init(const TAxis& etaAxis);
100 * @param fmd AliESDFMD object (possibly) corrected for sharing
101 * @param hists Histogram cache
102 * @param vtxBin Vertex bin
103 * @param lowFlux Low flux flag.
105 * @return true on successs
107 virtual Bool_t Calculate(const AliESDFMD& fmd,
108 AliForwardUtil::Histos& hists,
109 UShort_t vtxBin, Bool_t lowFlux, Double_t cent=-1);
111 * Scale the histograms to the total number of events
113 * @param dir where to put the output
114 * @param nEvents Number of events
116 virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
118 * Output diagnostic histograms to directory
120 * @param dir List to write in
122 virtual void DefineOutput(TList* dir);
124 * Set the debug level. The higher the value the more output
126 * @param dbg Debug level
128 void SetDebug(Int_t dbg=1) { fDebug = dbg; }
130 * Maximum particle weight to use
134 void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
136 * Set whether to use poisson statistics to estimate the
137 * number of particles that has hit within a region. If this is true,
138 * then the average charge particle density is given by
140 * \lambda = -\log\left(\frac{N_e}{N_t}\right)
142 * where $N_e$ is the number of strips within the region that has no
143 * hits over threshold, and $N_t$ is the total number of strips in the
146 * @param u Whether to use poisson statistics to estimate the
147 * number of particles that has hit within a region.
149 void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
151 * Set whether to use the phi acceptance correction.
153 * How the phi acceptance is used depends on the value passed.
154 * - 0: No phi acceptance
155 * - 1: Phi acceptance correction done to estimate of particles
156 * - 2: Phi acceptance correction done to energy deposited
158 * @param u If >0, use the phi acceptance (default is false)
160 void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; }
162 * Set the luming factors used in the Poisson method
164 * @param eta Must be 1 or larger
165 * @param phi Must be 1 or larger
167 void SetLumping(Int_t eta, Int_t phi) {
168 fEtaLumping = (eta < 1 ? 1 : eta);
169 fPhiLumping = (phi < 1 ? 1 : phi);
172 * Get the multiplicity cut. If the user has set fMultCut (via
173 * SetMultCut) then that value is used. If not, then the lower
174 * value of the fit range for the enery loss fits is returned.
176 * @return Lower cut on multiplicity
178 Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
179 Bool_t errors=true) const;
181 * Get the multiplicity cut. If the user has set fMultCut (via
182 * SetMultCut) then that value is used. If not, then the lower
183 * value of the fit range for the enery loss fits is returned.
185 * @return Lower cut on multiplicity
187 Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
188 Bool_t errors=true) const;
192 * @param option Print options
193 * - max Print max weights
195 void Print(Option_t* option="") const;
199 * @return Reference to cuts object
201 AliFMDMultCuts& GetCuts() { return fCuts; }
203 * Set the cuts to use
205 * @param c Cuts to use
207 void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
210 * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
212 * @param cor Correction
215 * @param iEta Eta bin
217 Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
218 UShort_t d, Char_t r, Int_t iEta) const;
221 * Find the max weights and cache them
224 void CacheMaxWeights();
226 * Find the (cached) maximum weight for FMD<i>dr</i> in
227 * @f$\eta@f$ bin @a iEta
231 * @param iEta Eta bin
233 * @return max weight or <= 0 in case of problems
235 Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
237 * Find the (cached) maximum weight for FMD<i>dr</i> iat
244 * @return max weight or <= 0 in case of problems
246 Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
249 * Get the number of particles corresponding to the signal mult
255 * @param t Strip (not used)
256 * @param v Vertex bin
257 * @param eta Pseudo-rapidity
258 * @param lowFlux Low-flux flag
260 * @return The number of particles
262 virtual Float_t NParticles(Float_t mult,
263 UShort_t d, Char_t r, UShort_t s, UShort_t t,
264 UShort_t v, Float_t eta, Bool_t lowFlux) const;
266 * Get the inverse correction factor. This consist of
268 * - acceptance correction (corners of sensors)
269 * - double hit correction (for low-flux events)
270 * - dead strip correction
275 * @param t Strip (not used)
276 * @param v Vertex bin
277 * @param eta Pseudo-rapidity
278 * @param lowFlux Low-flux flag
282 virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
283 UShort_t v, Float_t eta, Bool_t lowFlux) const;
285 * Get the acceptance correction for strip @a t in an ring of type @a r
287 * @param r Ring type ('I' or 'O')
288 * @param t Strip number
290 * @return Inverse acceptance correction
292 virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
294 * Generate the acceptance corrections
296 * @param r Ring to generate for
298 * @return Newly allocated histogram of acceptance corrections
300 virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const;
302 * Internal data structure to keep track of the histograms
304 struct RingHistos : public AliForwardUtil::RingHistos
316 RingHistos(UShort_t d, Char_t r);
320 * @param o Object to copy from
322 RingHistos(const RingHistos& o);
324 * Assignment operator
326 * @param o Object to assign from
328 * @return Reference to this
330 RingHistos& operator=(const RingHistos& o);
336 * Initialize the object
340 void Init(const TAxis& eAxis);
344 * @param dir Where to put it
346 void Output(TList* dir);
348 * Scale the histograms to the total number of events
350 * @param dir Where the output is
351 * @param nEvents Number of events
353 void ScaleHistograms(TList* dir, Int_t nEvents);
354 TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
355 TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
356 TProfile* fEtaVsN; // Average uncorrected Nch vs eta
357 TProfile* fEtaVsM; // Average corrected Nch vs eta
358 TProfile* fCorr; // Average correction vs eta
359 TH2D* fDensity; // Distribution inclusive Nch
360 TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
361 AliPoissonCalculator fPoisson; // Calculate density using Poisson method
362 TH1D* fELoss; // Energy loss as seen by this
363 TH1D* fELossUsed; // Energy loss in strips with signal
364 Double_t fMultCut; // If set, use this
366 ClassDef(RingHistos,6);
369 * Get the ring histogram container
374 * @return Ring histogram container
376 RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
377 TList fRingHistos; // List of histogram containers
378 TH1D* fSumOfWeights; // Histogram
379 TH1D* fWeightedSum; // Histogram
380 TH1D* fCorrections; // Histogram
381 UShort_t fMaxParticles; // Maximum particle weight to use
382 Bool_t fUsePoisson; // If true, then use poisson statistics
383 UShort_t fUsePhiAcceptance; // Whether to correct for corners
384 TH1D* fAccI; // Acceptance correction for inner rings
385 TH1D* fAccO; // Acceptance correction for outer rings
386 TArrayI fFMD1iMax; // Array of max weights
387 TArrayI fFMD2iMax; // Array of max weights
388 TArrayI fFMD2oMax; // Array of max weights
389 TArrayI fFMD3iMax; // Array of max weights
390 TArrayI fFMD3oMax; // Array of max weights
391 TH2D* fMaxWeights; // Histogram of max weights
392 TH2D* fLowCuts; // Histogram of low cuts
393 Int_t fEtaLumping; // How to lump eta bins for Poisson
394 Int_t fPhiLumping; // How to lump phi bins for Poisson
395 Int_t fDebug; // Debug level
396 AliFMDMultCuts fCuts; // Cuts
398 ClassDef(AliFMDDensityCalculator,7); // Calculate Nch density