Generate random value following a hard PDF

Fynal · January 14, 2021, 9:27pm

Hi everyone.

For a little data analysis I need to use a hard function as a PDF to generate random values. So I’ve created the function and would build an histogram with the random calls to reproduce the curve. I’ll join the code at the end. The problem is that, when I call this TF1 function I’ve created, I have some out of range and integral errors (exp > 709, integral = 0 …) .

I used many cout tests to see what happened and it looks like, when I call only once the function (with GetRandom), the function fPLP (function power law + peak) doesn’t run only once, but many times till the crash. For exemple I got this :

Processing distrib.C…
0
t0
t0
t0
t0
t0
t0
t0
t0
t0
t0
Error: exp param[0]=1022.15 up:709 low:-inf out of range distrib.C:15:
t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t2t0t2t0t2t0t2t0t2t0t2t0t2t0t2t0t2t0t2t0t2t0t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2Error in TF1::GetRandom: Integral of function is zero
*** Interpreter error recovered ***
t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2t2root [1]

If I modify the limits (5, 95) the value in the error (1022.15) also changes. I do not understand how the return doesn’t ‘break’ the call and so why it loop over it keeping in mind cumulative values apparently till the crash .

Here is the code :
distrib.C (2.5 KB)

Do someone have any solutions there ?
Thanks a lot.

Wile_E_Coyote · January 14, 2021, 10:11pm

First, fix all errors and warnings reported by:
root [0] .L distrib.C++

Fynal · January 15, 2021, 10:45am

Thanks for your answer Wile_E ,

Those Warning: Illegal numerical expression .L fPLP.C comes from all the returns ? Because I can’t see errors in the returns.
If I comment the .L fPLP.C ; I have no more warnings but still the same strange result (t0
t0
Error: exp param[0]=1022.1 up:709 low:-inf out of range distrib.C:18:
t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t0t …)

Wile_E_Coyote · January 15, 2021, 10:57am

root -l -q -e '.L distrib.C++'

Fynal · January 15, 2021, 11:05am

Warning in TApplication::GetOptions: macro .L distrib.C not found

Fynal · January 15, 2021, 11:06am

I’m not in local but on a server via ssh I think I can’t load any function maybe

Wile_E_Coyote · January 15, 2021, 11:18am

ROOT Manual -> Basics -> First steps with ROOT
ROOT Primer -> ROOT Macros

Fynal · January 15, 2021, 11:26am

Well, I created another .C containing only my fPLP fonction and made .L fPLP.C, no errors, then the same with .L distrib.C (removing in it the fPLP func), no errors, and when I run distrib.C (which calls fPLP) I have the same number of test (t0 t2 …) printed.

Wile_E_Coyote · January 15, 2021, 11:27am

You need “++”, e.g.:
.L distrib.C++

Fynal · January 15, 2021, 12:02pm

Ok so, here are the commands I’ve done :

