divaldo
/
BM1D
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

jump

This commit is contained in:
Szántó Géza 2017-11-03 15:45:39 +01:00
parent 3b9b3ba2e8
commit ae3084b1fd
6 changed files with 232 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
BM1D.cc
View File

@ -12,47 +12,30 @@ int main(int argc, char* argv[])
Int_t nSteps, nRuns;
Double_t p0,p1,x1,x2;
Double_t mu1, mu2, sigma1, sigma2;
Double_t j_mu1, j_sigma1;//////////////////////
Double_t rat;//////////////////////
const char* fileName="input.root";
char random_type='u';
char random_type='j';
Int_t vis, typeOfRun;
nSteps=nRuns=vis=typeOfRun=0;
mu1=mu2=sigma1=sigma2=p0=p1=x1=x2=0.0;
if(argc==15)
{
nSteps=atoi(argv[1]);
nRuns=atoi(argv[2]);
p0=atof(argv[3]);
p1=atof(argv[4]);
x1=atof(argv[5]);
x2=atof(argv[6]);
mu1=atof(argv[7]);
mu2=atof(argv[8]);
sigma1=atof(argv[9]);
sigma2=atof(argv[10]);
fileName=argv[11];
random_type=argv[12][0];
vis=atoi(argv[13]); if((vis!=0)&&(vis!=1)){vis=0;}
typeOfRun=atoi(argv[14]); if((typeOfRun!=0)&&(typeOfRun!=1)&&(typeOfRun!=2)){typeOfRun=0;}
}
else
{
//default runs with less parameters
random_type = 'g';
random_type = 'j';
nSteps = 1000;
nRuns = 1;
p0 = 0.5;
nRuns = 5;
p0 = 0;
p1 = 0;
x1 = 0;
x2 = 0;
x1 = 200000;
x2 = 200000;
mu1 = 1;
mu2 = 0;
sigma1 = 2;
sigma2 = 0;
mu2 = -6;
sigma1 = 0.8;
sigma2 = 0.2;
vis = 1;
typeOfRun = 1;
}
j_mu1=-25;
j_sigma1=0.3;
rat=0.9;
TApplication App("tapp", &argc, argv);
BM1DProcess *myBM1DProcess = new BM1DProcess();
@ -66,22 +49,43 @@ int main(int argc, char* argv[])
case 'l':
myBM1DProcess->Run(nRuns, nSteps, p0, x1, x2, mu1, sigma1, mu2, sigma2);
break;
case 'j':
myBM1DProcess->Run(nRuns, nSteps, p0, x1, x2, mu1, sigma1, mu2, sigma2, j_mu1, j_sigma1, rat);
//cout << "\nOK\n";
break;
default:
cout<<"ERROR! Wrong parameter for type of random generator! \n No run!"<<endl;
break;
}
TH1D *hJump = new TH1D("h", "test", 100, -30., 30.);
hJump->Sumw2();
std::vector<Double_t> x = myBM1DProcess->GetX();
hJump->Fill(0);
for(int i=1; i<x.size(); i++){
hJump->Fill(x[i]-x[i-1]);
}
hJump->Draw();
// std::cout << "Pontok szama: "<< x.size() << "\n";
Plotter* myPlotter = new Plotter(vis==1);
myPlotter->Plot(nRuns, nSteps, myBM1DProcess->GetT(), myBM1DProcess->GetX());
Analyse *myAnalyse = new Analyse();
switch(random_type){
case 'g' :
myAnalyse->AnalyseGaus(myBM1DProcess->GetT(),myBM1DProcess->GetX());
}
BM1DSave *save = new BM1DSave();
save->SaveToTree(myPlotter->GetTmultiGraph(), p0, p1, nSteps, nRuns, x1, x2, mu1, mu2, sigma1, sigma2, myBM1DProcess->GetT(), myBM1DProcess->GetX(),fileName);
// Analyse *myAnalyse = new Analyse();
// switch(random_type){
// case 'g' :
// myAnalyse->AnalyseGaus(myBM1DProcess->GetT(),myBM1DProcess->GetX());
// }
// BM1DSave *save = new BM1DSave();
// save->SaveToTree(myPlotter->GetTmultiGraph(), p0, p1, nSteps, nRuns, x1, x2, mu1, mu2, sigma1, sigma2, myBM1DProcess->GetT(), myBM1DProcess->GetX(),fileName);
App.Run();
return 0;
}

