#include <iostream>
#include <fstream>
#include "TH1D.h"
#include "TRandom.h"
#include "TFile.h"
#include "TTree.h"
#include "TVector3.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TF1.h"
#include "TH1D.h"
#include "TCanvas.h"
using namespace std;

double cxb_HE_spc(double x1, double x2){
	double ray1 = 0.54/0.4*(pow(0.015,-0.4)-pow(0.02,-0.4)); 
	double ray2 = 0.0117/1.38*(pow(0.02,-1.38)-pow(0.1,-1.38));
	double ray3 = 0.014/1.3*(pow(0.1,-1.3)-pow(100.,-1.3));
	double ray  = ray1+ray2+ray3;
	double WhichOne = x1;
	double Rand_energy = x2;
	double energy_output;
	if(WhichOne <= ray1/ray)
		energy_output = pow(Rand_energy *pow(0.02,-0.4)+(1.-Rand_energy )*pow(0.015,-0.4) ,-1./0.4);
	else if( ray1/ray< WhichOne && WhichOne <=(ray1+ray2)/ray)
		energy_output = pow(Rand_energy *pow(0.1,-1.38)+(1.-Rand_energy )*pow(0.02,-1.38) ,-1./1.38);
	else
		energy_output = pow(Rand_energy *pow(100.,-1.3)+(1.-Rand_energy )*pow(0.1,-1.3) ,-1./1.3);
	return energy_output;
}

double cxb_HE_spc_func(double *x, double *par){
	double f=0.;
	if(x[0]<0.02){
		f = 0.54*pow(x[0], -1.4);
	}
	else if(x[0]>=0.02&&x[0]<0.1){
		f = 0.0117*pow(x[0], -2.38);
	}
	else{
		f = 0.014*pow(x[0], -2.3);
	}
	return f;
}


int main()
{

	int Nsample = 1000000;
	int seed = 1;
	TString outfilename = "test.root";

	delete gRandom;
	gRandom = new TRandom3(seed);

	TH1D *h1_sample = new TH1D("h1_sample", "h1_sample", 2000, 0.015, 0.1);
	h1_sample->SetDirectory(NULL);
	TH1D *h1_func   = new TH1D("h1_func", "h1_func", 2000, 0.015, 0.1);
	h1_func->SetDirectory(NULL);
	TH1D *h1_func2   = new TH1D("h1_func2", "h1_func2", 2000, 0.015, 0.1);
	h1_func2->SetDirectory(NULL);

	TF1* fun = new TF1("fun",cxb_HE_spc_func,0.015,100.);
	fun->SetNpx(10000000);
	//fun->SetNpx( h1_func->GetNbinsX());
	h1_func->FillRandom("fun",1000000);


	TFile* outfile = new TFile(outfilename,"recreate");
	for(int i = 0; i < Nsample; i++){

		h1_func2->Fill(fun->GetRandom());

		double energy = cxb_HE_spc(gRandom->Rndm(), gRandom->Rndm());
		h1_sample->Fill(energy);
	}
	TCanvas* c = new TCanvas();
	h1_func->SetLineColor(kBlue);
	h1_sample->SetLineColor(kRed);
	h1_func2->SetLineColor(kGreen);
	h1_func->Draw();
	h1_func2->Draw("same");
	h1_sample->Draw("same");
	c->Write();
	h1_sample->Write();
	h1_func->Write();
	h1_func2->Write();
	outfile->Close();
	delete outfile;
	outfile = NULL;

	cout<< " run successfully " << endl;
	return 0;


}