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 pwglf_forward_aod
20 #include "AliForwardUtil.h"
21 #include "AliFMDMultCuts.h"
22 #include "AliPoissonCalculator.h"
27 class AliFMDCorrELossFit;
30 * This class calculates the inclusive charged particle density
31 * in each for the 5 FMD rings.
34 * - AliESDFMD object possibly corrected for sharing
37 * - 5 RingHistos objects - each with a number of vertex dependent
38 * 2D histograms of the inclusive charge particle density
40 * @par Corrections used:
41 * - AliFMDAnaCalibEnergyDistribution
42 * - AliFMDDoubleHitCorrection
43 * - AliFMDDeadCorrection
45 * @ingroup pwglf_forward_algo
46 * @ingroup pwglf_forward_aod
48 class AliFMDDensityCalculator : public TNamed
52 * How to correct for the missing phi coverage at the corners of the
59 /** Correct the calculated number charged particles */
61 /** Correct the energy loss */
67 static const char* fgkFolderName;
71 AliFMDDensityCalculator();
75 * @param name Name of object
77 AliFMDDensityCalculator(const char* name);
81 * @param o Object to copy from
83 AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
87 virtual ~AliFMDDensityCalculator();
89 * Assignement operator
91 * @param o Object to assign from
93 * @return Reference to this object
95 AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
97 * Initialize this sub-algorithm
99 * @param etaAxis Not used
101 virtual void SetupForData(const TAxis& etaAxis);
103 * Do the calculations
105 * @param fmd AliESDFMD object (possibly) corrected for sharing
106 * @param hists Histogram cache
107 * @param lowFlux Low flux flag.
108 * @param cent Centrality
109 * @param ip Coordinates of interaction point
111 * @return true on successs
113 virtual Bool_t Calculate(const AliESDFMD& fmd,
114 AliForwardUtil::Histos& hists,
117 const TVector3& ip=TVector3(1024,1024,0));
119 * Scale the histograms to the total number of events
121 * @param dir where to put the output
122 * @param output Output list
123 * @param nEvents Number of events
125 virtual void Terminate(const TList* dir, TList* output, Int_t nEvents);
127 * Output diagnostic histograms to directory
129 * @param dir List to write in
131 virtual void CreateOutputObjects(TList* dir);
133 * Set the debug level. The higher the value the more output
135 * @param dbg Debug level
137 void SetDebug(Int_t dbg=1) { fDebug = dbg; }
139 * Maximum particle weight to use
143 void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
145 * Set whether to use poisson statistics to estimate the
146 * number of particles that has hit within a region. If this is true,
147 * then the average charge particle density is given by
149 * \lambda = -\log\left(\frac{N_e}{N_t}\right)
151 * where $N_e$ is the number of strips within the region that has no
152 * hits over threshold, and $N_t$ is the total number of strips in the
155 * @param u Whether to use poisson statistics to estimate the
156 * number of particles that has hit within a region.
158 void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
160 * In case of a displaced vertices recalculate eta and angle correction
162 * @param use recalculate or not
165 void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; }
167 * In case of a displaced vertices recalculate eta and angle correction
169 * @param use recalculate or not
172 void SetRecalculatePhi(Bool_t use) { fRecalculatePhi = use; }
174 * Set whether to use the phi acceptance correction.
176 * How the phi acceptance is used depends on the value passed.
177 * - 0: No phi acceptance
178 * - 1: Phi acceptance correction done to estimate of particles
179 * - 2: Phi acceptance correction done to energy deposited
181 * @param u If >0, use the phi acceptance (default is false)
184 void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; }
186 * Set the luming factors used in the Poisson method
188 * @param eta Must be 1 or larger
189 * @param phi Must be 1 or larger
191 void SetLumping(Int_t eta, Int_t phi) {
192 fEtaLumping = (eta < 1 ? 1 : eta);
193 fPhiLumping = (phi < 1 ? 1 : phi);
196 * Set the minimum quality of the energy loss fits
198 * @param cut Cut value
200 void SetMinQuality(UShort_t cut=8) { fMinQuality = cut; }
202 * Get the multiplicity cut. If the user has set fMultCut (via
203 * SetMultCut) then that value is used. If not, then the lower
204 * value of the fit range for the enery loss fits is returned.
208 * @param eta Psuedo-rapidity
209 * @param errors Factor in errors
211 * @return Lower cut on multiplicity
213 Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
214 Bool_t errors=true) const;
216 * Get the multiplicity cut. If the user has set fMultCut (via
217 * SetMultCut) then that value is used. If not, then the lower
218 * value of the fit range for the enery loss fits is returned.
222 * @param ieta Psuedo-rapidity bin
223 * @param errors Factor in errors
225 * @return Lower cut on multiplicity
227 Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
228 Bool_t errors=true) const;
230 * Set the minimum quality of the energy loss fits
232 * @param cut Cut value
234 UShort_t GetMinQuality() const { return fMinQuality; }
238 * @param option Print options
239 * - max Print max weights
241 void Print(Option_t* option="") const;
245 * @return Reference to cuts object
247 AliFMDMultCuts& GetCuts() { return fCuts; }
249 * Set the cuts to use
251 * @param c Cuts to use
253 void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
256 * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
258 * @param cor Correction
261 * @param iEta Eta bin
263 * @return The maximum weight
265 Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
266 UShort_t d, Char_t r, Int_t iEta) const;
269 * Find the max weights and cache them
271 * @param axis Default @f$\eta@f$ axis from parent task
273 void CacheMaxWeights(const TAxis& axis);
275 * Find the (cached) maximum weight for FMD<i>dr</i> in
276 * @f$\eta@f$ bin @a iEta
280 * @param iEta Eta bin
282 * @return max weight or <= 0 in case of problems
284 Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
286 * Find the (cached) maximum weight for FMD<i>dr</i> iat
293 * @return max weight or <= 0 in case of problems
295 Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
298 * Get the number of particles corresponding to the signal mult
303 * @param eta Pseudo-rapidity
304 * @param lowFlux Low-flux flag
306 * @return The number of particles
308 virtual Float_t NParticles(Float_t mult,
312 Bool_t lowFlux) const;
314 * Get the inverse correction factor. This consist of
316 * - acceptance correction (corners of sensors)
317 * - double hit correction (for low-flux events)
318 * - dead strip correction
323 * @param eta Pseudo-rapidity
324 * @param lowFlux Low-flux flag
326 * @return the correction factor
328 virtual Float_t Correction(UShort_t d, Char_t r, UShort_t t,
329 Float_t eta, Bool_t lowFlux) const;
331 * Get the acceptance correction for strip @a t in an ring of type @a r
333 * @param r Ring type ('I' or 'O')
334 * @param t Strip number
336 * @return Inverse acceptance correction
338 virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
340 * Generate the acceptance corrections
342 * @param r Ring to generate for
344 * @return Newly allocated histogram of acceptance corrections
346 virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const;
348 * Internal data structure to keep track of the histograms
350 struct RingHistos : public AliForwardUtil::RingHistos
362 RingHistos(UShort_t d, Char_t r);
366 * @param o Object to copy from
368 RingHistos(const RingHistos& o);
370 * Assignment operator
372 * @param o Object to assign from
374 * @return Reference to this
376 RingHistos& operator=(const RingHistos& o);
382 * Initialize the object
386 void SetupForData(const TAxis& eAxis);
390 * @param dir Where to put it
392 void CreateOutputObjects(TList* dir);
394 * Scale the histograms to the total number of events
396 * @param dir Where the output is
397 * @param nEvents Number of events
399 void Terminate(TList* dir, Int_t nEvents);
401 // TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
402 // TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
403 // TProfile* fEtaVsN; // Average uncorrected Nch vs eta
404 // TProfile* fEtaVsM; // Average corrected Nch vs eta
405 TProfile* fCorr; // Average correction vs eta
406 TH2D* fDensity; // Distribution inclusive Nch
407 TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
408 TH1D* fDiffELossPoisson;// Relative difference to Poisson
409 AliPoissonCalculator fPoisson; // Calculate density using Poisson method
410 TH1D* fELoss; // Energy loss as seen by this
411 TH1D* fELossUsed; // Energy loss in strips with signal
412 Double_t fMultCut; // If set, use this
413 TH1D* fTotal; // Total number of strips per eta
414 TH1D* fGood; // Number of good strips per eta
415 TH2D* fPhiAcc; // Phi acceptance vs IpZ
416 TH1D* fPhiBefore; // Phi before re-calce
417 TH1D* fPhiAfter; // Phi after re-calc
418 ClassDef(RingHistos,9);
421 * Get the ring histogram container
426 * @return Ring histogram container
428 RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
429 TList fRingHistos; // List of histogram containers
430 TH1D* fSumOfWeights; // Histogram
431 TH1D* fWeightedSum; // Histogram
432 TH1D* fCorrections; // Histogram
433 UShort_t fMaxParticles; // Maximum particle weight to use
434 Bool_t fUsePoisson; // If true, then use poisson statistics
435 UShort_t fUsePhiAcceptance; // Whether to correct for corners
436 TH1D* fAccI; // Acceptance correction for inner rings
437 TH1D* fAccO; // Acceptance correction for outer rings
438 TArrayI fFMD1iMax; // Array of max weights
439 TArrayI fFMD2iMax; // Array of max weights
440 TArrayI fFMD2oMax; // Array of max weights
441 TArrayI fFMD3iMax; // Array of max weights
442 TArrayI fFMD3oMax; // Array of max weights
443 TH2D* fMaxWeights; // Histogram of max weights
444 TH2D* fLowCuts; // Histogram of low cuts
445 Int_t fEtaLumping; // How to lump eta bins for Poisson
446 Int_t fPhiLumping; // How to lump phi bins for Poisson
447 Int_t fDebug; // Debug level
448 AliFMDMultCuts fCuts; // Cuts
449 Bool_t fRecalculateEta; // Whether to recalc eta and angle correction (disp vtx)
450 Bool_t fRecalculatePhi; // Whether to correct for (X,Y) offset
451 UShort_t fMinQuality; // Least quality for fits
453 ClassDef(AliFMDDensityCalculator,12); // Calculate Nch density