2 // Utilities used in the forward multiplcity analysis
5 #ifndef ALIFORWARDUTIL_H
6 #define ALIFORWARDUTIL_H
8 * @file AliForwardUtil.h
9 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
10 * @date Wed Mar 23 14:06:54 2011
15 * @ingroup pwglf_forward
19 #include <TObjArray.h>
28 * Utilities used in the forward multiplcity analysis
30 * @ingroup pwglf_forward
32 class AliForwardUtil : public TObject
36 * Get the standard color for a ring
43 static Color_t RingColor(UShort_t d, Char_t r)
45 return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
46 + ((r == 'I' || r == 'i') ? 2 : -3));
48 //==================================================================
51 * @name Collision/run parameters
54 * Defined collision types
56 enum ECollisionSystem {
62 //__________________________________________________________________
64 * Parse a collision system spec given in a string. Known values are
66 * - "pp", "p-p" which returns kPP
67 * - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
68 * - "pPb", "p-Pb", "pA", p-A" which returns kPPb
69 * - Everything else gives kUnknown
71 * @param sys Collision system spec
73 * @return Collision system id
75 static UShort_t ParseCollisionSystem(const char* sys);
77 * Get a string representation of the collision system
79 * @param sys Collision system
83 * - anything else gives "unknown"
85 * @return String representation of the collision system
87 static const char* CollisionSystemString(UShort_t sys);
88 //__________________________________________________________________
90 * Parse the center of mass energy given as a float and return known
91 * values as a unsigned integer
93 * @param sys Collision system (needed for AA)
94 * @param cms Center of mass energy * total charge
96 * @return Center of mass energy per nucleon
98 static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
100 * Get a string representation of the center of mass energy per nuclean
102 * @param cms Center of mass energy per nucleon
104 * @return String representation of the center of mass energy per nuclean
106 static const char* CenterOfMassEnergyString(UShort_t cms);
107 //__________________________________________________________________
109 * Parse the magnetic field (in kG) as given by a floating point number
111 * @param field Magnetic field in kG
113 * @return Short integer value of magnetic field in kG
115 static Short_t ParseMagneticField(Float_t field);
119 * @param det, ring, sec, strip, zvtx
123 static Double_t GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx) ;
125 * Get a string representation of the magnetic field
127 * @param field Magnetic field in kG
129 * @return String representation of the magnetic field
131 static const char* MagneticFieldString(Short_t field);
136 * @name Energy stragling functions
138 //__________________________________________________________________
140 * Number of steps to do in the Landau, Gaussiam convolution
142 static Int_t fgConvolutionSteps; // Number of convolution steps
143 //------------------------------------------------------------------
145 * How many sigma's of the Gaussian in the Landau, Gaussian
146 * convolution to integrate over
148 static Double_t fgConvolutionNSigma; // Number of convolution sigmas
149 //------------------------------------------------------------------
151 * Calculate the shifted Landau
153 * f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
156 * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
157 * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
158 * @f$\Delta=0,\xi=1@f$.
160 * @param x Where to evaluate @f$ f'_{L}@f$
161 * @param delta Most probable value
162 * @param xi The 'width' of the distribution
164 * @return @f$ f'_{L}(x;\Delta,\xi) @f$
166 static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
168 //------------------------------------------------------------------
170 * Calculate the value of a Landau convolved with a Gaussian
173 * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
174 * \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
175 * \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
178 * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
179 * energy loss, @f$ \xi@f$ the width of the Landau, and
180 * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$. Here, @f$\sigma@f$ is the
181 * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
182 * noise in the detector.
184 * Note that this function uses the constants fgConvolutionSteps and
185 * fgConvolutionNSigma
188 * - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
189 * - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
190 * - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
192 * @param x where to evaluate @f$ f@f$
193 * @param delta @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
194 * @param xi @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
195 * @param sigma @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
196 * @param sigma_n @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
198 * @return @f$ f@f$ evaluated at @f$ x@f$.
200 static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi,
201 Double_t sigma, Double_t sigma_n);
203 //------------------------------------------------------------------
207 * f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
209 * corresponding to @f$ i@f$ particles i.e., with the substitutions
211 * \Delta \rightarrow \Delta_i &=& i(\Delta + \xi\log(i))\\
212 * \xi \rightarrow \xi_i &=& i \xi\\
213 * \sigma \rightarrow \sigma_i &=& \sqrt{i}\sigma\\
214 * \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
217 * @param x Where to evaluate
218 * @param delta @f$ \Delta@f$
219 * @param xi @f$ \xi@f$
220 * @param sigma @f$ \sigma@f$
221 * @param sigma_n @f$ \sigma_n@f$
224 * @return @f$ f_i @f$ evaluated
226 static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi,
227 Double_t sigma, Double_t sigma_n, Int_t i);
229 //------------------------------------------------------------------
231 * Numerically evaluate
233 * \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
235 * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter. The mapping
236 * of the parameters is given by
241 * - 3: @f$\sigma_n@f$
243 * This is the partial derivative with respect to the parameter of
244 * the response function corresponding to @f$ i@f$ particles i.e.,
245 * with the substitutions
247 * \Delta \rightarrow \Delta_i = i(\Delta + \xi\log(i))\\
248 * \xi \rightarrow \xi_i = i \xi\\
249 * \sigma \rightarrow \sigma_i = \sqrt{i}\sigma\\
250 * \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
253 * @param x Where to evaluate
254 * @param ipar Parameter number
255 * @param dp @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
256 * @param delta @f$ \Delta@f$
257 * @param xi @f$ \xi@f$
258 * @param sigma @f$ \sigma@f$
259 * @param sigma_n @f$ \sigma_n@f$
262 * @return @f$ f_i@f$ evaluated
264 static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
265 Double_t delta, Double_t xi,
266 Double_t sigma, Double_t sigma_n, Int_t i);
268 //------------------------------------------------------------------
272 * f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
275 * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
276 * Landau with a Gaussian (see LandauGaus). Note that
277 * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
278 * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
281 * - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
282 * - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
283 * - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
285 * @param x Where to evaluate @f$ f_N@f$
286 * @param delta @f$ \Delta_1@f$
287 * @param xi @f$ \xi_1@f$
288 * @param sigma @f$ \sigma_1@f$
289 * @param sigma_n @f$ \sigma_n@f$
290 * @param n @f$ N@f$ in the sum above.
291 * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
294 * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
296 static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
297 Double_t sigma, Double_t sigma_n, Int_t n,
300 * Generate a TF1 object of @f$ f_I@f$
303 * @param delta @f$ \Delta@f$
304 * @param xi @f$ \xi_1@f$
305 * @param sigma @f$ \sigma_1@f$
306 * @param sigma_n @f$ \sigma_n@f$
307 * @param i @f$ i@f$ - the number of particles
308 * @param xmin Least value of range
309 * @param xmax Largest value of range
311 * @return Newly allocated TF1 object
313 static TF1* MakeILandauGaus(Double_t c,
314 Double_t delta, Double_t xi,
315 Double_t sigma, Double_t sigma_n,
317 Double_t xmin, Double_t xmax);
319 * Generate a TF1 object of @f$ f_N@f$
322 * @param delta @f$ \Delta@f$
323 * @param xi @f$ \xi_1@f$
324 * @param sigma @f$ \sigma_1@f$
325 * @param sigma_n @f$ \sigma_n@f$
326 * @param n @f$ N@f$ - how many particles to sum to
327 * @param a Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
329 * @param xmin Least value of range
330 * @param xmax Largest value of range
332 * @return Newly allocated TF1 object
334 static TF1* MakeNLandauGaus(Double_t c,
335 Double_t delta, Double_t xi,
336 Double_t sigma, Double_t sigma_n,
337 Int_t n, const Double_t* a,
338 Double_t xmin, Double_t xmax);
340 //__________________________________________________________________
342 * Structure to do fits to the energy loss spectrum
344 * @ingroup pwglf_forward
360 * @param lowCut Lower cut of spectrum - data below this cuts is ignored
361 * @param maxRange Maximum range to fit to
362 * @param minusBins The number of bins below maximum to use
364 ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins);
369 virtual ~ELossFitter();
371 * Clear internal arrays
376 * Fit a 1-particle signal to the passed energy loss distribution
378 * Note that this function clears the internal arrays first
380 * @param dist Data to fit the function to
381 * @param sigman If larger than zero, the initial guess of the
382 * detector induced noise. If zero or less, then this
383 * parameter is ignored in the fit (fixed at 0)
385 * @return The function fitted to the data
387 TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
389 * Fit a N-particle signal to the passed energy loss distribution
391 * If there's no 1-particle fit present, it does that first
393 * @param dist Data to fit the function to
394 * @param n Number of particle signals to fit
395 * @param sigman If larger than zero, the initial guess of the
396 * detector induced noise. If zero or less, then this
397 * parameter is ignored in the fit (fixed at 0)
399 * @return The function fitted to the data
401 TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
403 * Get Lower cut on data
405 * @return Lower cut on data
407 Double_t GetLowCut() const { return fLowCut; }
409 * Get Maximum range to fit
411 * @return Maximum range to fit
413 Double_t GetMaxRange() const { return fMaxRange; }
415 * Get Number of bins from maximum to fit 1st peak
417 * @return Number of bins from maximum to fit 1st peak
419 UShort_t GetMinusBins() const { return fMinusBins; }
421 * Get Array of fit results
423 * @return Array of fit results
425 const TObjArray& GetFitResults() const { return fFitResults; }
427 * Get Array of fit results
429 * @return Array of fit results
431 TObjArray& GetFitResults() { return fFitResults; }
433 * Get Array of functions
435 * @return Array of functions
437 const TObjArray& GetFunctions() const { return fFunctions; }
439 * Get Array of functions
441 * @return Array of functions
443 TObjArray& GetFunctions() { return fFunctions; }
445 const Double_t fLowCut; // Lower cut on data
446 const Double_t fMaxRange; // Maximum range to fit
447 const UShort_t fMinusBins; // Number of bins from maximum to fit 1st peak
448 TObjArray fFitResults; // Array of fit results
449 TObjArray fFunctions; // Array of functions
454 //==================================================================
457 * @name Convenience containers
460 * Structure to hold histograms
462 * @ingroup pwglf_forward
464 struct Histos : public TObject
471 Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
475 * @param o Object to copy from
477 Histos(const Histos& o)
486 * Assignement operator
488 * @return Reference to this
490 Histos& operator=(const Histos&) { return *this;}
496 * Initialize the object
498 * @param etaAxis Eta axis to use
500 void Init(const TAxis& etaAxis);
506 * @param etaAxis Eta axis to use
508 * @return Newly allocated histogram
510 TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
514 * @param option Not used
516 void Clear(Option_t* option="");
517 // const TH2D* Get(UShort_t d, Char_t r) const;
519 * Get the histogram for a particular detector,ring
524 * @return Histogram for detector,ring or nul
526 TH2D* Get(UShort_t d, Char_t r) const;
527 TH2D* fFMD1i; // Histogram for FMD1i
528 TH2D* fFMD2i; // Histogram for FMD2i
529 TH2D* fFMD2o; // Histogram for FMD2o
530 TH2D* fFMD3i; // Histogram for FMD3i
531 TH2D* fFMD3o; // Histogram for FMD3o
536 //__________________________________________________________________
538 * Base class for structure holding ring specific histograms
540 * @ingroup pwglf_forward
542 struct RingHistos : public TObject
548 RingHistos() : fDet(0), fRing('\0'), fName("") {}
555 RingHistos(UShort_t d, Char_t r)
556 : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r))
561 * @param o Object to copy from
563 RingHistos(const RingHistos& o)
564 : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
569 virtual ~RingHistos() {}
571 * Assignement operator
573 * @param o Object to assign from
575 * @return Reference to this
577 RingHistos& operator=(const RingHistos& o)
579 if (&o == this) return *this;
580 TObject::operator=(o);
587 * Define the outout list in @a d
589 * @param d Where to put the output list
591 * @return Newly allocated TList object or null
593 TList* DefineOutputList(TList* d) const;
595 * Get our output list from the container @a d
597 * @param d where to get the output list from
599 * @return The found TList or null
601 TList* GetOutputList(const TList* d) const;
603 * Find a specific histogram in the source list @a d
605 * @param d (top)-container
606 * @param name Name of histogram
608 * @return Found histogram or null
610 TH1* GetOutputHist(const TList* d, const char* name) const;
617 Color_t Color() const
619 return AliForwardUtil::RingColor(fDet, fRing);
621 const char* GetName() const { return fName.Data(); }
622 UShort_t fDet; // Detector
623 Char_t fRing; // Ring
624 TString fName; // Name
626 ClassDef(RingHistos,1)
630 //__________________________________________________________________
633 DebugGuard(Int_t lvl, Int_t msgLvl, const char* format, ...);
636 void Format(bool in);
641 AliForwardUtil(const AliForwardUtil& o) : TObject(o) {}
642 AliForwardUtil& operator=(const AliForwardUtil&) { return *this; }
646 ClassDef(AliForwardUtil,1) // Utilities - do not make object
650 # define DGUARD(L,N,F,...) do {} while(false)
652 # define DGUARD(L,N,F,...) \
653 AliForwardUtil::DebugGuard _GUARD(L,N,F, ## __VA_ARGS__)