diff --git a/BM1D.cc b/BM1D.cc index 7ae7687..93d0813 100644 --- a/BM1D.cc +++ b/BM1D.cc @@ -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; - const char* fileName="input.root"; - char random_type='u'; + Double_t j_mu1, j_sigma1;////////////////////// + Double_t rat;////////////////////// + const char* fileName="input.root"; + 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!"<Sumw2(); + std::vector x = myBM1DProcess->GetX(); + hJump->Fill(0); + for(int i=1; iFill(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(); - 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); + // 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; } diff --git a/BM1D.cc_ b/BM1D.cc_ new file mode 100644 index 0000000..7ae7687 --- /dev/null +++ b/BM1D.cc_ @@ -0,0 +1,87 @@ +#include +#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!"<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; +} diff --git a/include/BM1DProcess.hh b/include/BM1DProcess.hh index 6b52b96..ee6c71c 100644 --- a/include/BM1DProcess.hh +++ b/include/BM1DProcess.hh @@ -7,6 +7,8 @@ #include #include "TRandom.h" +#include "TF1.h" +#include "BM1DRandomGenerator.hh" class BM1DProcess { public: @@ -16,13 +18,15 @@ 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 GetT(){return t;} std::vector GetX(){return x;} private: TRandom* randomGenerator; TRandom* randomGeneratorGauss; + //TF1* jump; Double_t rand1; std::vector t; std::vector x; diff --git a/include/BM1DRandomGenerator.hh b/include/BM1DRandomGenerator.hh new file mode 100644 index 0000000..4caa4b4 --- /dev/null +++ b/include/BM1DRandomGenerator.hh @@ -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; +}; diff --git a/src/BM1DProcess.cc b/src/BM1DProcess.cc index 2e5e9ca..6f00a7e 100644 --- a/src/BM1DProcess.cc +++ b/src/BM1DProcess.cc @@ -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; } @@ -119,4 +122,55 @@ void BM1DProcess::Run(int nRuns, int nSteps, double p0, double x1, double x2, do } } } -} \ No newline at end of file +} + + +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 + } + } + } +} diff --git a/src/BM1DRandomGenerator.cc b/src/BM1DRandomGenerator.cc new file mode 100644 index 0000000..0fbcf33 --- /dev/null +++ b/src/BM1DRandomGenerator.cc @@ -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(); +}