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