]>
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 | * | |
bd6f5206 | 14 | * @ingroup pwglf_forward_aod |
ffca499d | 15 | */ |
7e4038b5 | 16 | #include <TNamed.h> |
17 | #include <TList.h> | |
1174780f | 18 | #include <TArrayI.h> |
5ca83fee | 19 | #include <TVector3.h> |
7e4038b5 | 20 | #include "AliForwardUtil.h" |
d2638bb7 | 21 | #include "AliFMDMultCuts.h" |
821ffd28 | 22 | #include "AliPoissonCalculator.h" |
7e4038b5 | 23 | class AliESDFMD; |
24 | class TH2D; | |
dd497217 | 25 | class TH1D; |
0bd4b00f | 26 | class TProfile; |
1174780f | 27 | class AliFMDCorrELossFit; |
7e4038b5 | 28 | |
29 | /** | |
30 | * This class calculates the inclusive charged particle density | |
31 | * in each for the 5 FMD rings. | |
32 | * | |
33 | * @par Input: | |
34 | * - AliESDFMD object possibly corrected for sharing | |
35 | * | |
36 | * @par Output: | |
37 | * - 5 RingHistos objects - each with a number of vertex dependent | |
38 | * 2D histograms of the inclusive charge particle density | |
39 | * | |
40 | * @par Corrections used: | |
41 | * - AliFMDAnaCalibEnergyDistribution | |
42 | * - AliFMDDoubleHitCorrection | |
43 | * - AliFMDDeadCorrection | |
44 | * | |
bd6f5206 | 45 | * @ingroup pwglf_forward_algo |
46 | * @ingroup pwglf_forward_aod | |
7e4038b5 | 47 | */ |
48 | class AliFMDDensityCalculator : public TNamed | |
49 | { | |
50 | public: | |
56236b95 | 51 | /** |
52 | * How to correct for the missing phi coverage at the corners of the | |
53 | * sensors | |
54 | * | |
55 | */ | |
56 | enum { | |
57 | /** No correction */ | |
58 | kPhiNoCorrect, | |
59 | /** Correct the calculated number charged particles */ | |
60 | kPhiCorrectNch, | |
61 | /** Correct the energy loss */ | |
62 | kPhiCorrectELoss | |
63 | }; | |
5ca83fee | 64 | /** |
65 | * Folder name | |
66 | */ | |
67 | static const char* fgkFolderName; | |
7e4038b5 | 68 | /** |
69 | * Constructor | |
70 | */ | |
71 | AliFMDDensityCalculator(); | |
72 | /** | |
73 | * Constructor | |
74 | * | |
75 | * @param name Name of object | |
76 | */ | |
77 | AliFMDDensityCalculator(const char* name); | |
78 | /** | |
79 | * Copy constructor | |
80 | * | |
81 | * @param o Object to copy from | |
82 | */ | |
83 | AliFMDDensityCalculator(const AliFMDDensityCalculator& o); | |
84 | /** | |
85 | * Destructor | |
86 | */ | |
87 | virtual ~AliFMDDensityCalculator(); | |
88 | /** | |
89 | * Assignement operator | |
90 | * | |
91 | * @param o Object to assign from | |
92 | * | |
93 | * @return Reference to this object | |
94 | */ | |
7c1a1f1d | 95 | AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o); |
1174780f | 96 | /** |
97 | * Initialize this sub-algorithm | |
98 | * | |
99 | * @param etaAxis Not used | |
100 | */ | |
5934a3e3 | 101 | virtual void SetupForData(const TAxis& etaAxis); |
7e4038b5 | 102 | /** |
103 | * Do the calculations | |
104 | * | |
105 | * @param fmd AliESDFMD object (possibly) corrected for sharing | |
106 | * @param hists Histogram cache | |
7e4038b5 | 107 | * @param lowFlux Low flux flag. |
290052e7 | 108 | * @param cent Centrality |
33438b4c | 109 | * @param ip Coordinates of interaction point |
7e4038b5 | 110 | * |
111 | * @return true on successs | |
112 | */ | |
290052e7 | 113 | virtual Bool_t Calculate(const AliESDFMD& fmd, |
7e4038b5 | 114 | AliForwardUtil::Histos& hists, |
290052e7 | 115 | Bool_t lowFlux, |
116 | Double_t cent=-1, | |
5ca83fee | 117 | const TVector3& ip=TVector3(1024,1024,0)); |
7e4038b5 | 118 | /** |
119 | * Scale the histograms to the total number of events | |
120 | * | |
c389303e | 121 | * @param dir where to put the output |
33438b4c | 122 | * @param output Output list |
7e4038b5 | 123 | * @param nEvents Number of events |
124 | */ | |
5934a3e3 | 125 | virtual void Terminate(const TList* dir, TList* output, Int_t nEvents); |
7e4038b5 | 126 | /** |
127 | * Output diagnostic histograms to directory | |
128 | * | |
129 | * @param dir List to write in | |
130 | */ | |
5934a3e3 | 131 | virtual void CreateOutputObjects(TList* dir); |
ea3e5d95 | 132 | /** |
133 | * Set the debug level. The higher the value the more output | |
134 | * | |
135 | * @param dbg Debug level | |
136 | */ | |
1f7aa5c7 | 137 | void SetDebug(Int_t dbg=1) { fDebug = dbg; } |
138 | void SetDoTiming(Bool_t enable=true) { fDoTiming = enable; } | |
0bd4b00f | 139 | /** |
140 | * Maximum particle weight to use | |
141 | * | |
142 | * @param m | |
143 | */ | |
144 | void SetMaxParticles(UShort_t m) { fMaxParticles = m; } | |
9d05ffeb | 145 | /** |
146 | * Set whether to use poisson statistics to estimate the | |
147 | * number of particles that has hit within a region. If this is true, | |
148 | * then the average charge particle density is given by | |
149 | * @f[ | |
ffca499d | 150 | * \lambda = -\log\left(\frac{N_e}{N_t}\right) |
9d05ffeb | 151 | * @f] |
152 | * where $N_e$ is the number of strips within the region that has no | |
153 | * hits over threshold, and $N_t$ is the total number of strips in the | |
154 | * region/ | |
155 | * | |
156 | * @param u Whether to use poisson statistics to estimate the | |
157 | * number of particles that has hit within a region. | |
158 | */ | |
159 | void SetUsePoisson(Bool_t u) { fUsePoisson = u; } | |
6f4a5c0d | 160 | /** |
161 | * In case of a displaced vertices recalculate eta and angle correction | |
162 | * | |
163 | * @param use recalculate or not | |
164 | * | |
165 | */ | |
166 | void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; } | |
5ca83fee | 167 | /** |
168 | * In case of a displaced vertices recalculate eta and angle correction | |
169 | * | |
170 | * @param use recalculate or not | |
171 | * | |
172 | */ | |
173 | void SetRecalculatePhi(Bool_t use) { fRecalculatePhi = use; } | |
5bb5d1f6 | 174 | /** |
175 | * Set whether to use the phi acceptance correction. | |
176 | * | |
0082a8fc | 177 | * How the phi acceptance is used depends on the value passed. |
178 | * - 0: No phi acceptance | |
179 | * - 1: Phi acceptance correction done to estimate of particles | |
180 | * - 2: Phi acceptance correction done to energy deposited | |
181 | * | |
182 | * @param u If >0, use the phi acceptance (default is false) | |
5bb5d1f6 | 183 | */ |
6f4a5c0d | 184 | |
56236b95 | 185 | void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; } |
b6b35c77 | 186 | /** |
187 | * Set the luming factors used in the Poisson method | |
188 | * | |
189 | * @param eta Must be 1 or larger | |
190 | * @param phi Must be 1 or larger | |
191 | */ | |
192 | void SetLumping(Int_t eta, Int_t phi) { | |
193 | fEtaLumping = (eta < 1 ? 1 : eta); | |
194 | fPhiLumping = (phi < 1 ? 1 : phi); | |
195 | } | |
8449e3e0 | 196 | /** |
197 | * Set the minimum quality of the energy loss fits | |
198 | * | |
199 | * @param cut Cut value | |
200 | */ | |
81775aba | 201 | void SetMinQuality(UShort_t cut=10) { fMinQuality = cut; } |
77f97e3f CHC |
202 | /** |
203 | * Set the maximum ratio of outlier bins to the total number of bins | |
204 | * with data. | |
205 | * | |
206 | * @param ratio Maximum ratio (number between 0 and 1) | |
207 | */ | |
208 | void SetMaxOutliers(Double_t ratio=0.10) { fMaxOutliers = ratio; } | |
209 | /** | |
210 | * Set the maximum relative diviation between @f$N_{ch}^{Poisson}@f$ | |
211 | * and @f$N_{ch}^{\Delta}@f$ | |
212 | * | |
213 | * @param cut Relative cut (number between 0 and 1) | |
214 | */ | |
215 | void SetOutlierCut(Double_t cut=0.50) { fOutlierCut = cut; } | |
5bb5d1f6 | 216 | /** |
217 | * Get the multiplicity cut. If the user has set fMultCut (via | |
218 | * SetMultCut) then that value is used. If not, then the lower | |
219 | * value of the fit range for the enery loss fits is returned. | |
220 | * | |
290052e7 | 221 | * @param d Detector |
222 | * @param r Ring | |
223 | * @param eta Psuedo-rapidity | |
224 | * @param errors Factor in errors | |
225 | * | |
5bb5d1f6 | 226 | * @return Lower cut on multiplicity |
227 | */ | |
228 | Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta, | |
229 | Bool_t errors=true) const; | |
0bd4b00f | 230 | /** |
231 | * Get the multiplicity cut. If the user has set fMultCut (via | |
232 | * SetMultCut) then that value is used. If not, then the lower | |
233 | * value of the fit range for the enery loss fits is returned. | |
234 | * | |
290052e7 | 235 | * @param d Detector |
236 | * @param r Ring | |
237 | * @param ieta Psuedo-rapidity bin | |
238 | * @param errors Factor in errors | |
239 | * | |
0bd4b00f | 240 | * @return Lower cut on multiplicity |
241 | */ | |
5bb5d1f6 | 242 | Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta, |
243 | Bool_t errors=true) const; | |
8449e3e0 | 244 | /** |
245 | * Set the minimum quality of the energy loss fits | |
246 | * | |
c8b1a7db | 247 | * @return Cut value |
8449e3e0 | 248 | */ |
249 | UShort_t GetMinQuality() const { return fMinQuality; } | |
0bd4b00f | 250 | /** |
251 | * Print information | |
252 | * | |
1174780f | 253 | * @param option Print options |
254 | * - max Print max weights | |
0bd4b00f | 255 | */ |
256 | void Print(Option_t* option="") const; | |
e18cb8bd | 257 | /** |
258 | * Get the cuts used | |
259 | * | |
260 | * @return Reference to cuts object | |
261 | */ | |
d2638bb7 | 262 | AliFMDMultCuts& GetCuts() { return fCuts; } |
e18cb8bd | 263 | /** |
264 | * Set the cuts to use | |
265 | * | |
266 | * @param c Cuts to use | |
267 | */ | |
d2638bb7 | 268 | void SetCuts(const AliFMDMultCuts& c) { fCuts = c; } |
7e4038b5 | 269 | protected: |
1174780f | 270 | /** |
271 | * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta | |
272 | * | |
273 | * @param cor Correction | |
274 | * @param d Detector | |
275 | * @param r Ring | |
276 | * @param iEta Eta bin | |
290052e7 | 277 | * |
278 | * @return The maximum weight | |
1174780f | 279 | */ |
fb3430ac | 280 | Int_t FindMaxWeight(const AliFMDCorrELossFit* cor, |
1174780f | 281 | UShort_t d, Char_t r, Int_t iEta) const; |
282 | ||
7095962e CHC |
283 | /** |
284 | * Find the max weight to use for FMD<i>dr</i> in eta @a eta | |
285 | * | |
286 | * @param cor Correction | |
287 | * @param d Detector | |
288 | * @param r Ring | |
289 | * @param eta Eta | |
290 | * | |
291 | * @return The maximum weight | |
292 | */ | |
293 | Int_t FindMaxWeight(const AliFMDCorrELossFit* cor, | |
294 | UShort_t d, Char_t r, Double_t iEta) const; | |
295 | ||
1174780f | 296 | /** |
297 | * Find the max weights and cache them | |
298 | * | |
e2ebf8c4 | 299 | * @param axis Default @f$\eta@f$ axis from parent task |
1174780f | 300 | */ |
e2ebf8c4 | 301 | void CacheMaxWeights(const TAxis& axis); |
1174780f | 302 | /** |
303 | * Find the (cached) maximum weight for FMD<i>dr</i> in | |
304 | * @f$\eta@f$ bin @a iEta | |
305 | * | |
306 | * @param d Detector | |
307 | * @param r Ring | |
308 | * @param iEta Eta bin | |
309 | * | |
310 | * @return max weight or <= 0 in case of problems | |
311 | */ | |
312 | Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const; | |
313 | /** | |
314 | * Find the (cached) maximum weight for FMD<i>dr</i> iat | |
315 | * @f$\eta@f$ | |
316 | * | |
317 | * @param d Detector | |
318 | * @param r Ring | |
319 | * @param eta Eta bin | |
320 | * | |
321 | * @return max weight or <= 0 in case of problems | |
322 | */ | |
323 | Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const; | |
324 | ||
7e4038b5 | 325 | /** |
326 | * Get the number of particles corresponding to the signal mult | |
327 | * | |
328 | * @param mult Signal | |
329 | * @param d Detector | |
330 | * @param r Ring | |
7e4038b5 | 331 | * @param eta Pseudo-rapidity |
332 | * @param lowFlux Low-flux flag | |
333 | * | |
334 | * @return The number of particles | |
335 | */ | |
5ca83fee | 336 | virtual Float_t NParticles(Float_t mult, |
337 | UShort_t d, | |
338 | Char_t r, | |
339 | Float_t eta, | |
340 | Bool_t lowFlux) const; | |
7e4038b5 | 341 | /** |
342 | * Get the inverse correction factor. This consist of | |
343 | * | |
344 | * - acceptance correction (corners of sensors) | |
345 | * - double hit correction (for low-flux events) | |
346 | * - dead strip correction | |
347 | * | |
348 | * @param d Detector | |
349 | * @param r Ring | |
33438b4c | 350 | * @param t Strip |
7e4038b5 | 351 | * @param eta Pseudo-rapidity |
352 | * @param lowFlux Low-flux flag | |
353 | * | |
33438b4c | 354 | * @return the correction factor |
7e4038b5 | 355 | */ |
5ca83fee | 356 | virtual Float_t Correction(UShort_t d, Char_t r, UShort_t t, |
357 | Float_t eta, Bool_t lowFlux) const; | |
7e4038b5 | 358 | /** |
359 | * Get the acceptance correction for strip @a t in an ring of type @a r | |
360 | * | |
361 | * @param r Ring type ('I' or 'O') | |
362 | * @param t Strip number | |
363 | * | |
364 | * @return Inverse acceptance correction | |
365 | */ | |
366 | virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const; | |
0bd4b00f | 367 | /** |
368 | * Generate the acceptance corrections | |
369 | * | |
370 | * @param r Ring to generate for | |
371 | * | |
372 | * @return Newly allocated histogram of acceptance corrections | |
373 | */ | |
374 | virtual TH1D* GenerateAcceptanceCorrection(Char_t r) const; | |
77f97e3f CHC |
375 | /** |
376 | * Check if, for a given region, whether this is an outlier | |
377 | * | |
378 | * The condition for an outlier event are | |
379 | * @f[ | |
380 | * |N_{ch}^{Poisson} - N_{ch}^{\Delta}| / N_{ch}^{\Delta} > c | |
381 | * @f] | |
382 | * | |
383 | * @param eloss @f$ N_{ch}^{\Delta}@f$ - number of charged particles | |
384 | * @param poisson @f$ N_{ch}^{Poisson}@f$ - number of charged particles | |
385 | * @param cut @f$ c@f$ - the cut | |
386 | * | |
387 | * @return true if the region reflects an outlier event | |
388 | */ | |
389 | virtual Bool_t CheckOutlier(Double_t eloss, | |
390 | Double_t poisson, | |
391 | Double_t cut=0.5) const; | |
7e4038b5 | 392 | /** |
393 | * Internal data structure to keep track of the histograms | |
394 | */ | |
9d99b0dd | 395 | struct RingHistos : public AliForwardUtil::RingHistos |
7e4038b5 | 396 | { |
397 | /** | |
398 | * Default CTOR | |
399 | */ | |
400 | RingHistos(); | |
401 | /** | |
402 | * Constructor | |
403 | * | |
404 | * @param d detector | |
405 | * @param r ring | |
406 | */ | |
407 | RingHistos(UShort_t d, Char_t r); | |
408 | /** | |
409 | * Copy constructor | |
410 | * | |
411 | * @param o Object to copy from | |
412 | */ | |
413 | RingHistos(const RingHistos& o); | |
414 | /** | |
415 | * Assignment operator | |
416 | * | |
417 | * @param o Object to assign from | |
418 | * | |
419 | * @return Reference to this | |
420 | */ | |
421 | RingHistos& operator=(const RingHistos& o); | |
422 | /** | |
423 | * Destructor | |
424 | */ | |
425 | ~RingHistos(); | |
9d05ffeb | 426 | /** |
427 | * Initialize the object | |
428 | * | |
429 | * @param eAxis | |
430 | */ | |
5934a3e3 | 431 | void SetupForData(const TAxis& eAxis); |
c389303e | 432 | /** |
433 | * Make output | |
434 | * | |
435 | * @param dir Where to put it | |
436 | */ | |
5934a3e3 | 437 | void CreateOutputObjects(TList* dir); |
9d99b0dd | 438 | /** |
439 | * Scale the histograms to the total number of events | |
440 | * | |
c389303e | 441 | * @param dir Where the output is |
9d99b0dd | 442 | * @param nEvents Number of events |
443 | */ | |
5934a3e3 | 444 | void Terminate(TList* dir, Int_t nEvents); |
5ca83fee | 445 | TList* fList; |
8449e3e0 | 446 | // TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch |
447 | // TH2D* fEvsM; // Correlation of Eloss vs corrected Nch | |
448 | // TProfile* fEtaVsN; // Average uncorrected Nch vs eta | |
449 | // TProfile* fEtaVsM; // Average corrected Nch vs eta | |
b6b35c77 | 450 | TProfile* fCorr; // Average correction vs eta |
451 | TH2D* fDensity; // Distribution inclusive Nch | |
9d05ffeb | 452 | TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch |
5ca83fee | 453 | TH1D* fDiffELossPoisson;// Relative difference to Poisson |
77f97e3f CHC |
454 | TH2D* fELossVsPoissonOut; // Correlation of energy loss vs Poisson N_ch |
455 | TH1D* fDiffELossPoissonOut;// Relative difference to Poisson | |
456 | TH1D* fOutliers; // Fraction of outliers per event | |
821ffd28 | 457 | AliPoissonCalculator fPoisson; // Calculate density using Poisson method |
f7cfc454 | 458 | TH1D* fELoss; // Energy loss as seen by this |
459 | TH1D* fELossUsed; // Energy loss in strips with signal | |
460 | Double_t fMultCut; // If set, use this | |
5ca83fee | 461 | TH1D* fTotal; // Total number of strips per eta |
462 | TH1D* fGood; // Number of good strips per eta | |
463 | TH2D* fPhiAcc; // Phi acceptance vs IpZ | |
464 | TH1D* fPhiBefore; // Phi before re-calce | |
465 | TH1D* fPhiAfter; // Phi after re-calc | |
77f97e3f | 466 | ClassDef(RingHistos,10); |
7e4038b5 | 467 | }; |
468 | /** | |
469 | * Get the ring histogram container | |
470 | * | |
471 | * @param d Detector | |
472 | * @param r Ring | |
473 | * | |
474 | * @return Ring histogram container | |
475 | */ | |
476 | RingHistos* GetRingHistos(UShort_t d, Char_t r) const; | |
ea3e5d95 | 477 | TList fRingHistos; // List of histogram containers |
1174780f | 478 | TH1D* fSumOfWeights; // Histogram |
479 | TH1D* fWeightedSum; // Histogram | |
480 | TH1D* fCorrections; // Histogram | |
0bd4b00f | 481 | UShort_t fMaxParticles; // Maximum particle weight to use |
9d05ffeb | 482 | Bool_t fUsePoisson; // If true, then use poisson statistics |
0082a8fc | 483 | UShort_t fUsePhiAcceptance; // Whether to correct for corners |
0bd4b00f | 484 | TH1D* fAccI; // Acceptance correction for inner rings |
485 | TH1D* fAccO; // Acceptance correction for outer rings | |
1174780f | 486 | TArrayI fFMD1iMax; // Array of max weights |
487 | TArrayI fFMD2iMax; // Array of max weights | |
488 | TArrayI fFMD2oMax; // Array of max weights | |
489 | TArrayI fFMD3iMax; // Array of max weights | |
490 | TArrayI fFMD3oMax; // Array of max weights | |
5bb5d1f6 | 491 | TH2D* fMaxWeights; // Histogram of max weights |
492 | TH2D* fLowCuts; // Histogram of low cuts | |
b6b35c77 | 493 | Int_t fEtaLumping; // How to lump eta bins for Poisson |
494 | Int_t fPhiLumping; // How to lump phi bins for Poisson | |
ea3e5d95 | 495 | Int_t fDebug; // Debug level |
e18cb8bd | 496 | AliFMDMultCuts fCuts; // Cuts |
8449e3e0 | 497 | Bool_t fRecalculateEta; // Whether to recalc eta and angle correction (disp vtx) |
498 | Bool_t fRecalculatePhi; // Whether to correct for (X,Y) offset | |
499 | UShort_t fMinQuality; // Least quality for fits | |
1f7aa5c7 | 500 | AliForwardUtil::Histos fCache; |
501 | Bool_t fDoTiming; | |
502 | TProfile* fHTiming; | |
77f97e3f CHC |
503 | Double_t fMaxOutliers; // Maximum ratio of outlier bins |
504 | Double_t fOutlierCut; // Maximum relative diviation | |
7e4038b5 | 505 | |
77f97e3f | 506 | ClassDef(AliFMDDensityCalculator,14); // Calculate Nch density |
7e4038b5 | 507 | }; |
508 | ||
509 | #endif | |
510 | // Local Variables: | |
511 | // mode: C++ | |
512 | // End: | |
513 |