f2728ea53b46c1bbd45324d78dfedbf5a367bf87
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.h
1 // 
2 // Utilities used in the forward multiplcity analysis 
3 // 
4 //
5 #ifndef ALIFORWARDUTIL_H
6 #define ALIFORWARDUTIL_H
7 /**
8  * @file   AliForwardUtil.h
9  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
10  * @date   Wed Mar 23 14:06:54 2011
11  * 
12  * @brief  
13  * 
14  * 
15  * @ingroup pwg2_forward 
16  */
17 #include <TObject.h>
18 #include <TString.h>
19 #include <TObjArray.h>
20 class TH2D;
21 class TH1I;
22 class TH1;
23 class TF1;
24 class TAxis;
25 class AliESDEvent;
26
27 /** 
28  * Utilities used in the forward multiplcity analysis 
29  * 
30  * @ingroup pwg2_forward 
31  */
32 class AliForwardUtil : public TObject
33 {
34 public:
35   /** 
36    * Get the standard color for a ring  
37    *
38    * @param d Detector
39    * @param r Ring 
40    * 
41    * @return 
42    */
43   static Color_t RingColor(UShort_t d, Char_t r)
44   { 
45     return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
46             + ((r == 'I' || r == 'i') ? 2 : -3));
47   }
48   //==================================================================
49   /** 
50    * @{ 
51    * @name Collision/run parameters 
52    */
53   /**                                           
54    * Defined collision types 
55    */
56   enum ECollisionSystem {
57     kUnknown, 
58     kPP, 
59     kPbPb,
60     kPPb
61   };
62   //__________________________________________________________________
63   /** 
64    * Parse a collision system spec given in a string.   Known values are 
65    * 
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 
70    * 
71    * @param sys Collision system spec 
72    * 
73    * @return Collision system id 
74    */
75   static UShort_t ParseCollisionSystem(const char* sys);
76   /** 
77    * Get a string representation of the collision system 
78    * 
79    * @param sys  Collision system 
80    * - kPP -> "pp"
81    * - kPbPb -> "PbPb" 
82    * - kPPb -> "pPb"
83    * - anything else gives "unknown"
84    * 
85    * @return String representation of the collision system 
86    */
87   static const char* CollisionSystemString(UShort_t sys);
88   //__________________________________________________________________
89   /** 
90    * Parse the center of mass energy given as a float and return known 
91    * values as a unsigned integer
92    * 
93    * @param sys  Collision system (needed for AA)
94    * @param cms  Center of mass energy * total charge 
95    * 
96    * @return Center of mass energy per nucleon
97    */
98   static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
99   /** 
100    * Get a string representation of the center of mass energy per nuclean
101    * 
102    * @param cms  Center of mass energy per nucleon
103    * 
104    * @return String representation of the center of mass energy per nuclean
105    */
106   static const char* CenterOfMassEnergyString(UShort_t cms);
107   //__________________________________________________________________
108   /** 
109    * Parse the magnetic field (in kG) as given by a floating point number
110    * 
111    * @param field  Magnetic field in kG 
112    * 
113    * @return Short integer value of magnetic field in kG 
114    */
115   static Short_t ParseMagneticField(Float_t field);
116   /** 
117    * Get a string representation of the magnetic field
118    * 
119    * @param field Magnetic field in kG
120    * 
121    * @return String representation of the magnetic field
122    */
123   static const char* MagneticFieldString(Short_t field);
124   /* @} */
125
126   /** 
127    * @{ 
128    * @name Energy stragling functions 
129    */
130   //__________________________________________________________________
131   /**
132    * Number of steps to do in the Landau, Gaussiam convolution 
133    */
134   static Int_t fgConvolutionSteps; // Number of convolution steps
135   //------------------------------------------------------------------
136   /** 
137    * How many sigma's of the Gaussian in the Landau, Gaussian
138    * convolution to integrate over
139    */
140   static Double_t fgConvolutionNSigma; // Number of convolution sigmas 
141   //------------------------------------------------------------------
142   /** 
143    * Calculate the shifted Landau
144    * @f[
145    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
146    * @f]
147    *
148    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
149    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
150    * @f$\Delta=0,\xi=1@f$. 
151    *
152    * @param x      Where to evaluate @f$ f'_{L}@f$ 
153    * @param delta  Most probable value 
154    * @param xi     The 'width' of the distribution 
155    *
156    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
157    */
158   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
159   
160   //------------------------------------------------------------------
161   /** 
162    * Calculate the value of a Landau convolved with a Gaussian 
163    * 
164    * @f[ 
165    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
166    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
167    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
168    * @f]
169    * 
170    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
171    * energy loss, @f$ \xi@f$ the width of the Landau, and 
172    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
173    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
174    * noise in the detector.  
175    *
176    * Note that this function uses the constants fgConvolutionSteps and
177    * fgConvolutionNSigma
178    * 
179    * References: 
180    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
181    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
182    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
183    * 
184    * @param x         where to evaluate @f$ f@f$
185    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
186    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
187    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
188    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
189    * 
190    * @return @f$ f@f$ evaluated at @f$ x@f$.  
191    */
192   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
193                              Double_t sigma, Double_t sigma_n);
194
195   //------------------------------------------------------------------
196   /** 
197    * Evaluate 
198    * @f[ 
199    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
200    * @f] 
201    * corresponding to @f$ i@f$ particles i.e., with the substitutions 
202    * @f{eqnarray*}{ 
203    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
204    *    \xi       \rightarrow \xi_i       &=& i \xi\\
205    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
206    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
207    * @f} 
208    * 
209    * @param x        Where to evaluate 
210    * @param delta    @f$ \Delta@f$ 
211    * @param xi       @f$ \xi@f$ 
212    * @param sigma    @f$ \sigma@f$ 
213    * @param sigma_n  @f$ \sigma_n@f$
214    * @param i        @f$ i @f$
215    * 
216    * @return @f$ f_i @f$ evaluated
217    */  
218   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
219                               Double_t sigma, Double_t sigma_n, Int_t i);
220
221   //------------------------------------------------------------------
222   /** 
223    * Numerically evaluate 
224    * @f[ 
225    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
226    * @f] 
227    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
228    * of the parameters is given by 
229    *
230    * - 0: @f$\Delta@f$ 
231    * - 1: @f$\xi@f$ 
232    * - 2: @f$\sigma@f$ 
233    * - 3: @f$\sigma_n@f$ 
234    *
235    * This is the partial derivative with respect to the parameter of
236    * the response function corresponding to @f$ i@f$ particles i.e.,
237    * with the substitutions
238    * @f[ 
239    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
240    *    \xi       \rightarrow \xi_i       = i \xi\\
241    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
242    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
243    * @f] 
244    * 
245    * @param x        Where to evaluate 
246    * @param ipar     Parameter number 
247    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
248    * @param delta    @f$ \Delta@f$ 
249    * @param xi       @f$ \xi@f$ 
250    * @param sigma    @f$ \sigma@f$ 
251    * @param sigma_n  @f$ \sigma_n@f$
252    * @param i        @f$ i@f$
253    * 
254    * @return @f$ f_i@f$ evaluated
255    */  
256   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
257                                    Double_t delta, Double_t xi, 
258                                    Double_t sigma, Double_t sigma_n, Int_t i);
259
260   //------------------------------------------------------------------
261   /** 
262    * Evaluate 
263    * @f[ 
264    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
265    * @f] 
266    * 
267    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
268    * Landau with a Gaussian (see LandauGaus).  Note that 
269    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
270    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
271    *  
272    * References: 
273    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
274    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
275    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
276    * 
277    * @param x        Where to evaluate @f$ f_N@f$
278    * @param delta    @f$ \Delta_1@f$ 
279    * @param xi       @f$ \xi_1@f$
280    * @param sigma    @f$ \sigma_1@f$ 
281    * @param sigma_n  @f$ \sigma_n@f$ 
282    * @param n        @f$ N@f$ in the sum above.
283    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
284    *                 @f$ i > 1@f$ 
285    * 
286    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
287    */
288   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
289                               Double_t sigma, Double_t sigma_n, Int_t n, 
290                               const Double_t* a);
291   /** 
292    * Generate a TF1 object of @f$ f_I@f$ 
293    * 
294    * @param c        Constant
295    * @param delta    @f$ \Delta@f$ 
296    * @param xi       @f$ \xi_1@f$              
297    * @param sigma    @f$ \sigma_1@f$           
298    * @param sigma_n  @f$ \sigma_n@f$           
299    * @param i        @f$ i@f$ - the number of particles
300    * @param xmin     Least value of range
301    * @param xmax     Largest value of range
302    * 
303    * @return Newly allocated TF1 object
304    */
305   static TF1* MakeILandauGaus(Double_t c, 
306                               Double_t delta, Double_t xi, 
307                               Double_t sigma, Double_t sigma_n,
308                               Int_t    i, 
309                               Double_t xmin,  Double_t  xmax);
310   /** 
311    * Generate a TF1 object of @f$ f_N@f$ 
312    * 
313    * @param c         Constant                         
314    * @param delta     @f$ \Delta@f$                    
315    * @param xi        @f$ \xi_1@f$                     
316    * @param sigma     @f$ \sigma_1@f$                  
317    * @param sigma_n   @f$ \sigma_n@f$                  
318    * @param n         @f$ N@f$ - how many particles to sum to
319    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
320    *                  @f$ i > 1@f$ 
321    * @param xmin      Least value of range  
322    * @param xmax      Largest value of range
323    * 
324    * @return Newly allocated TF1 object
325    */
326   static TF1* MakeNLandauGaus(Double_t c, 
327                               Double_t delta, Double_t  xi, 
328                               Double_t sigma, Double_t  sigma_n,
329                               Int_t    n,     const Double_t* a, 
330                               Double_t xmin,  Double_t  xmax);
331                                                     
332   //__________________________________________________________________
333   /** 
334    * Structure to do fits to the energy loss spectrum 
335    * 
336    * @ingroup pwg2_forward 
337    */
338   struct ELossFitter 
339   {
340     enum { 
341       kC     = 0,
342       kDelta, 
343       kXi, 
344       kSigma, 
345       kSigmaN, 
346       kN, 
347       kA
348     };
349     /** 
350      * Constructor 
351      * 
352      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
353      * @param maxRange   Maximum range to fit to 
354      * @param minusBins  The number of bins below maximum to use 
355      */
356     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
357     /** 
358      * Destructor
359      * 
360      */
361     virtual ~ELossFitter();
362     /** 
363      * Clear internal arrays 
364      * 
365      */
366     void Clear();
367     /** 
368      * Fit a 1-particle signal to the passed energy loss distribution 
369      * 
370      * Note that this function clears the internal arrays first 
371      *
372      * @param dist    Data to fit the function to 
373      * @param sigman If larger than zero, the initial guess of the
374      *               detector induced noise. If zero or less, then this 
375      *               parameter is ignored in the fit (fixed at 0)
376      * 
377      * @return The function fitted to the data 
378      */
379     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
380     /** 
381      * Fit a N-particle signal to the passed energy loss distribution 
382      *
383      * If there's no 1-particle fit present, it does that first 
384      *
385      * @param dist   Data to fit the function to 
386      * @param n      Number of particle signals to fit 
387      * @param sigman If larger than zero, the initial guess of the
388      *               detector induced noise. If zero or less, then this 
389      *               parameter is ignored in the fit (fixed at 0)
390      * 
391      * @return The function fitted to the data 
392      */
393     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
394     /**
395      * Get Lower cut on data 
396      *
397      * @return Lower cut on data 
398      */
399     Double_t GetLowCut() const { return fLowCut; }
400     /**
401      * Get Maximum range to fit 
402      *
403      * @return Maximum range to fit 
404      */
405     Double_t GetMaxRange() const { return fMaxRange; }
406     /**
407      * Get Number of bins from maximum to fit 1st peak
408      *
409      * @return Number of bins from maximum to fit 1st peak
410      */
411     UShort_t GetMinusBins() const { return fMinusBins; }
412     /**
413      * Get Array of fit results 
414      *
415      * @return Array of fit results 
416      */
417     const TObjArray& GetFitResults() const { return fFitResults; }
418     /** 
419      * Get Array of fit results  
420      *
421      * @return Array of fit results 
422      */
423     TObjArray& GetFitResults() { return fFitResults; }
424     /**
425      * Get Array of functions 
426      *
427      * @return Array of functions 
428      */
429     const TObjArray& GetFunctions() const { return fFunctions; }
430     /** 
431      * Get Array of functions  
432      *
433      * @return Array of functions 
434      */
435     TObjArray& GetFunctions() { return fFunctions; }
436   private:
437     const Double_t fLowCut;     // Lower cut on data 
438     const Double_t fMaxRange;   // Maximum range to fit 
439     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
440     TObjArray fFitResults;      // Array of fit results 
441     TObjArray fFunctions;       // Array of functions 
442   };
443   /* @} */
444       
445
446   //==================================================================
447   /** 
448    * @{
449    * @name Convenience containers 
450    */
451   /** 
452    * Structure to hold histograms 
453    *
454    * @ingroup pwg2_forward 
455    */
456   struct Histos : public TObject
457   {     
458     /** 
459      * Constructor 
460      * 
461      * 
462      */
463     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
464     /** 
465      * Copy constructor 
466      * 
467      * @param o Object to copy from 
468      */
469     Histos(const Histos& o) 
470       : TObject(o), 
471         fFMD1i(o.fFMD1i), 
472         fFMD2i(o.fFMD2i), 
473         fFMD2o(o.fFMD2o), 
474         fFMD3i(o.fFMD3i), 
475         fFMD3o(o.fFMD3o)
476     {}
477     /** 
478      * Assignement operator 
479      * 
480      * @return Reference to this 
481      */
482     Histos& operator=(const Histos&) { return *this;}
483     /** 
484      * Destructor
485      */
486     ~Histos();
487     /** 
488      * Initialize the object 
489      * 
490      * @param etaAxis Eta axis to use 
491      */
492     void Init(const TAxis& etaAxis);
493     /** 
494      * Make a histogram 
495      * 
496      * @param d        Detector
497      * @param r        Ring 
498      * @param etaAxis  Eta axis to use
499      * 
500      * @return Newly allocated histogram 
501      */
502     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
503     /** 
504      * Clear data 
505      * 
506      * @param option Not used 
507      */
508     void  Clear(Option_t* option="");
509     // const TH2D* Get(UShort_t d, Char_t r) const;
510     /** 
511      * Get the histogram for a particular detector,ring
512      * 
513      * @param d Detector 
514      * @param r Ring 
515      * 
516      * @return Histogram for detector,ring or nul 
517      */
518     TH2D* Get(UShort_t d, Char_t r) const;
519     TH2D* fFMD1i; // Histogram for FMD1i
520     TH2D* fFMD2i; // Histogram for FMD2i
521     TH2D* fFMD2o; // Histogram for FMD2o
522     TH2D* fFMD3i; // Histogram for FMD3i
523     TH2D* fFMD3o; // Histogram for FMD3o
524
525     ClassDef(Histos,1) 
526   };
527
528   //__________________________________________________________________
529   /**
530    * Base class for structure holding ring specific histograms
531    * 
532    * @ingroup pwg2_forward 
533    */
534   struct RingHistos : public TObject
535   {
536     /** 
537      * Constructor
538      * 
539      */
540     RingHistos() : fDet(0), fRing('\0'), fName("") {}
541     /** 
542      * 
543      * 
544      * @param d Detector
545      * @param r Ring 
546      */
547     RingHistos(UShort_t d, Char_t r) 
548       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
549     {}
550     /** 
551      * Copy constructor
552      * 
553      * @param o Object to copy from 
554      */
555     RingHistos(const RingHistos& o) 
556       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
557     {}
558     /** 
559      * 
560      */
561     virtual ~RingHistos() {}
562     /** 
563      * Assignement operator
564      * 
565      * @param o Object to assign from
566      * 
567      * @return Reference to this
568      */
569     RingHistos& operator=(const RingHistos& o) 
570     {
571       TObject::operator=(o);
572       fDet  = o.fDet;
573       fRing = o.fRing;
574       fName = o.fName;
575       return *this;
576     }
577     /** 
578      * Define the outout list in @a d
579      * 
580      * @param d Where to put the output list
581      * 
582      * @return Newly allocated TList object or null
583      */
584     TList* DefineOutputList(TList* d) const;
585     /** 
586      * Get our output list from the container @a d
587      * 
588      * @param d where to get the output list from 
589      * 
590      * @return The found TList or null
591      */
592     TList* GetOutputList(const TList* d) const;
593     /** 
594      * Find a specific histogram in the source list @a d
595      * 
596      * @param d     (top)-container 
597      * @param name  Name of histogram
598      * 
599      * @return Found histogram or null
600      */
601     TH1* GetOutputHist(const TList* d, const char* name) const;
602     /** 
603      * 
604      * 
605      * 
606      * @return 
607      */
608     Color_t Color() const 
609     { 
610       return AliForwardUtil::RingColor(fDet, fRing);
611     }
612     const char* GetName() const { return fName.Data(); } 
613     UShort_t fDet;   // Detector
614     Char_t   fRing;  // Ring
615     TString  fName;  // Name
616
617     ClassDef(RingHistos,1) 
618   };
619   /* @} */
620 private:
621   AliForwardUtil() {}
622   AliForwardUtil(const AliForwardUtil& o) : TObject(o) {}
623   AliForwardUtil& operator=(const AliForwardUtil&) { return *this; }
624   ~AliForwardUtil() {}
625   
626
627   ClassDef(AliForwardUtil,1) // Utilities - do not make object
628 };
629
630 #endif
631 // Local Variables:
632 //  mode: C++
633 // End:
634