Inserting TMath.h where required by the new version of ROOT
[u/mrichter/AliRoot.git] / ITS / AliITSPidParItem.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /////////////////////////////////////////////////////////
19 //Class for PID in the ITS                             //
20 //                                                     //
21 //                                                     //
22 /////////////////////////////////////////////////////////
23
24 #include <Riostream.h>
25 #include <TF1.h>
26 #include <TMath.h>
27
28 #include "AliITSPidParItem.h"
29
30 ClassImp(AliITSPidParItem)
31 //____________________________________________________________________
32 AliITSPidParItem::AliITSPidParItem():
33 fPCenter(0),
34 fPWidth(0)
35 {
36   // default constructor
37   for(Int_t i=0;i<39;i++){
38   fBuff[i]=0;
39   }
40 }//____________________________________________________________________
41 AliITSPidParItem::AliITSPidParItem(Float_t center,Float_t width,Double_t *buff):
42 fPCenter(center),
43 fPWidth(width)
44 {
45   // standard constructor
46   for (Int_t i=0;i<39;i++) fBuff[i]=buff[i];
47
48 }
49
50 //____________________________________________________________________
51 void AliITSPidParItem::GetParameters(Double_t *buff) const{
52   //get all the parameters
53   for (Int_t i=0;i<39;i++) buff[i]=fBuff[i];
54
55 }
56 //____________________________________________________________________
57 void AliITSPidParItem::GetProtonPar(Double_t *buffp) const{
58   //get the protons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
59   buffp[0]=fBuff[0];
60   buffp[1]=fBuff[1];
61   buffp[2]=fBuff[2];
62   buffp[3]=fBuff[3];
63   buffp[4]=fBuff[8];
64   buffp[5]=fBuff[9];
65   buffp[6]=fBuff[10];
66 }
67 //____________________________________________________________________
68 void AliITSPidParItem::GetKaonPar(Double_t *buffk) const{
69   //get the kaons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
70   buffk[0]=fBuff[13];
71   buffk[1]=fBuff[14];
72   buffk[2]=fBuff[15];
73   buffk[3]=fBuff[16];
74   buffk[4]=fBuff[21];
75   buffk[5]=fBuff[22];
76   buffk[6]=fBuff[23];
77 }
78 //____________________________________________________________________
79 void AliITSPidParItem::GetPionPar(Double_t *buffpi) const{
80   //get the pions' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
81   buffpi[0]=fBuff[26];
82   buffpi[1]=fBuff[27];
83   buffpi[2]=fBuff[28];
84   buffpi[3]=fBuff[29];
85   buffpi[4]=fBuff[34];
86   buffpi[5]=fBuff[35];
87   buffpi[6]=fBuff[36];
88 }
89 //____________________________________________________________________
90 void AliITSPidParItem::GetPar0(Double_t *buff0) const{
91   //Width Landau for protons, kaons, pions.
92   buff0[0]=fBuff[0];
93   buff0[1]=fBuff[13];
94   buff0[2]=fBuff[26];
95 }
96 //____________________________________________________________________
97 void AliITSPidParItem::GetPar1(Double_t *buff1) const{
98   //Most Probable for protons, kaons, pions.
99   buff1[0]=fBuff[1];
100   buff1[1]=fBuff[14];
101   buff1[2]=fBuff[27];
102 }
103 //____________________________________________________________________
104 void AliITSPidParItem::GetPar2(Double_t *buff2) const{
105   //Area for protons, kaons, pions.
106   buff2[0]=fBuff[2];
107   buff2[1]=fBuff[15];
108   buff2[2]=fBuff[28];
109 }
110 //____________________________________________________________________
111 void AliITSPidParItem::GetPar3(Double_t *buff3) const{
112   //Width Gaussian for protons, kaons, pions.
113   buff3[0]=fBuff[3];
114   buff3[1]=fBuff[16];
115   buff3[2]=fBuff[29];
116 }
117 //____________________________________________________________________
118 void AliITSPidParItem::GetChisquare(Double_t *buffchi) const{
119   //Chi2 of the fit for protons, kaons, pions.
120   buffchi[0]=fBuff[8];
121   buffchi[1]=fBuff[21];
122   buffchi[2]=fBuff[34];
123 }
124 //____________________________________________________________________
125 void AliITSPidParItem::GetNDF(Double_t *buffndf) const{
126   //NDF of the fit for protons, kaons, pions.
127   buffndf[0]=fBuff[9];
128   buffndf[1]=fBuff[22];
129   buffndf[2]=fBuff[35];
130 }
131 //____________________________________________________________________
132 void AliITSPidParItem::GetProParFun(Double_t *pfun) const{
133   //some Protons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
134   pfun[0]=fBuff[0];
135   pfun[1]=fBuff[1];
136   pfun[2]=fBuff[2];
137   pfun[3]=fBuff[3];
138   pfun[4]=fBuff[10];
139 }
140 //____________________________________________________________________
141 void AliITSPidParItem::GetKaoParFun(Double_t *kfun) const{
142   //some Kaons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
143   kfun[0]=fBuff[13];
144   kfun[1]=fBuff[14];
145   kfun[2]=fBuff[15];
146   kfun[3]=fBuff[16];
147   kfun[4]=fBuff[23];
148 }
149 //____________________________________________________________________
150 void AliITSPidParItem::GetPiParFun(Double_t *pifun) const{
151   //some Pions parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
152   pifun[0]=fBuff[26];
153   pifun[1]=fBuff[27];
154   pifun[2]=fBuff[28];
155   pifun[3]=fBuff[29];
156   pifun[4]=fBuff[36];
157 }
158 //____________________________________________________________________
159 void AliITSPidParItem::GetRangeLim(Double_t *range) const{
160   //Range limits for the response functions
161   range[0]=fBuff[11];//proton low
162   range[1]=fBuff[12];//proton high
163   range[2]=fBuff[24];//kaon low
164   range[3]=fBuff[25];//kaon high
165   range[4]=fBuff[37];//pion low
166   range[5]=fBuff[38];//pion high
167 }
168 //____________________________________________________________________
169 void AliITSPidParItem::GetProtonParErr(Double_t *bufferp)const{
170   //errors on the protons' parameters
171   bufferp[0]=fBuff[4];
172   bufferp[1]=fBuff[5];
173   bufferp[2]=fBuff[6];
174   bufferp[3]=fBuff[7];
175 }
176 //____________________________________________________________________
177 void AliITSPidParItem::GetKaonParErr(Double_t *bufferk)const{
178   //errors on the kaons' parameters
179   bufferk[0]=fBuff[17];
180   bufferk[1]=fBuff[18];
181   bufferk[2]=fBuff[19];
182   bufferk[3]=fBuff[20];
183 }
184 //____________________________________________________________________
185 void AliITSPidParItem::GetPionParErr(Double_t *bufferpi)const{
186   //errors on the pions' parameters
187   bufferpi[0]=fBuff[30];
188   bufferpi[1]=fBuff[31];
189   bufferpi[2]=fBuff[32];
190   bufferpi[3]=fBuff[33];
191 }
192 //____________________________________________________________________
193 void AliITSPidParItem::PrintParameters() const {
194   // Prints the data members of this class
195 cout<<"==============================***************"<<endl;
196 cout<<"Momentum (GeV/c) - center of the bin - "<<fPCenter<<endl;
197 cout<<" Width of momentum bin (GeV/c) "<<fPWidth<<endl;
198 for (Int_t i=0;i<39;i++) cout<<"Parameter"<<i<<" = "<<fBuff[i]<<endl;
199
200 }
201
202 //_______________________________________________________________________
203 TF1* AliITSPidParItem::CookFunIts(TString namefun,Double_t *par,Double_t rangei,Double_t rangef,TString comment){
204   //
205   TF1 *fun;
206     if (par[4]!=0) {
207       fun=new TF1(comment,Langaufun2,rangei,rangef,5);
208       fun->SetParameters(par);
209      
210     }
211     else {fun=new TF1(namefun,"0");}
212     return fun;
213 }
214
215 //_________________________________________________________________________
216 Double_t AliITSPidParItem::Langaufun(Double_t *x, Double_t *par) {
217
218   //Fit parameters:
219   //par[0]=Width (scale) parameter of Landau density
220   //par[1]=Most Probable (MP, location) parameter of Landau density
221   //par[2]=Total area (integral -inf to inf, normalization constant)
222   //par[3]=Width (sigma) of convoluted Gaussian function
223   //
224   //In the Landau distribution (represented by the CERNLIB approximation), 
225   //the maximum is located at x=-0.22278298 with the location parameter=0.
226   //This shift is corrected within this function, so that the actual
227   //maximum is identical to the MP parameter.
228
229   // Numeric constants
230   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
231   Double_t mpshift  = -0.22278298;       // Landau maximum location
232
233   // Control constants
234   Double_t np = 100.0;      // number of convolution steps
235   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
236
237   // Variables
238   Double_t xx;
239   Double_t mpc;
240   Double_t fland;
241   Double_t sum = 0.0;
242   Double_t xlow,xupp;
243   Double_t step;
244   Double_t i;
245
246
247   // MP shift correction
248   mpc = par[1] - mpshift * par[0]; 
249
250   // Range of convolution integral
251   xlow = x[0] - sc * par[3];
252   xupp = x[0] + sc * par[3];
253
254   step = (xupp-xlow) / np;
255
256   // Convolution integral of Landau and Gaussian by sum
257   for(i=1.0; i<=np/2; i++) {
258     xx = xlow + (i-.5) * step;
259     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
260     sum += fland * TMath::Gaus(x[0],xx,par[3]);
261
262     xx = xupp - (i-.5) * step;
263     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
264     sum += fland * TMath::Gaus(x[0],xx,par[3]);
265   }
266
267   return (par[2] * step * sum * invsq2pi / par[3]);
268 }
269 //_______________________________________________________________________
270 Double_t AliITSPidParItem::Langaufun2(Double_t *x, Double_t *par){
271   //
272   return 1/par[4]*Langaufun(x,par);
273 }