Corrected protection.
[u/mrichter/AliRoot.git] / STEER / AliITSPidParams.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 //                                                               //
20 // Implementation of the class to store parameters of ITS        //
21 // response funcions for dE/dx based PID                         //
22 // Origin: F.Prino, Torino, prino@to.infn.it                     //
23 //                                                               //
24 ///////////////////////////////////////////////////////////////////
25
26 #include <TFormula.h>
27 #include <TNamed.h>
28 #include <TMath.h>
29 #include "AliITSPidParams.h"
30 #include "AliPID.h"
31
32 ClassImp(AliITSPidParams)
33
34 //______________________________________________________________________
35 AliITSPidParams::AliITSPidParams():
36   TNamed("default",""),
37   fSDDPionMPV(0),
38   fSDDPionLandauWidth(0),
39   fSDDPionGaussWidth(0),
40   fSSDPionMPV(0),
41   fSSDPionLandauWidth(0),
42   fSSDPionGaussWidth(0),
43   fSDDKaonMPV(0),
44   fSDDKaonLandauWidth(0),
45   fSDDKaonGaussWidth(0),
46   fSSDKaonMPV(0),
47   fSSDKaonLandauWidth(0),
48   fSSDKaonGaussWidth(0),
49   fSDDProtMPV(0),
50   fSDDProtLandauWidth(0),
51   fSDDProtGaussWidth(0),
52   fSSDProtMPV(0),
53   fSSDProtLandauWidth(0),
54   fSSDProtGaussWidth(0)
55 {
56   // default constructor
57   InitDefaults();
58 }
59 //______________________________________________________________________
60 AliITSPidParams::AliITSPidParams(Char_t * name):
61   TNamed(name,""),
62   fSDDPionMPV(0),
63   fSDDPionLandauWidth(0),
64   fSDDPionGaussWidth(0),
65   fSSDPionMPV(0),
66   fSSDPionLandauWidth(0),
67   fSSDPionGaussWidth(0),
68   fSDDKaonMPV(0),
69   fSDDKaonLandauWidth(0),
70   fSDDKaonGaussWidth(0),
71   fSSDKaonMPV(0),
72   fSSDKaonLandauWidth(0),
73   fSSDKaonGaussWidth(0),
74   fSDDProtMPV(0),
75   fSDDProtLandauWidth(0),
76   fSDDProtGaussWidth(0),
77   fSSDProtMPV(0),
78   fSSDProtLandauWidth(0),
79   fSSDProtGaussWidth(0)
80 {
81   // standard constructor
82   InitDefaults();
83 }
84 //______________________________________________________________________
85 AliITSPidParams::~AliITSPidParams(){
86   // 
87   if(fSDDPionMPV) delete fSDDPionMPV;
88   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
89   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
90
91   if(fSSDPionMPV) delete fSSDPionMPV;
92   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
93   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
94
95   if(fSDDKaonMPV) delete fSDDKaonMPV;
96   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
97   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
98
99   if(fSSDKaonMPV) delete fSSDKaonMPV;
100   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
101   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
102
103   if(fSDDProtMPV) delete fSDDProtMPV;
104   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
105   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
106
107   if(fSSDProtMPV) delete fSSDProtMPV;
108   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
109   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
110 }
111
112 //______________________________________________________________________
113 void AliITSPidParams::InitDefaults(){
114   // initialize TFormulas to default values (=p-p simulations PYTHIA+GEANT)
115   // parameter values from Emanuele Biolcati
116
117   // pions
118   if(fSDDPionMPV) delete fSDDPionMPV;
119   fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]*TMath::Log(x*x)+[3]");
120   fSDDPionMPV->SetParameters(1.4831933,-0.000403,3.9722756,92.710680);
121
122   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
123   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
124   fSDDPionLandauWidth->SetParameters(-0.091098,-0.001355,8.3019280);
125
126   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
127   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
128   fSDDPionGaussWidth->SetParameters(-0.129570,-0.002686,15.287701);
129
130   if(fSSDPionMPV) delete fSSDPionMPV;
131   fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]*TMath::Log(x*x)+[3]");
132   fSSDPionMPV->SetParameters(1.2455667,-0.000743,3.2260119,84.237030);
133
134   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
135   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
136   fSSDPionLandauWidth->SetParameters(-0.083918,-0.001246,7.4750478);
137
138   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
139   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
140   fSSDPionGaussWidth->SetParameters(-0.140441,-0.002984,10.936747);
141
142   // kaons
143   if(fSDDKaonMPV) delete fSDDKaonMPV;
144   fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
145   fSDDKaonMPV->SetParameters(16.998674,-0.058006,89.660709);
146
147   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
148   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
149   fSDDKaonLandauWidth->SetParameters(1.0014235,-5.61e-15,9.2672327);
150
151   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
152   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
153   fSDDKaonGaussWidth->SetParameters(0.7919025,-0.103579,14.016803);
154
155   if(fSSDKaonMPV) delete fSSDKaonMPV;
156   fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
157   fSSDKaonMPV->SetParameters(14.090845,-0.087253,81.782088);
158
159   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
160   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
161   fSSDKaonLandauWidth->SetParameters(1.0769127,-9.06e-13,7.5221492);
162
163   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
164   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
165   fSSDKaonGaussWidth->SetParameters(0.5155878,-0.098696,10.771975);
166
167   // protons
168   if(fSDDProtMPV) delete fSDDProtMPV;
169   fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
170   fSDDProtMPV->SetParameters(70.325146,0.0386808,87.797052);
171
172   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
173   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
174   fSDDProtLandauWidth->SetParameters(4.0476840,-3.77e-14,10.294707);
175
176   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
177   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
178   fSDDProtGaussWidth->SetParameters(5.9780498,-2.16e-11,13.357744);
179
180   if(fSSDProtMPV) delete fSSDProtMPV;
181   fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
182   fSSDProtMPV->SetParameters(58.918831,-0.303855,80.341765);
183                              
184   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
185   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
186   fSSDProtLandauWidth->SetParameters(3.0814273,-1.26e-13,8.8353833);
187
188   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
189   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x*x)+[2]");
190   fSSDProtGaussWidth->SetParameters(5.6177229,-1.67e-13,10.538921);
191
192 }
193 //_______________________________________________________________________
194 Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
195   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
196   if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
197   else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
198   else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
199   else return 0.;
200 }
201 //_______________________________________________________________________
202   Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
203   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
204
205   Double_t par[3];
206   Bool_t isSet=kFALSE;
207   if(partType==AliPID::kPion){
208     if(lay==3 || lay==4){
209       par[0]=GetSDDPionLandauWidth(mom);
210       par[1]=GetSDDPionMPV(mom);
211       par[2]=GetSDDPionGaussWidth(mom);
212       isSet=kTRUE;
213     }
214     else if(lay==5 || lay==6){
215       par[0]=GetSSDPionLandauWidth(mom);
216       par[1]=GetSSDPionMPV(mom);
217       par[2]=GetSSDPionGaussWidth(mom);
218       isSet=kTRUE;
219     }
220   }else if(partType==AliPID::kKaon){
221     if(lay==3 || lay==4){
222       par[0]=GetSDDKaonLandauWidth(mom);
223       par[1]=GetSDDKaonMPV(mom);
224       par[2]=GetSDDKaonGaussWidth(mom);
225       isSet=kTRUE;
226     }
227     else if(lay==5 || lay==6){
228       par[0]=GetSSDKaonLandauWidth(mom);
229       par[1]=GetSSDKaonMPV(mom);
230       par[2]=GetSSDKaonGaussWidth(mom);
231       isSet=kTRUE;
232     }
233   }else if(partType==AliPID::kProton){
234     if(lay==3 || lay==4){
235       par[0]=GetSDDProtLandauWidth(mom);
236       par[1]=GetSDDProtMPV(mom);
237       par[2]=GetSDDProtGaussWidth(mom);
238       isSet=kTRUE;
239     }
240     else if(lay==5 || lay==6){
241       par[0]=GetSSDProtLandauWidth(mom);
242       par[1]=GetSSDProtMPV(mom);
243       par[2]=GetSSDProtGaussWidth(mom);
244       isSet=kTRUE;
245     }
246   }
247   if(!isSet) return 0.;
248   // Numeric constants
249   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
250   Double_t mpshift  = -0.22278298;       // Landau maximum location
251   // Control constants
252   Double_t np = 100.0;      // number of convolution steps
253   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
254   // Variables
255   Double_t xx;
256   Double_t mpc;
257   Double_t fland;
258   Double_t sum = 0.0;
259   Double_t xlow,xupp;
260   Double_t step;
261   Double_t i;
262   
263   // MP shift correction
264   mpc = par[1] - mpshift * par[0];
265   // Range of convolution integral
266   xlow = dedx - sc * par[2];
267   xupp = dedx + sc * par[2];
268   if(np!=0) step = (xupp-xlow) / np;
269   
270   // Convolution integral of Landau and Gaussian by sum
271   for(i=1.0; i<=np/2; i++) {
272     xx = xlow + (i-.5) * step;
273    
274     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
275     sum += fland * TMath::Gaus(dedx,xx,par[2]);
276     
277     xx = xupp - (i-.5) * step;
278     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
279     sum += fland * TMath::Gaus(dedx,xx,par[2]);
280   }
281   
282   return (step * sum * invsq2pi / par[2]);
283 }
284