[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDNDFast.cxx

// Author: Daniel.Lohner@cern.ch

#include "AliTRDNDFast.h"
#include "TCanvas.h"

extern Double_t langaufunc(Double_t *x, Double_t *par) {

    // needs protection, parameter [3]>0!!!!!
    // needs protection, parameter [4]>0!!!!!

    //Fit parameters:
    //par[0]=Width (scale) parameter of Landau density
    //par[1]=Most Probable (MP, location) parameter of Landau density
    //par[2]=Total area (integral -inf to inf, normalization constant)
    //par[3]=Width (sigma) of convoluted Gaussian function
    //par[4]=Exponential Slope Parameter
    //
    //In the Landau distribution (represented by the CERNLIB approximation),
    //the maximum is located at x=-0.22278298 with the location parameter=0.
    //This shift is corrected within this function, so that the actual
    //maximum is identical to the MP parameter.

    // Numeric constants
    Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
    Double_t mpshift  = -0.22278298;       // Landau maximum location

    // Control constants
    Double_t np = 100.0;      // number of convolution stpdf
    Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

    // Variables
    Double_t xx;
    Double_t mpc;
    Double_t fland;
    Double_t sum = 0.0;
    Double_t xlow,xupp;
    Double_t step;
    Double_t i;

    // MP shift correction
    mpc = par[1] - mpshift * par[0];

    // Range of convolution integral
    xlow = x[0] - sc * par[3];
    xupp = x[0] + sc * par[3];

    if(xlow<0)xlow=0;
    if(xupp<xlow)cout<<"ERROR: Convolution around Negative MPV"<<endl;

    step = (xupp-xlow) / np;

    // Convolution integral of Landau and Gaussian by sum
    for(i=1.0; i<=np/2; i++) {
	xx = xlow + (i-.5) * step;
	fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];	// mpc taken out
	sum += fland * TMath::Gaus(x[0],xx,par[3]);

	xx = xupp - (i-.5) * step;
	fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];	// mpc taken out
	sum += fland * TMath::Gaus(x[0],xx,par[3]);
    }

    return (par[2] * step * sum * invsq2pi / par[3]);
}


ClassImp(AliTRDNDFast);

AliTRDNDFast::AliTRDNDFast(): TObject(),
    fNDim(0),
    fTitle(""),
    fFunc(NULL),
    fHistos(NULL),
    fPars()
{
    // default constructor
}

AliTRDNDFast::AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup): TObject(),
fNDim(ndim),
fTitle(name),
fFunc(NULL),
fHistos(NULL),
fPars()
{
    Init();
    for(Int_t idim=0;idim<fNDim;idim++){
	fHistos[idim]=new TH1F(Form("%s_axis_%d",fTitle.Data(),idim),Form("%s_axis_%d",fTitle.Data(),idim),nbins,xlow,xup);
	fHistos[idim]->SetDirectory(0);
    }
}

AliTRDNDFast::AliTRDNDFast(const AliTRDNDFast &ref) : TObject(ref),
fNDim(ref.fNDim),
fTitle(ref.fTitle),
fFunc(NULL),
fHistos(NULL),
fPars()
{
    Cleanup();
    Init();
    for(Int_t idim=0;idim<fNDim;idim++){
	fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim));
	fHistos[idim]->SetDirectory(0);
	for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
    }
}

AliTRDNDFast &AliTRDNDFast::operator=(const AliTRDNDFast &ref){
    //
    // Assignment operator
    //
    if(this != &ref){
	TObject::operator=(ref);
	fNDim=ref.fNDim;
	fTitle=ref.fTitle;
	fFunc = ref.fFunc;
	for(Int_t idim=0;idim<fNDim;idim++){
	  fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim));
	  fHistos[idim]->SetDirectory(0);
	  for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
	}
    }
    return *this;
}

AliTRDNDFast::~AliTRDNDFast(){
    Cleanup();

}

void AliTRDNDFast::Init(){

    for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].Set(fNDim);
    fFunc=new TF1*[fNDim];
    fHistos=new TH1F*[fNDim];
    for(Int_t idim=0;idim<fNDim;idim++){
        fHistos[idim]=NULL;
	for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(0,idim);
	fFunc[idim]=NULL;
    }
}

