Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDNDFast.cxx
CommitLineData
239fe91c 1// Author: Daniel.Lohner@cern.ch
2
3#include "AliTRDNDFast.h"
4#include "TCanvas.h"
5
6extern 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
68ClassImp(AliTRDNDFast);
69
70AliTRDNDFast::AliTRDNDFast(): TObject(),
71 fNDim(0),
72 fTitle(""),
73 fFunc(NULL),
74 fHistos(NULL),
75 fPars()
76{
77 // default constructor
78}
79
80AliTRDNDFast::AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup): TObject(),
81fNDim(ndim),
82fTitle(name),
83fFunc(NULL),
84fHistos(NULL),
85fPars()
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
94AliTRDNDFast::AliTRDNDFast(const AliTRDNDFast &ref) : TObject(ref),
95fNDim(ref.fNDim),
96fTitle(ref.fTitle),
97fFunc(NULL),
98fHistos(NULL),
99fPars()
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
110AliTRDNDFast &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
128AliTRDNDFast::~AliTRDNDFast(){
129 Cleanup();
130
131}
132
133void 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
145void 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
173void 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
179void 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//---------------------------------------------------------------
203TF1 *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
221TF1 *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
235void 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
254void 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
265void 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
296Double_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
305void AliTRDNDFast::Random(Double_t *point) const{
306 for(Int_t idim=0;idim<fNDim;idim++){
307 point[idim]=fHistos[idim]->GetRandom();
308 }
309}
310
311void 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
317Int_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
326Double_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