]>
Commit | Line | Data |
---|---|---|
239fe91c | 1 | // Author: Daniel.Lohner@cern.ch |
2 | ||
3 | #include "AliTRDNDFast.h" | |
4 | #include "TCanvas.h" | |
5 | ||
6 | extern Double_t langaufunc(Double_t *x, Double_t *par) { | |
7 | ||
8 | // needs protection, parameter [3]>0!!!!! | |
9 | // needs protection, parameter [4]>0!!!!! | |
10 | ||
11 | //Fit parameters: | |
12 | //par[0]=Width (scale) parameter of Landau density | |
13 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
14 | //par[2]=Total area (integral -inf to inf, normalization constant) | |
15 | //par[3]=Width (sigma) of convoluted Gaussian function | |
16 | //par[4]=Exponential Slope Parameter | |
17 | // | |
18 | //In the Landau distribution (represented by the CERNLIB approximation), | |
19 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
20 | //This shift is corrected within this function, so that the actual | |
21 | //maximum is identical to the MP parameter. | |
22 | ||
23 | // Numeric constants | |
24 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
25 | Double_t mpshift = -0.22278298; // Landau maximum location | |
26 | ||
27 | // Control constants | |
28 | Double_t np = 100.0; // number of convolution stpdf | |
29 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
30 | ||
31 | // Variables | |
32 | Double_t xx; | |
33 | Double_t mpc; | |
34 | Double_t fland; | |
35 | Double_t sum = 0.0; | |
36 | Double_t xlow,xupp; | |
37 | Double_t step; | |
38 | Double_t i; | |
39 | ||
40 | // MP shift correction | |
41 | mpc = par[1] - mpshift * par[0]; | |
42 | ||
43 | // Range of convolution integral | |
44 | xlow = x[0] - sc * par[3]; | |
45 | xupp = x[0] + sc * par[3]; | |
46 | ||
47 | if(xlow<0)xlow=0; | |
48 | if(xupp<xlow)cout<<"ERROR: Convolution around Negative MPV"<<endl; | |
49 | ||
50 | step = (xupp-xlow) / np; | |
51 | ||
52 | // Convolution integral of Landau and Gaussian by sum | |
53 | for(i=1.0; i<=np/2; i++) { | |
54 | xx = xlow + (i-.5) * step; | |
55 | fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0]; // mpc taken out | |
56 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
57 | ||
58 | xx = xupp - (i-.5) * step; | |
59 | fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0]; // mpc taken out | |
60 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
61 | } | |
62 | ||
63 | return (par[2] * step * sum * invsq2pi / par[3]); | |
64 | } | |
65 | ||
66 | ||
67 | ||
68 | ClassImp(AliTRDNDFast); | |
69 | ||
70 | AliTRDNDFast::AliTRDNDFast(): TObject(), | |
71 | fNDim(0), | |
72 | fTitle(""), | |
73 | fFunc(NULL), | |
74 | fHistos(NULL), | |
75 | fPars() | |
76 | { | |
77 | // default constructor | |
78 | } | |
79 | ||
80 | AliTRDNDFast::AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup): TObject(), | |
81 | fNDim(ndim), | |
82 | fTitle(name), | |
83 | fFunc(NULL), | |
84 | fHistos(NULL), | |
85 | fPars() | |
86 | { | |
87 | Init(); | |
88 | for(Int_t idim=0;idim<fNDim;idim++){ | |
89 | fHistos[idim]=new TH1F(Form("%s_axis_%d",fTitle.Data(),idim),Form("%s_axis_%d",fTitle.Data(),idim),nbins,xlow,xup); | |
90 | fHistos[idim]->SetDirectory(0); | |
91 | } | |
92 | } | |
93 | ||
94 | AliTRDNDFast::AliTRDNDFast(const AliTRDNDFast &ref) : TObject(ref), | |
95 | fNDim(ref.fNDim), | |
96 | fTitle(ref.fTitle), | |
97 | fFunc(NULL), | |
98 | fHistos(NULL), | |
99 | fPars() | |
100 | { | |
101 | Cleanup(); | |
102 | Init(); | |
103 | for(Int_t idim=0;idim<fNDim;idim++){ | |
104 | fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim)); | |
105 | fHistos[idim]->SetDirectory(0); | |
106 | for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar]; | |
107 | } | |
108 | } | |
109 | ||
110 | AliTRDNDFast &AliTRDNDFast::operator=(const AliTRDNDFast &ref){ | |
111 | // | |
112 | // Assignment operator | |
113 | // | |
114 | if(this != &ref){ | |
115 | TObject::operator=(ref); | |
116 | fNDim=ref.fNDim; | |
117 | fTitle=ref.fTitle; | |
118 | fFunc = ref.fFunc; | |
239fe91c | 119 | for(Int_t idim=0;idim<fNDim;idim++){ |
84f3d859 | 120 | fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim)); |
121 | fHistos[idim]->SetDirectory(0); | |
122 | for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar]; | |
239fe91c | 123 | } |
124 | } | |
125 | return *this; | |
126 | } | |
127 | ||
128 | AliTRDNDFast::~AliTRDNDFast(){ | |
129 | Cleanup(); | |
130 | ||
131 | } | |
132 | ||
133 | void AliTRDNDFast::Init(){ | |
134 | ||
135 | for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].Set(fNDim); | |
136 | fFunc=new TF1*[fNDim]; | |
137 | fHistos=new TH1F*[fNDim]; | |
138 | for(Int_t idim=0;idim<fNDim;idim++){ | |
139 | fHistos[idim]=NULL; | |
140 | for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(0,idim); | |
141 | fFunc[idim]=NULL; | |
142 | } | |
143 | } | |
144 | ||
145 | void AliTRDNDFast::Cleanup(){ | |
146 | if(fHistos){ | |
147 | for(Int_t idim=0;idim<fNDim;idim++){ | |
148 | if(fHistos[idim]){ | |
149 | delete fHistos[idim]; | |
150 | fHistos[idim]=NULL; | |
151 | } | |
152 | } | |
153 | delete[] fHistos; | |
154 | fHistos=NULL; | |
155 | } | |
156 | if(fPars){ | |
157 | for(Int_t ipar=0;ipar<kNpar;ipar++){ | |
158 | fPars[ipar].Reset(); | |
159 | } | |
160 | } | |
161 | if(fFunc){ | |
162 | for(Int_t idim=0;idim<fNDim;idim++){ | |
163 | if(fFunc[idim]){ | |
164 | delete fFunc[idim]; | |
165 | fFunc[idim]=NULL; | |
166 | } | |
167 | } | |
168 | delete[] fFunc; | |
169 | fFunc=NULL; | |
170 | } | |
171 | } | |
172 | ||
173 | void AliTRDNDFast::PrintPars(){ | |
174 | for(Int_t idim=0;idim<fNDim;idim++){ | |
175 | for(Int_t ipar=0;ipar<kNpar;ipar++)cout<<idim<<" "<<ipar<<" "<<GetParam(idim,ipar)<<endl; | |
176 | } | |
177 | } | |
178 | ||
179 | void AliTRDNDFast::ScaleLangauFun(TF1 *func,Double_t mpv){ | |
180 | ||
181 | Double_t start[kNpar],low[kNpar],high[kNpar]; | |
182 | for(Int_t ii=0;ii<kNpar;ii++){ | |
183 | start[ii]=func->GetParameter(ii); | |
184 | func->GetParLimits(ii,low[ii],high[ii]); | |
185 | } | |
186 | ||
187 | Double_t scalefactor=mpv/start[1]; | |
188 | ||
189 | for(Int_t ii=0;ii<kNpar;ii++){ | |
190 | Double_t scale=1; | |
191 | if(ii==0||ii==1||ii==3)scale=scalefactor; | |
192 | if(ii==4)scale=1./scalefactor; | |
193 | start[ii]*=scale; | |
194 | low[ii]*=scale; | |
195 | high[ii]*=scale; | |
196 | func->SetParLimits(ii, low[ii], high[ii]); | |
197 | } | |
198 | func->SetParameters(start); | |
199 | } | |
200 | ||
201 | //--------------------------------------------------------------- | |
202 | //--------------------------------------------------------------- | |
203 | TF1 *AliTRDNDFast::GetLangauFun(TString funcname,Double_t range[2],Double_t scalefactor){ | |
204 | ||
205 | Double_t start[kNpar] = {120, 1000, 1, 100, 1e-5}; | |
206 | Double_t low[kNpar] = {10, 300, 0.01, 1, 0.}; | |
207 | Double_t high[kNpar] = {1000, 3000, 100, 1000, 1.}; | |
208 | ||
209 | TF1 *hlandauE=new TF1(funcname.Data(),langaufunc,0,range[1],kNpar); | |
210 | hlandauE->SetParameters(start); | |
211 | hlandauE->SetParNames("Width","MP","Area","GSigma","delta"); | |
212 | for (int i=0; i<kNpar; i++) { | |
213 | hlandauE->SetParLimits(i, low[i], high[i]); | |
214 | } | |
215 | ||
216 | if(scalefactor!=1){ScaleLangauFun(hlandauE,scalefactor*start[1]);} | |
217 | ||
218 | return hlandauE; | |
219 | } | |
220 | ||
221 | TF1 *AliTRDNDFast::FitLandau(TString name,TH1F *htemp,Double_t range[2],TString option){ | |
222 | ||
223 | TF1 *fitlandau1D=GetLangauFun(name,range); | |
224 | TF1 fit("land","landau"); | |
225 | Double_t max=htemp->GetXaxis()->GetBinCenter(htemp->GetMaximumBin()); | |
226 | fit.SetParameter(1,max); | |
227 | fit.SetParLimits(1,0,htemp->GetXaxis()->GetXmax()); | |
228 | fit.SetParameter(2,0.3*max); // MPV may never become negative!!!!! | |
229 | htemp->Fit("land",option.Data(),"",range[0],range[1]); | |
230 | ScaleLangauFun(fitlandau1D,fit.GetParameter(1)); | |
231 | htemp->Fit(fitlandau1D,option.Data(),"",range[0],range[1]); // range for references | |
232 | return fitlandau1D; | |
233 | } | |
234 | ||
235 | void AliTRDNDFast::BuildHistos(){ | |
236 | ||
237 | for(Int_t idim=0;idim<fNDim;idim++){ | |
238 | fHistos[idim]->Reset(); | |
239 | // Fill Histo | |
240 | Double_t pars[kNpar]; | |
241 | for(Int_t ipar=0;ipar<kNpar;ipar++)pars[ipar]=GetParam(idim,ipar); | |
242 | // Also Fill overflow bins!!! | |
d2a461e9 | 243 | for(Int_t ii=1;ii<=fHistos[idim]->GetNbinsX()+1;ii++){ |
239fe91c | 244 | Double_t xval=fHistos[idim]->GetXaxis()->GetBinCenter(ii); |
245 | Double_t val=langaufunc(&xval,&pars[0]); | |
246 | //Double_t val=fFunc[idim]->Eval(xval); | |
247 | fHistos[idim]->SetBinContent(ii,val); | |
248 | } | |
249 | // Normalization | |
250 | fHistos[idim]->Scale(1./fHistos[idim]->Integral(1,fHistos[idim]->GetNbinsX(),"width")); | |
251 | } | |
252 | } | |
253 | ||
254 | void AliTRDNDFast::Build(Double_t **pars){ | |
255 | // pars[ndim][npar] | |
256 | for(Int_t idim=0;idim<fNDim;idim++){ | |
257 | for(Int_t ipar=0;ipar<kNpar;ipar++){ | |
258 | fPars[ipar].SetAt(pars[idim][ipar],idim); | |
259 | } | |
260 | } | |
261 | BuildHistos(); | |
262 | } | |
263 | ||
264 | ||
265 | void AliTRDNDFast::Build(TH1F **hdEdx,TString path){ | |
266 | ||
267 | Double_t range[2]; | |
268 | ||
269 | TCanvas *CANV=new TCanvas("a","a"); | |
270 | if(fNDim==2)CANV->Divide(2,1); | |
271 | if(fNDim==3)CANV->Divide(2,2); | |
272 | if(fNDim>6)CANV->Divide(3,3); | |
273 | // Fit and Extract Parameters | |
274 | for(Int_t idim=0;idim<fNDim;idim++){ | |
275 | CANV->cd(idim+1); | |
276 | gPad->SetLogy(); | |
277 | range[0]=hdEdx[idim]->GetXaxis()->GetXmin(); | |
278 | range[1]=hdEdx[idim]->GetXaxis()->GetXmax(); | |
279 | // Norm Histogram | |
280 | hdEdx[idim]->Scale(1./hdEdx[idim]->Integral(1,hdEdx[idim]->GetNbinsX(),"width")); | |
281 | // Fit Histogram | |
282 | fFunc[idim]=FitLandau(Form("fit%d",idim),hdEdx[idim],range,"RMB0"); | |
283 | // Norm Landau | |
284 | fFunc[idim]->SetParameter(2,fFunc[idim]->GetParameter(2)/fFunc[idim]->Integral(range[0],range[1])); | |
285 | hdEdx[idim]->DrawCopy(); | |
286 | fFunc[idim]->DrawCopy("same"); | |
287 | // Set Pars | |
288 | for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(fFunc[idim]->GetParameter(ipar),idim); | |
289 | } | |
290 | if(path!="")CANV->Print(Form("%s/%s_Build.pdf",path.Data(),fTitle.Data())); | |
291 | delete CANV; | |
292 | ||
293 | BuildHistos(); | |
294 | } | |
295 | ||
296 | Double_t AliTRDNDFast::Eval(Double_t *point) const{ | |
297 | Double_t val=1; | |
298 | for(Int_t idim=0;idim<fNDim;idim++){ | |
299 | Int_t bin=fHistos[idim]->GetXaxis()->FindBin(point[idim]); | |
300 | val*=fHistos[idim]->GetBinContent(bin); | |
301 | } | |
302 | return val; | |
303 | } | |
304 | ||
305 | void AliTRDNDFast::Random(Double_t *point) const{ | |
306 | for(Int_t idim=0;idim<fNDim;idim++){ | |
307 | point[idim]=fHistos[idim]->GetRandom(); | |
308 | } | |
309 | } | |
310 | ||
311 | void AliTRDNDFast::Random(Double_t *point,AliTRDNDFast *nd0,AliTRDNDFast *nd1,Double_t w0,Double_t w1){ | |
312 | for(Int_t idim=0;idim<nd0->fNDim;idim++){ | |
313 | point[idim]=GetRandomInterpolation(nd0->fHistos[idim],nd1->fHistos[idim],w0,w1); | |
314 | } | |
315 | } | |
316 | ||
317 | Int_t AliTRDNDFast::BinarySearchInterpolation(Int_t start,Int_t end,Double_t *a0,Double_t *a1,Double_t w0,Double_t w1,Double_t val){ | |
318 | ||
319 | if((end-start)==1)return start; | |
320 | Int_t mid=Int_t((end+start)/2); | |
321 | Double_t valmid=(w0*a0[mid]+w1*a1[mid])/(w0+w1); | |
322 | if(val>=valmid)return BinarySearchInterpolation(mid,end,a0,a1,w0,w1,val); | |
323 | return BinarySearchInterpolation(start,mid,a0,a1,w0,w1,val); | |
324 | } | |
325 | ||
326 | Double_t AliTRDNDFast::GetRandomInterpolation(TH1F *hist0,TH1F *hist1,Double_t w0,Double_t w1){ | |
327 | ||
328 | // Draw Random Variable | |
329 | Double_t rand=gRandom->Rndm(); | |
330 | ||
331 | // Get Integral Arrays | |
332 | Double_t *integral0=hist0->GetIntegral(); | |
333 | Double_t *integral1=hist1->GetIntegral(); | |
334 | ||
335 | Int_t nbinsX=hist0->GetNbinsX(); | |
336 | ||
337 | Int_t ibin=BinarySearchInterpolation(1,nbinsX+1,integral0,integral1,w0,w1,rand); | |
338 | return hist0->GetXaxis()->GetBinCenter(ibin); | |
339 | } | |
340 | ||
341 | ||
342 |