Skip to content

Commit

Permalink
Merge pull request #2 from baltzell/Project
Browse files Browse the repository at this point in the history 
Fix ambiguous namespace/local reference
  • Loading branch information
@maksaska
maksaska authored Jan 20, 2023
2 parents 62c0591 + c5bc2ac commit e812c6f
Show file tree
Hide file tree
Showing 5 changed files with 70 additions and 81 deletions.
20 changes: 9 additions & 11 deletions event.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,7 @@
#include <deque>
#include <fstream>

using namespace std;

vector<vector<double>> LUND_OUTPUT;
std::vector<std::vector<double>> LUND_OUTPUT;

struct Particle
{
Expand All @@ -17,7 +15,7 @@ class Event
{
private:
double x, y, z, cross_section, beam_energy;
deque<Particle> bunch;
std::deque<Particle> bunch;
int beam_polarization;
bool cm_system;
public:
Expand All @@ -41,7 +39,7 @@ class Event
void add_particle(Particle&);
void clear_event(){bunch.clear();}
void print_vector();
void print_lund(string&);
void print_lund(std::string&);
void set_cm_system();
void set_lab_system();
void cm_to_lab(double&, double&);
Expand Down Expand Up @@ -109,7 +107,7 @@ void Event::add_particle(Particle& buff)

void Event::print_vector()
{
vector<double> buff;
std::vector<double> buff;
buff.push_back(bunch.size());
buff.push_back(1); buff.push_back(1); buff.push_back(0);
buff.push_back(beam_polarization); buff.push_back(11);
Expand All @@ -131,18 +129,18 @@ void Event::print_vector()
}
}

void Event::print_lund(string& Path)
void Event::print_lund(std::string& Path)
{
ofstream File;
std::ofstream File;

File.open(Path,fstream::in | fstream::out | fstream::app);
File.open(Path,std::fstream::in | std::fstream::out | std::fstream::app);

File << bunch.size() << "\t1\t1\t0\t" << beam_polarization << "\t11\t" << beam_energy << "\t2212\t0\t" << cross_section << endl;
File << bunch.size() << "\t1\t1\t0\t" << beam_polarization << "\t11\t" << beam_energy << "\t2212\t0\t" << cross_section << std::endl;

for(long unsigned int i = 0; i < bunch.size(); i++)
{
File << i+1 << "\t0\t1\t" << bunch[i].id << "\t0\t0\t" << (bunch[i].p).Px() << "\t" << (bunch[i].p).Py() << "\t"
<< (bunch[i].p).Pz() << "\t" << (bunch[i].p).E() << "\t" << bunch[i].mass << "\t" << x << "\t" << y << "\t" << z << endl;
<< (bunch[i].p).Pz() << "\t" << (bunch[i].p).E() << "\t" << bunch[i].mass << "\t" << x << "\t" << y << "\t" << z << std::endl;
}

File.close();
Expand Down
86 changes: 42 additions & 44 deletions general.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,15 @@
#include <getopt.h>
#include <fstream>

using namespace std;

double E0(6.5), R(1), L(10), W_min(1.08), W_max(2.0), Q2_min(0.05), Q2_max(5.0);
int h(0), N(1000000); double seed_(0);
bool channel(true), method(false), histogram(false), rad_corr(false);
string path = "./MCEGENpiN_radcorr.dat";
string source = "pi0p.csv";
string source_interp = "pi0p_int.csv";
vector<vector<double>> data, data_interp;
std::string path = "./MCEGENpiN_radcorr.dat";
std::string source = "pi0p.csv";
std::string source_interp = "pi0p_int.csv";
std::vector<std::vector<double>> data, data_interp;
double m_p(0.93827), m_n(0.93957), m_pip(0.13957), m_pi0(0.13498), m_e(0.000511), delta(0.01);
vector<double> values_rad{9999, 1, 9999, 1};
std::vector<double> values_rad{9999, 1, 9999, 1};

void add_hadrons(const double& W_, const double& Q2_, const double& theta, const double& phi, Event& bob); /* Adds hadron system in cm into Event object */
void add_electron(const double& W, const double& Q2, Event& bob); /* Adds outgoing electron in cm into Event object */
Expand All @@ -27,17 +25,17 @@ bool first_HM(true), truncate_out(false), extra_Q2(false);
int Total_Number(0), Accepted_Number(0);
/* ------------------- */

void Reading(string Path, vector<vector<double>>& V)
void Reading(std::string Path, std::vector<std::vector<double>>& V)
{
string line; stringstream ss;
ifstream File;
double dub; int w(1); vector<double> Numbers;
std::string line; std::stringstream ss;
std::ifstream File;
double dub; int w(1); std::vector<double> Numbers;

File.open(Path);

if (!File.is_open())
{
cout << "Can't open " << Path << " !" << endl;
std::cout << "Can't open " << Path << " !" << std::endl;
}
else
{
Expand Down Expand Up @@ -204,33 +202,33 @@ void input_check(int argc, char* argv[])
break;
};
case '?': default: {
cerr << "Unkhown option" << endl;
std::cerr << "Unkhown option" << std::endl;
break;
};
};
};

//path = data_path + "/" + path;

cout << " ------------------------------------------------------------------- " << endl;
cout << "| Monte Carlo event generator for exclusive pion electroproduction | \n| with radiative corrections \"MCEGENpiN_radcorr V9b\" | \n| |\n| Authors: Davydov M. - MSU, Physics dep. |\n| Isupov E. - MSU, SINP |\n| |\n| https://github.com/Maksaska/MCEGENpiN_radcorr |\n ------------------------------------------------------------------- " << endl;

cout << endl;

cout << "Beam energy is E = " << E0 << " GeV with";
if(h != 0){cout << " polarization h = " << h << endl;}
else{cout << " no polarization" << endl;}
cout << "Chosen kinematic area of (W, Q^2) values is W: " << W_min << " - " << W_max << " GeV , Q2: " << Q2_min << " - " << Q2_max << " GeV^2" << endl;
if(channel){cout << "Channel: pi0p with decay" << endl;}
else{cout << "Channel: pi+n" << endl;}
if(method) cout << "Uniform distribution with weights" << endl;
else cout << "Metropolis–Hastings MCMC algorithm" << endl;
cout << "Number of events: " << N << endl;
if(extra_Q2) cout << "All multiples (amplitudes) extrapolated as 1/Q4 (except M1+ ~ 1/Q6)" << endl;
else cout << "All multiples (amplitudes) extrapolated as 1/Q2 (except M1+ ~ 1/Q6)" << endl;
if(histogram){cout << "Histograms will be created\n" << endl;}
if(rad_corr){cout << "Radiative corrections: Enabled\n" << endl;}
else{cout << "Radiative corrections: Disabled\n" << endl;}
std::cout << " ------------------------------------------------------------------- " << std::endl;
std::cout << "| Monte Carlo event generator for exclusive pion electroproduction | \n| with radiative corrections \"MCEGENpiN_radcorr V9b\" | \n| |\n| Authors: Davydov M. - MSU, Physics dep. |\n| Isupov E. - MSU, SINP |\n| |\n| https://github.com/Maksaska/MCEGENpiN_radcorr |\n ------------------------------------------------------------------- " << std::endl;

std::cout << std::endl;

std::cout << "Beam energy is E = " << E0 << " GeV with";
if(h != 0){std::cout << " polarization h = " << h << std::endl;}
else{std::cout << " no polarization" << std::endl;}
std::cout << "Chosen kinematic area of (W, Q^2) values is W: " << W_min << " - " << W_max << " GeV , Q2: " << Q2_min << " - " << Q2_max << " GeV^2" << std::endl;
if(channel){std::cout << "Channel: pi0p with decay" << std::endl;}
else{std::cout << "Channel: pi+n" << std::endl;}
if(method) std::cout << "Uniform distribution with weights" << std::endl;
else std::cout << "Metropolis–Hastings MCMC algorithm" << std::endl;
std::cout << "Number of events: " << N << std::endl;
if(extra_Q2) std::cout << "All multiples (amplitudes) extrapolated as 1/Q4 (except M1+ ~ 1/Q6)" << std::endl;
else std::cout << "All multiples (amplitudes) extrapolated as 1/Q2 (except M1+ ~ 1/Q6)" << std::endl;
if(histogram){std::cout << "Histograms will be created\n" << std::endl;}
if(rad_corr){std::cout << "Radiative corrections: Enabled\n" << std::endl;}
else{std::cout << "Radiative corrections: Disabled\n" << std::endl;}

Reading(source, data);
Reading(source_interp, data_interp);
Expand Down Expand Up @@ -289,10 +287,10 @@ double P(const int& der,const int& n, const double& theta)
return Value;
}

