Dear ROOTer,
At first, I must say that I am using root v5.22 on a Mac OS X Tiger 10.5.6
I am trying to randomly filled a 2 dimensional histogram from a personalized TF2 function [Convolution of a Gaussian(X,Y) with exp(-|x|/h-|y|/h’)].
Here, is my programs, embedded into an header file:
[code]//*************************************************************
// FUNCTION TO CREATE HISTOGRAMM RANDOMLY FILLED BY The CONFOLUTION OF DECAYING EXPONANTIAL WITH GAUSSIAN
//*************************************************************
TH2D *FillExpG2D(TH2 *IN, double Xmean, double XSig0, double XSig1, double Ymean, double YSig0, double YSig1)
{
// IN -> the histogram to be filled (LIMIT MUST BE [0,360] or [-180,180] for the X axis and [-90,90] for the Y axis).
// Xmean -> the averaged X position
// XSig0 -> Sigma of the gaussian
// XSig1 -> Scale height of the decaying exponential exp(-0.5|x|/XSig1)
// IDEM FOR Ymean, YSig0, YSig1
bool IsCentered0 = false;
if(IN->GetXaxis()->GetXmin() < 0) IsCentered0 = true;
double Xmax, Xmin, NXbin;
double Ymax, Ymin, NYbin;
double X,Y;
Xmax = IN->GetXaxis()->GetXmax();
Xmin = IN->GetXaxis()->GetXmin();
NXbin= IN->GetNbinsX();
Ymax = IN->GetYaxis()->GetXmax();
Ymin = IN->GetYaxis()->GetXmin();
NYbin= IN->GetNbinsY();
TF2 *Exp2 = new TF2("Exp2",ExpGausConvolution,
Xmin-10*(XSig0+XSig1),Xmax+10*(XSig0+XSig1),Ymin-10*(YSig0+YSig1),Ymax+10*(YSig0+YSig1), 7);
//// IN CASE OF THE 2D Gaussian :
//// TF2 Exp2 = new TF2(“G2”,“xygaus”,Xmin-10XSig,Xmax+10XSig,Ymin-10YSig,Ymax+10YSig);
Exp2->SetNpx(int((Xmax-Xmin+2(10*(XSig0+XSig1)))/(Xmax-Xmin)NXbin));
Exp2->SetNpy(int((Ymax-Ymin+2(10*(YSig0+YSig1)))/(Ymax-Ymin)*NYbin));
Exp2->SetParameters(1.,Xmean,XSig0, XSig1,Ymean,YSig0,YSig1);
//// IN CASE OF THE 2D Gaussian :
//// Exp2->SetParameters(1.,Xmean,XSig0,Ymean,YSig0);
TH2D *TEMP= (TH2D*) IN->Clone("TEMP");
TEMP->Scale(0);
TEMP->SetEntries(0);
for(int i=1; i<= 1000; i++)
{
Exp2->GetRandom2(X,Y);
if(!IsCentered0)
{
if(fabs(Y) > 90) {Y = 180.-Y; X+=IN->GetXaxis()->GetXmax()/2.;}
while(X > IN->GetXaxis()->GetXmax()) X -=IN->GetXaxis()->GetXmax();
while(X < 0) X +=IN->GetXaxis()->GetXmax();
}
else
{
if(fabs(Y) > 90) {Y = 180.-Y; X+=IN->GetXaxis()->GetXmax();}
while(X > IN->GetXaxis()->GetXmax()) X =IN->GetXaxis()->GetXmin() + (IN->GetXaxis()->GetXmin() + X);
while(X < IN->GetXaxis()->GetXmin()) X =IN->GetXaxis()->GetXmax() + (IN->GetXaxis()->GetXmax() + X);
}
TEMP->Fill(X,Y);
}
TEMP->Scale(1./TEMP->Integral("W"));
TEMP->SetEntries(1);
delete Exp2;
return TEMP;
}
// AN APPROXIMATION OF THE ERROR FUNCTION
double ERF(double x)
{
double Pi=acos(-1);
double a=-8./(3.Pi)(Pi-3.)/(Pi-4.);
double erf = 1. - exp(-pow(x,2.)(4./Pi+apow(x,2.))/(1+apow(x,2.)));
if(fabs(x) > 0) return fabs(x)/xsqrt(erf);
else return sqrt(erf);
}
// THE CONVOLUTION OF THE GAUSSIAN WITH exp(-|x|/h-|y|/h’)
double ExpGausConvolution(double *x, double *par)
{
double ConvX, ConvY;
double Norm=par[0], X0=par[1], Sx=par[2], Hx=par[3],
Y0=par[4], Sy=par[5], Hy=par[6];
double X=x[0]-X0, Y=x[1]-Y0;
ConvX = exp((4*Hx*X+pow(Sx,2.))/(8*pow(Hx,2)))*(1-ERF(sqrt(2)/4.*(2*Hx*X+pow(Sx,2.))/(Sx*Hx)))
+ exp(-(4*Hx*X-pow(Sx,2.))/(8*pow(Hx,2)))*(1+ERF(sqrt(2)/4.*(2*Hx*X-pow(Sx,2.))/(Sx*Hx)));
ConvY = exp((4*Hy*Y+pow(Sy,2.))/(8*pow(Hy,2)))*(1-ERF(sqrt(2)/4.*(2*Hy*Y+pow(Sy,2.))/(Sy*Hy)))
+ exp(-(4*Hy*Y-pow(Sy,2.))/(8*pow(Hy,2)))*(1+ERF(sqrt(2)/4.*(2*Hy*Y-pow(Sy,2.))/(Sy*Hy)));
if(ConvX*ConvY > 0) return Norm*ConvX*ConvY;
else return 0;
}
[/code]
when I run it,
#include<MyHeader.h>
TH2D *IN = new TH2D("TETS","TEST",183,0,360,82,-90,90);
TH2D *OUT=(TH2D*) FillExpG2D(IN,130,1.,1.5,-10,1,1.5);it crashes while it try to get the random number at the line Exp2->GetRandom2(X,Y). The error message is quite long but it start as follow:
*** Break *** bus error
/Volumes/USERS-DATA/gillard/PAMELA/CR-MAP/ELECTRON/DAYLY/32103: No such file or directory.
Attaching to process 32103.
Reading symbols for shared libraries . done
Reading symbols for shared libraries ............. done
0x906c0189 in wait4 ()
...
It seems the problem occur when at the Exp2->GetRandom(X,Y) function and, the Exp2 TF2 object seems to be responsible of the problem. If now I change my personalized function by a simple “gaussxy”, the program run fine.
However, if I create the Exp2 function from the ROOT shell directly (see code below) and try to plot it, or obtain 1000 random number X and Y, I do not have problems.
double X, Y;
TF2 *Exp2 = new TF2("Exp2",ExpGausConvolution,0,360,-90,90,7);
Exp2->SetParameters(1.,180,1, 5.5,0,1,5.5);
Exp2->Draw("COLZ");
for(int i=0;i<1000;i++) {Exp2->GetRandom2(X,Y); cout<<" "<<X<<" "<<Y<<endl;}So, if someone can help me to solve my problem, I will be very grateful.
).