87
BM1D.cc_ Normal file
View File

@ -0,0 +1,87 @@
#include <iostream>
#include "BM1DProcess.hh"
#include "Plotter.hh"
#include "TApplication.h"
#include "Analyse.hh"
#include "BM1DSave.hh"
using namespace std;
int main(int argc, char* argv[])
{
Int_t nSteps, nRuns;
Double_t p0,p1,x1,x2;
Double_t mu1, mu2, sigma1, sigma2;
const char* fileName="input.root";
char random_type='u';
Int_t vis, typeOfRun;
nSteps=nRuns=vis=typeOfRun=0;
mu1=mu2=sigma1=sigma2=p0=p1=x1=x2=0.0;
if(argc==15)
{
nSteps=atoi(argv[1]);
nRuns=atoi(argv[2]);
p0=atof(argv[3]);
p1=atof(argv[4]);
x1=atof(argv[5]);
x2=atof(argv[6]);
mu1=atof(argv[7]);
mu2=atof(argv[8]);
sigma1=atof(argv[9]);
sigma2=atof(argv[10]);
fileName=argv[11];
random_type=argv[12][0];
vis=atoi(argv[13]); if((vis!=0)&&(vis!=1)){vis=0;}
typeOfRun=atoi(argv[14]); if((typeOfRun!=0)&&(typeOfRun!=1)&&(typeOfRun!=2)){typeOfRun=0;}
}
else
{
//default runs with less parameters
random_type = 'g';
nSteps = 1000;
nRuns = 1;
p0 = 0.5;
p1 = 0;
x1 = 0;
x2 = 0;
mu1 = 1;
mu2 = 0;
sigma1 = 2;
sigma2 = 0;
vis = 1;
typeOfRun = 1;
}
TApplication App("tapp", &argc, argv);
BM1DProcess *myBM1DProcess = new BM1DProcess();
switch(random_type){
case 'u':
myBM1DProcess->Run(nRuns, nSteps, p0, p1);
break;
case 'g':
myBM1DProcess->Run(nRuns, nSteps, p0, mu1, sigma1);
break;
case 'l':
myBM1DProcess->Run(nRuns, nSteps, p0, x1, x2, mu1, sigma1, mu2, sigma2);
break;
default:
cout<<"ERROR! Wrong parameter for type of random generator! \n No run!"<<endl;
break;
}
Plotter* myPlotter = new Plotter(vis==1);
myPlotter->Plot(nRuns, nSteps, myBM1DProcess->GetT(), myBM1DProcess->GetX());
Analyse *myAnalyse = new Analyse();
switch(random_type){
case 'g' :
myAnalyse->AnalyseGaus(myBM1DProcess->GetT(),myBM1DProcess->GetX());
}
BM1DSave *save = new BM1DSave();
save->SaveToTree(myPlotter->GetTmultiGraph(), p0, p1, nSteps, nRuns, x1, x2, mu1, mu2, sigma1, sigma2, myBM1DProcess->GetT(), myBM1DProcess->GetX(),fileName);
App.Run();
return 0;
}

4
include/BM1DProcess.hh
View File