vector<complex<double>> Finder(const double& W, const double& Q2)
std::vector<std::complex<double>> Finder(const double& W, const double& Q2)
{
vector<complex<double>> mult;
complex<double> buff;
std::vector<std::complex<double>> mult;
std::complex<double> buff;

for(long unsigned int i = 0; i < data.size(); i++)
{
Expand Down Expand Up @@ -339,10 +337,10 @@ vector<complex<double>> Finder(const double& W, const double& Q2)
return mult;
}

vector<complex<double>> Helicity_amplitudes(const double& W, const double& Q2, const double& theta)
std::vector<std::complex<double>> Helicity_amplitudes(const double& W, const double& Q2, const double& theta)
{
vector<complex<double>> AMP, H; double l;
complex<double> buff; buff = 0; bool extra_check(false);
std::vector<std::complex<double>> AMP, H; double l;
std::complex<double> buff; buff = 0; bool extra_check(false);

double factor_all(1), factor_M1(1), Q2_1(0);
Q2_1 = Q2;
Expand Down Expand Up @@ -427,7 +425,7 @@ vector<complex<double>> Helicity_amplitudes(const double& W, const double& Q2,
double Section(const double& W, const double& Q2, const double& theta, const double& phi)
{
double S_t, S_l, S_tt, S_lt, S_lt_pr, S, Gamma_flux, eps, nu, L(0.14817);
vector<complex<double>> H;
std::vector<std::complex<double>> H;

nu = (W*W + Q2 - m_p*m_p)/(2*m_p);
eps = 1/(1 + 2*(nu*nu + Q2)/(4*(E0 - nu)*E0 - Q2));
Expand All @@ -448,9 +446,9 @@ double Section(const double& W, const double& Q2, const double& theta, const d
return S;
}

vector<double> Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2)
std::vector<double> Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2)
{
vector<double> Abc(2);
std::vector<double> Abc(2);

Abc[0] = (y1 - y2)/(x1 - x2);
Abc[1] = (x1*y2 - y1*x2)/(x1 - x2);
Expand All @@ -460,7 +458,7 @@ vector<double> Coefficients_lin(const double& x1, const double& y1, const double

double Linear(const double& W, const double& Q2, const double& theta, const double& phi)
{
double S; vector<double> ab;
double S; std::vector<double> ab;
double W_min_d, W_max_d, Q2_min_d, Q2_max_d, S1, S2, S3, S4, S_1, S_2;

W_min_d = floor(W*100)/100;
Expand Down Expand Up @@ -519,7 +517,7 @@ double Section_int(const double& W, const double& Q2, const double& E_ini)

double Section_interp_int(const double& W, const double& Q2, const double& E_ini)
{
double S; vector<double> ab;
double S; std::vector<double> ab;
double W_min_d, W_max_d, Q2_min_d, Q2_max_d, S1, S2, S3, S4, S_1, S_2;

W_min_d = floor(W*100)/100;
Expand Down Expand Up @@ -943,7 +941,7 @@ bool check_ps(const double& W_, const double& Q2_, const double& E_beam)
return check;
}

double Normal(const double& mean, const double& sigma, const string& str)
double Normal(const double& mean, const double& sigma, const std::string& str)
{
double func(-1), rand_v(0);
double value = fRand(0, 1);
Expand Down
16 changes: 7 additions & 9 deletions header.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,14 @@
#include "TH2F.h"
#include <fstream>

using namespace std;

extern double E0, R, L, W_min, W_max, Q2_min, Q2_max, delta;
extern vector<vector<double>> data, data_interp, LUND_OUTPUT;
extern std::vector<std::vector<double>> data, data_interp, LUND_OUTPUT;
extern int h, N;
extern double seed_;
extern string path, source, source_interp;
extern std::string path, source, source_interp;
extern bool channel, method, histogram, rad_corr, truncate_out, extra_Q2;
extern double m_p, m_n, m_pip, m_pi0, m_e;
extern vector<double> values_rad;
extern std::vector<double> values_rad;

extern double Q2_line, W_line;

Expand Down Expand Up @@ -50,9 +48,9 @@ double sin_2(const double& W, const double& Q2); /* sin(theta)_2 - for polar a
/* -------- general.cpp -------- */

void input_check(int argc, char* argv[]); /* This function checks if there any option were passed in main() */
void Reading(string Path,vector<vector<double>>&V); /* Reads the data from csv */
void Reading(std::string Path,std::vector<std::vector<double>>&V); /* Reads the data from csv */
double P(const int& der,const int& n, const double& theta); /* Legendre polynomials */
vector<complex<double>> Finder(const double& W, const double& Q2); /* it searches for multipoles EMS */
std::vector<std::complex<double>> Finder(const double& W, const double& Q2); /* it searches for multipoles EMS */
void generate_particle(const int& k); /*
* This function generates the event
* It creates kin. parameters from given E0 [W_min, W_max] && [Q2_min, Q2_max]
Expand All @@ -62,13 +60,13 @@ void generate_particle(const int& k); /*
* Writes the event in output file
*/

vector<double> Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2); /* Coef. for linear interp. */
std::vector<double> Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2); /* Coef. for linear interp. */
double Linear(const double& W, const double& Q2, const double& theta, const double& phi); /* linear interpolation */
double Section(const double& W, const double& Q2, const double& theta, const double& phi); /*
Cross_section evaluation from Helicity ampl.
*/

vector<complex<double>> Helicity_amplitudes(const double& W, const double& Q2, const double& theta);/* Helicity ampl. eval. */
std::vector<std::complex<double>> Helicity_amplitudes(const double& W, const double& Q2, const double& theta);/* Helicity ampl. eval. */
double Section_int(const double& W, const double& Q2, const double& E_beam); /* Integral cross section dS/dOmegadE from dataset */
double Section_interp_int(const double& W, const double& Q2, const double& E_beam); /* Integral cross section dS/dOmegadE for random W, Q2 */
double Spence(const double& x); /* Spence function */
Expand Down
3 changes: 0 additions & 3 deletions kinematics.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,5 @@
#include <iostream>
#include "header.h"

using namespace std;

double fRand(const double& fMin, const double& fMax)
{
double f = (double)rand() / RAND_MAX;
Expand Down
26 changes: 12 additions & 14 deletions main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
#include "header.h"

using namespace std;

auto c1 = new TCanvas("c1", "Histogram", 1280, 1080);
TH2* h1 = new TH2F("h1", "Histogram (W,Q^{2})_2d", 9600, 0.0, 40.0, 3200, 0.0, 40.0);
TH1F* h3 = new TH1F("h3", "Histogram W", 9600, 0.0, 40.0);
Expand All @@ -21,8 +19,8 @@ int main(int argc, char* argv[])
srand(seed_);
if(rad_corr){counter_ = 5; inter = 3;}

ofstream File;
string PATH;
std::ofstream File;
std::string PATH;
int id_dot = path.find_last_of('.');

if(truncate_out)
Expand All @@ -39,20 +37,20 @@ int main(int argc, char* argv[])
PATH = path;
PATH.insert(id_dot,std::to_string(counter));

File.open(PATH,fstream::in | fstream::out | fstream::app);
File.open(PATH,std::fstream::in | std::fstream::out | std::fstream::app);

for(auto i:LUND_OUTPUT)
{
for(auto j:i)
{
File << j << "\t";
}
File << endl;
File << std::endl;
}

finish_ = std::chrono::high_resolution_clock::now();
elapsed_ = double(ceil(double(N)/10000.0) - 1 - counter)*(finish_ - start)/double(counter);
cout << "Stand by ..." << ceil(counter*10000.0/(ceil(double(N)/10000.0) - 1))/100.0 << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << flush;
std::cout << "Stand by ..." << ceil(counter*10000.0/(ceil(double(N)/10000.0) - 1))/100.0 << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << std::flush;

File.close(); LUND_OUTPUT.clear();
}
Expand All @@ -65,28 +63,28 @@ int main(int argc, char* argv[])
{
finish_ = std::chrono::high_resolution_clock::now();
elapsed_ = double(N - k)*(finish_ - start)/double(k);
if(!method) cout << "Acceptance Rate: " << floor(10000*(double(Accepted_Number)/double(Total_Number)))/100 << "%\tStand by ..." << k*100/N << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << flush;
else cout << "Stand by ..." << k*100/N << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << flush;
if(!method) std::cout << "Acceptance Rate: " << floor(10000*(double(Accepted_Number)/double(Total_Number)))/100 << "%\tStand by ..." << k*100/N << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << std::flush;
else std::cout << "Stand by ..." << k*100/N << "% Time remaining: " << floor(elapsed_.count()/3600) << " h " << floor((elapsed_.count() - 3600*floor(elapsed_.count()/3600))/60) << " min " << elapsed_.count() - 60*floor(elapsed_.count()/60) << " s \r" << std::flush;
}
}

File.open(path,fstream::in | fstream::out | fstream::app);
File.open(path,std::fstream::in | std::fstream::out | std::fstream::app);

for(int i = 0; i < LUND_OUTPUT.size(); i++)
{
for(auto j:LUND_OUTPUT[i])
{
File << j << "\t";
}
File << endl;
File << std::endl;

if(i % counter_ == 0) cout << "Recording LUND ... " << ceil((i+1)*10000.0/LUND_OUTPUT.size())/100.0 << "% \r" << flush;
if(i % counter_ == 0) std::cout << "Recording LUND ... " << ceil((i+1)*10000.0/LUND_OUTPUT.size())/100.0 << "% \r" << std::flush;
}

File.close();
}

cout << "Recording complete " << endl;
std::cout << "Recording complete " << std::endl;

if(histogram)
{
Expand Down Expand Up @@ -135,7 +133,7 @@ int main(int argc, char* argv[])

auto finish = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = finish - start;
cout << "Elapsed time: " << floor(elapsed.count()/3600) << " h " << floor((elapsed.count() - 3600*floor(elapsed.count()/3600))/60) << " min " << elapsed.count() - 60*floor(elapsed.count()/60) << " s\n";
std::cout << "Elapsed time: " << floor(elapsed.count()/3600) << " h " << floor((elapsed.count() - 3600*floor(elapsed.count()/3600))/60) << " min " << elapsed.count() - 60*floor(elapsed.count()/60) << " s\n";

data.clear(); data_interp.clear(); values_rad.clear(), LUND_OUTPUT.clear();

Expand Down

0 comments on commit e812c6f

Please sign in to comment.