void AliTRDNDFast::Cleanup(){
    if(fHistos){
	for(Int_t idim=0;idim<fNDim;idim++){
	    if(fHistos[idim]){
		delete fHistos[idim];
		fHistos[idim]=NULL;
	    }
	}
	delete[] fHistos;
	fHistos=NULL;
    }
    if(fPars){
	for(Int_t ipar=0;ipar<kNpar;ipar++){
	    fPars[ipar].Reset();
	}
    }
    if(fFunc){
	for(Int_t idim=0;idim<fNDim;idim++){
	    if(fFunc[idim]){
		delete fFunc[idim];
		fFunc[idim]=NULL;
	    }
	}
	delete[] fFunc;
	fFunc=NULL;
    }
}

void AliTRDNDFast::PrintPars(){
    for(Int_t idim=0;idim<fNDim;idim++){
	for(Int_t ipar=0;ipar<kNpar;ipar++)cout<<idim<<" "<<ipar<<" "<<GetParam(idim,ipar)<<endl;
    }
}

void AliTRDNDFast::ScaleLangauFun(TF1 *func,Double_t mpv){

    Double_t start[kNpar],low[kNpar],high[kNpar];
    for(Int_t ii=0;ii<kNpar;ii++){
	start[ii]=func->GetParameter(ii);
	func->GetParLimits(ii,low[ii],high[ii]);
    }

    Double_t scalefactor=mpv/start[1];

    for(Int_t ii=0;ii<kNpar;ii++){
	Double_t scale=1;
	if(ii==0||ii==1||ii==3)scale=scalefactor;
	if(ii==4)scale=1./scalefactor;
	start[ii]*=scale;
	low[ii]*=scale;
	high[ii]*=scale;
	func->SetParLimits(ii, low[ii], high[ii]);
    }
    func->SetParameters(start);
}

//---------------------------------------------------------------
//---------------------------------------------------------------
TF1 *AliTRDNDFast::GetLangauFun(TString funcname,Double_t range[2],Double_t scalefactor){

    Double_t start[kNpar] = {120, 1000, 1, 100, 1e-5};
    Double_t low[kNpar] = {10, 300, 0.01, 1, 0.};
    Double_t high[kNpar] = {1000, 3000, 100, 1000, 1.};

    TF1 *hlandauE=new TF1(funcname.Data(),langaufunc,0,range[1],kNpar);
    hlandauE->SetParameters(start);
    hlandauE->SetParNames("Width","MP","Area","GSigma","delta");
    for (int i=0; i<kNpar; i++) {
	hlandauE->SetParLimits(i, low[i], high[i]);
    }

    if(scalefactor!=1){ScaleLangauFun(hlandauE,scalefactor*start[1]);}

    return hlandauE;
}

TF1 *AliTRDNDFast::FitLandau(TString name,TH1F *htemp,Double_t range[2],TString option){

    TF1 *fitlandau1D=GetLangauFun(name,range);
    TF1 fit("land","landau");
    Double_t max=htemp->GetXaxis()->GetBinCenter(htemp->GetMaximumBin());
    fit.SetParameter(1,max);
    fit.SetParLimits(1,0,htemp->GetXaxis()->GetXmax());
    fit.SetParameter(2,0.3*max); // MPV may never become negative!!!!!
    htemp->Fit("land",option.Data(),"",range[0],range[1]);
    ScaleLangauFun(fitlandau1D,fit.GetParameter(1));
    htemp->Fit(fitlandau1D,option.Data(),"",range[0],range[1]); // range for references
    return fitlandau1D;
}

void AliTRDNDFast::BuildHistos(){

    for(Int_t idim=0;idim<fNDim;idim++){
        fHistos[idim]->Reset();
	// Fill Histo
	Double_t pars[kNpar];
	for(Int_t ipar=0;ipar<kNpar;ipar++)pars[ipar]=GetParam(idim,ipar);
        // Also Fill overflow bins!!!
	for(Int_t ii=1;ii<=fHistos[idim]->GetNbinsX()+1;ii++){
	    Double_t xval=fHistos[idim]->GetXaxis()->GetBinCenter(ii);
	    Double_t val=langaufunc(&xval,&pars[0]);
	    //Double_t val=fFunc[idim]->Eval(xval);
	    fHistos[idim]->SetBinContent(ii,val);
	}
	// Normalization
	fHistos[idim]->Scale(1./fHistos[idim]->Integral(1,fHistos[idim]->GetNbinsX(),"width"));
    }
}

