目录
1.问题描述
2.解决方法
3.代码实现
1.问题描述
更具自定义的数据点,生成能量的随机抽样;
需要的数据分布点;
2.解决方法
先积分然后归一化得到概率函数,采用二分法找到对应的区间的X;
生成的随机抽样分布如上;
3. 代码实现
#include "iostream"
using namespace std;
void GenerateTest()
{
vector<double> bins,vals,nvals;
//读取数据
double x=0.,y=0.;
ifstream input("Alt100.txt");
bins.clear();vals.clear();
bins.push_back(0);
vals.push_back(0);
for(int i=1;i<191;i++)
{
input>>x>>y;
bins.push_back(x);
vals.push_back(y);
}
//求和和归一化
double sum=0;
vals.clear();
for(int ii=1;ii<191;ii++)
{
// value=vals[ii]*(bins[ii]-bins[ii-1]);
// sum = sum+value;
vals[ii]=vals[ii];//*(bins[ii]-bins[ii-1]);
sum = sum+vals[ii];
// cout<<" Sum"<<sum<<endl;
// cout<<""<<vals[ii]<<endl;
}
// for(int ii=1;ii<190;ii++)
// {
// vals[ii]=vals[ii]/(bins[ii]-bins[ii-1]);
// // sum = sum+vals[ii];
// // cout<<" Sum"<<sum<<endl;
// // cout<<""<<vals[ii]<<endl;
// }
nvals.clear();
//
double norm=0.;
for(int iii=1;iii<191;iii++)
{
norm = vals[iii]/sum+norm;
cout<<" norm "<<norm<<endl;
// cout<<" vals"<<vals[iii]<<endl;
nvals.push_back(norm);
}
// for(int j=0;j<189;j++) cout<<nvals[j]<<endl;//test nvals
//
// double Emin=1.0e-8;
// double Emax=1000;
TRandom3 rndm;
//--------对数坐标变化--------
double low=1.e-8;
double up=1e4;
double nbin=192;
Double_t xAxis[193];
for ( int i=0; i <= nbin; i++) {
double val_bin=low * std::pow(10., i * std::log10(up/low)/nbin);
double exp_10=4.-int(std::log10(val_bin));
double factor =std::pow(10., exp_10);
val_bin=int(factor*val_bin)/factor;
xAxis[i] = val_bin;
}
TH1D* h1 = new TH1D("h1","h1",192,xAxis);
// TH1D* h1 = new TH1D("h1","h1",400,0.0,1.0);
//--------------------------
for(int j=0;j<1000000;j++){
int numberOfBin=nvals.size();
int nmin=0;
int nmid=numberOfBin/2;
int nmax=numberOfBin-1;
double a=rndm.Uniform();
while(nmin!=nmax-1)
{
if(a>nvals[nmid])
nmin=nmid;
else
nmax=nmid;
nmid=nmin+(nmax-nmin+1)/2;
}
double weight=(bins[nmid]-bins[nmid-1])/(nvals[nmid]-nvals[nmid-1]);
double rand=bins[nmid-1]+weight*(a-nvals[nmid-1]);
h1->Fill(rand);
}
TCanvas *c1 = new TCanvas();
c1->cd();
TGraph *gr=new TGraph("./Alt100.txt","%lg%lg");
gr->Draw();
TCanvas *c2 = new TCanvas();
c2->cd();
//
h1->Draw();
}