root [0] .L fPLP.C++
Info in TUnixSystem::ACLiC: creating shared library /data/dev/cbc1/andres/ada/./fPLP_C.so
In file included from /data/dev/cbc1/andres/ada/fPLP_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/fPLP_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./fPLP.C: In function ‘Double_t fPLP(Double_t*, Double_t*)’:
/data/dev/cbc1/andres/ada/./fPLP.C:33:1: warning: control reaches end of non-void function [-Wreturn-type]
}
^
root [1] .L distrib.C++
Info in TUnixSystem::ACLiC: creating shared library /data/dev/cbc1/andres/ada/./distrib_C.so
In file included from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./distrib.C: In function ‘void distrib()’:
/data/dev/cbc1/andres/ada/./distrib.C:63:46: error: invalid use of incomplete type ‘class TF1’
TF1 fp = new TF1(“fct plp”,fPLP,4.53,86.73,7);
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:2:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/virgoApp/root/v5r34p050/Linux-x86_64-CL7/include/TH1.h:66:7: error: forward declaration of ‘class TF1’
class TF1;
^
In file included from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./distrib.C:64:3: error: invalid use of incomplete type ‘class TF1’
fp->SetParNames(“Mmin”,“Mmax”,“Alpha”,“Mu”,“Sigma”,“DeltaM”,“LambdaP”);
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:2:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/virgoApp/root/v5r34p050/Linux-x86_64-CL7/include/TH1.h:66:7: error: forward declaration of ‘class TF1’
class TF1;
^
In file included from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./distrib.C:65:3: error: invalid use of incomplete type ‘class TF1’
fp->SetParameters(4.53,86.73,2.62,33.5,5.09,4.88,0.08);
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:2:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/virgoApp/root/v5r34p050/Linux-x86_64-CL7/include/TH1.h:66:7: error: forward declaration of ‘class TF1’
class TF1;
^
In file included from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./distrib.C:71:17: error: invalid use of incomplete type ‘class TF1’
Double_t x = fp->GetRandom(4.5,90.); // mon fit fplp
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:2:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/virgoApp/root/v5r34p050/Linux-x86_64-CL7/include/TH1.h:66:7: error: forward declaration of ‘class TF1’
class TF1;
^
In file included from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./distrib.C:81:3: error: invalid use of incomplete type ‘class TF1’
fp->Draw(“same”);
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:2:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/virgoApp/root/v5r34p050/Linux-x86_64-CL7/include/TH1.h:66:7: error: forward declaration of ‘class TF1’
class TF1;
^
In file included from /data/dev/cbc1/andres/ada/./distrib.C:14:0,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.h:34,
from /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.cxx:17:
/data/dev/cbc1/andres/ada/./fPLP.C: In function ‘Double_t fPLP(Double_t, Double_t*)’:
/data/dev/cbc1/andres/ada/./fPLP.C:33:1: warning: control reaches end of non-void function [-Wreturn-type]
}
^
g++: error: /data/dev/cbc1/andres/ada/distrib_C_ACLiC_dict.o: No such file or directory
Error in : Compilation failed!
root [2] .L distrib.C
root [3] .x distrib.C
root [4]

Still remains many errors apparently but I obtain my histogram, it seems to work. I’m confused

Wile_E_Coyote · January 15, 2021, 12:10pm

For the “control reaches end of non-void function”, you need to fix this function so that it always returns a valid value.
For the “incomplete type 'class TF1'”, you obviously need: #include "TF1.h"

Fynal · January 15, 2021, 12:14pm

Yes I forget the TF1 thanks … For the fPLP.C function I have :

#include <TRandom.h>
#include
#define M_PI 3.14159265358979323846
Double_t fPLP(Double_t *x, Double_t *par)
{
Double_t xx=x[0];
if (xx < par[0])
{
return 0;
}
else if (xx >= par[0] && xx < par[0] + par[5])
{
if (xx <= par[1])
{
return ((1-par[6])pow(xx,-par[2])+par[6](1/(par[4]sqrt(2M_PI))exp(-pow((xx-par[3]),2)/(2par[4]*par[4]))))pow(exp(par[5]/(xx-par[0])+par[5]/(xx-par[0]-par[5]))+1,-1);
}
else
{
return par[6](1/(par[4]sqrt(2M_PI))exp(-pow((xx-par[3]),2)/(2par[4]*par[4])))*pow(exp(par[5]/(xx-par[0])+par[5]/(xx-par[0]-par[5]))+1,-1);
}
}
else if (x[0] >= par[0] + par[5])
{
if (x[0] <= par[1])
{
return (1-par[6])pow(xx,-par[2])+par[6](1/(par[4]sqrt(2M_PI))exp(-pow((xx-par[3]),2)/(2par[4]par[4])));
}
else
{
return par[6](1/(par[4]sqrt(2M_PI))exp(-pow((xx-par[3]),2)/(2par[4]*par[4])));
}
}
}

Why there should be invalid values if I don’t have any errors ?

Wile_E_Coyote · January 15, 2021, 12:17pm

For example, you do not return any value when none of the “if” statements fired (the compiler is not really able to verify that something like this may never happen).

distrib.C (2.9 KB)

Fynal · January 15, 2021, 12:25pm

Well I verified , every value in the generated interval must verify at least 1 statement (and values out of the interval will return the first : 0 or the last return), anyway, it seems to work. Maybe I have to declare this function as void ? (control reaches end of non-void function [-Wreturn-type])

system · January 29, 2021, 12:25pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.