@ -7,6 +7,8 @@
#include <cmath>
#include "TRandom.h"
#include "TF1.h"
#include "BM1DRandomGenerator.hh"
class BM1DProcess {
public:
@ -16,6 +18,7 @@ public:
void Run(int nSteps, int nRuns, double p0, double p1);
void Run(int nSteps, int nRuns, double p0, double mu, double sigma);
void Run(int nSteps, int nRuns, double p0, double x1, double x2, double mu1, double sigma1, double mu2, double sigma2);
void Run(int nRuns, int nSteps, double p0, double x1, double x2, double mu1, double sigma1, double mu2, double sigma2, double j_mu1, double j_sigma1, double rat);
std::vector<Double_t> GetT(){return t;}
std::vector<Double_t> GetX(){return x;}
@ -23,6 +26,7 @@ public:
private:
TRandom* randomGenerator;
TRandom* randomGeneratorGauss;
//TF1* jump;
Double_t rand1;
std::vector<Double_t> t;
std::vector<Double_t> x;

23
include/BM1DRandomGenerator.hh Normal file
View File

@ -0,0 +1,23 @@
#ifndef BM1DRandomGenerator_h
#define BM1DRandomGenerator_h 1
#endif
#include "TRandom3.h"
#include "TF1.h"
class BM1DRandomGenerator: public TRandom3
{
public:
BM1DRandomGenerator();
// ~BM1DRandomGenerator();
//void SetTF1();
void SetTF1();
void CallTF1();
void SetJump(Double_t min_, Double_t max_); //set limits
Double_t Jump(Double_t mu, Double_t sigma, Double_t mu_jump, Double_t sigma_jump, Double_t rat);
//Double_t Cauchy();
private:
TF1* tf1;
};

54
src/BM1DProcess.cc
View File

@ -4,12 +4,15 @@ BM1DProcess::BM1DProcess()
{
randomGenerator = new TRandom();
randomGeneratorGauss = new TRandom();
//jump = new TF1("jump","gaus(0)+gaus(3)",-30,30);
}
BM1DProcess::~BM1DProcess()
{
delete randomGenerator;
delete randomGeneratorGauss;
// delete jump;
}
@ -120,3 +123,54 @@ void BM1DProcess::Run(int nRuns, int nSteps, double p0, double x1, double x2, do
}
}
}
void BM1DProcess::Run(int nRuns, int nSteps, double p0, double x1, double x2, double mu1, double sigma1, double mu2, double sigma2, double j_mu1, double j_sigma1, double rat)
{
t.clear();
x.clear();
BM1DRandomGenerator* jump = new BM1DRandomGenerator();
Double_t minim = (mu1-5*sigma1 < j_mu1-5*j_sigma1) ? mu1-5*sigma1 : j_mu1-5*j_sigma1;
Double_t maxim = (mu1+5*sigma1 > j_mu1+5*j_sigma1) ? mu1+5*sigma1 : j_mu1+5*j_sigma1;
jump->SetJump(minim, maxim);
for(int i = 0; i < nRuns; i++) //multiple runs
{
t.push_back(0.0); //let's start at t=0, x=0, you can change it if you vant, please use Set methods
x.push_back(0.0);
double randJ;
for(int ii = 1; ii < nSteps ; ii++)
{
rand1 = randomGenerator->Uniform();
if(rand1 < p0) //step in time, but no step in x
{
t.push_back(t.back() + 1);
x.push_back(x.back());
}
else //step left or right
{
if(x.back() < x2 && x.back() > x1)
{
// par[1] = mu2; par[2] = sigma2;
// jump->SetParameters(par);
// randJ = jump->GetRandom();
randJ=jump->Jump(mu2, sigma2, j_mu1, j_sigma1, rat );
}
else
{
// par[1] = mu1; par[2] = sigma1;
// jump->SetParameters(par);
// randJ = jump->GetRandom();
randJ=jump->Jump(mu1, sigma1, j_mu1, j_sigma1, rat );
}
t.push_back(t.back() + 1);
x.push_back(x.back() + randJ); //one step right
}
}
}
}

18
src/BM1DRandomGenerator.cc Normal file
View File

@ -0,0 +1,18 @@
#include "BM1DRandomGenerator.hh"
BM1DRandomGenerator::BM1DRandomGenerator(): TRandom3(){
tf1 = new TF1("jump","gaus(0)+gaus(3)",-30,30);
}
void BM1DRandomGenerator::SetJump(Double_t min_, Double_t max_){ //set limits
tf1 = new TF1("jump","gaus(0)+gaus(3)", min_, max_);
}
Double_t BM1DRandomGenerator::Jump(Double_t mu, Double_t sigma, Double_t mu_jump, Double_t sigma_jump, Double_t rat){
Double_t par[6];
par[0] = rat; par[1]=mu; par[2]=sigma; par[3]=1-rat; par[4]=mu_jump; par[5]=sigma_jump;
tf1->SetParameters(par);
return tf1->GetRandom();
}