void AliTRDNDFast::Build(Double_t **pars){
    // pars[ndim][npar]
    for(Int_t idim=0;idim<fNDim;idim++){
	for(Int_t ipar=0;ipar<kNpar;ipar++){
	    fPars[ipar].SetAt(pars[idim][ipar],idim);
	}
    }
    BuildHistos();
}


void AliTRDNDFast::Build(TH1F **hdEdx,TString path){

    Double_t range[2];

    TCanvas *CANV=new TCanvas("a","a");
    if(fNDim==2)CANV->Divide(2,1);
    if(fNDim==3)CANV->Divide(2,2);
    if(fNDim>6)CANV->Divide(3,3);
    // Fit and Extract Parameters
    for(Int_t idim=0;idim<fNDim;idim++){
	CANV->cd(idim+1);
        gPad->SetLogy();
	range[0]=hdEdx[idim]->GetXaxis()->GetXmin();
	range[1]=hdEdx[idim]->GetXaxis()->GetXmax();
	// Norm Histogram
	hdEdx[idim]->Scale(1./hdEdx[idim]->Integral(1,hdEdx[idim]->GetNbinsX(),"width"));
        // Fit Histogram
	fFunc[idim]=FitLandau(Form("fit%d",idim),hdEdx[idim],range,"RMB0");
	// Norm Landau
	fFunc[idim]->SetParameter(2,fFunc[idim]->GetParameter(2)/fFunc[idim]->Integral(range[0],range[1]));
	hdEdx[idim]->DrawCopy();
        fFunc[idim]->DrawCopy("same");
	// Set Pars
	for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(fFunc[idim]->GetParameter(ipar),idim);
    }
    if(path!="")CANV->Print(Form("%s/%s_Build.pdf",path.Data(),fTitle.Data()));
    delete CANV;

    BuildHistos();
}

Double_t AliTRDNDFast::Eval(Double_t *point) const{
    Double_t val=1;
    for(Int_t idim=0;idim<fNDim;idim++){
	Int_t bin=fHistos[idim]->GetXaxis()->FindBin(point[idim]);
	val*=fHistos[idim]->GetBinContent(bin);
    }
    return val;
}

void AliTRDNDFast::Random(Double_t *point) const{
    for(Int_t idim=0;idim<fNDim;idim++){
	point[idim]=fHistos[idim]->GetRandom();
    }
}

void AliTRDNDFast::Random(Double_t *point,AliTRDNDFast *nd0,AliTRDNDFast *nd1,Double_t w0,Double_t w1){
    for(Int_t idim=0;idim<nd0->fNDim;idim++){
	point[idim]=GetRandomInterpolation(nd0->fHistos[idim],nd1->fHistos[idim],w0,w1);
    }
}

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){

    if((end-start)==1)return start;
    Int_t mid=Int_t((end+start)/2);
    Double_t valmid=(w0*a0[mid]+w1*a1[mid])/(w0+w1);
    if(val>=valmid)return BinarySearchInterpolation(mid,end,a0,a1,w0,w1,val);
    return BinarySearchInterpolation(start,mid,a0,a1,w0,w1,val);
}

Double_t AliTRDNDFast::GetRandomInterpolation(TH1F *hist0,TH1F *hist1,Double_t w0,Double_t w1){

    // Draw Random Variable
    Double_t rand=gRandom->Rndm();

    // Get Integral Arrays
    Double_t *integral0=hist0->GetIntegral();
    Double_t *integral1=hist1->GetIntegral();

    Int_t nbinsX=hist0->GetNbinsX();

    Int_t ibin=BinarySearchInterpolation(1,nbinsX+1,integral0,integral1,w0,w1,rand);
    return hist0->GetXaxis()->GetBinCenter(ibin);
}
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=(w0a0[mid]+w1a1[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