Added V0A23 (V0 rings 2-3), V0C01 (V0 rings 0-1) and V0S = V0A23+V0C01
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDNDFast.cxx
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;
119         for(Int_t idim=0;idim<fNDim;idim++){
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];
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!!!
243         for(Int_t ii=1;ii<=fHistos[idim]->GetNbinsX()+1;ii++){
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