diff --git a/event.h b/event.h index 72057f7..a9f0e9a 100644 --- a/event.h +++ b/event.h @@ -2,9 +2,7 @@ #include #include -using namespace std; - -vector> LUND_OUTPUT; +std::vector> LUND_OUTPUT; struct Particle { @@ -17,7 +15,7 @@ class Event { private: double x, y, z, cross_section, beam_energy; - deque bunch; + std::deque bunch; int beam_polarization; bool cm_system; public: @@ -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&); @@ -109,7 +107,7 @@ void Event::add_particle(Particle& buff) void Event::print_vector() { - vector buff; + std::vector 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); @@ -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(); diff --git a/general.cpp b/general.cpp index b48ae0e..e650ed2 100644 --- a/general.cpp +++ b/general.cpp @@ -3,17 +3,15 @@ #include #include -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> data, data_interp; +std::string path = "./MCEGENpiN_radcorr.dat"; +std::string source = "pi0p.csv"; +std::string source_interp = "pi0p_int.csv"; +std::vector> 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 values_rad{9999, 1, 9999, 1}; +std::vector 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 */ @@ -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>& V) +void Reading(std::string Path, std::vector>& V) { - string line; stringstream ss; - ifstream File; - double dub; int w(1); vector Numbers; + std::string line; std::stringstream ss; + std::ifstream File; + double dub; int w(1); std::vector Numbers; File.open(Path); if (!File.is_open()) { - cout << "Can't open " << Path << " !" << endl; + std::cout << "Can't open " << Path << " !" << std::endl; } else { @@ -204,7 +202,7 @@ void input_check(int argc, char* argv[]) break; }; case '?': default: { - cerr << "Unkhown option" << endl; + std::cerr << "Unkhown option" << std::endl; break; }; }; @@ -212,25 +210,25 @@ void input_check(int argc, char* argv[]) //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); @@ -289,10 +287,10 @@ double P(const int& der,const int& n, const double& theta) return Value; } -vector> Finder(const double& W, const double& Q2) +std::vector> Finder(const double& W, const double& Q2) { - vector> mult; - complex buff; + std::vector> mult; + std::complex buff; for(long unsigned int i = 0; i < data.size(); i++) { @@ -339,10 +337,10 @@ vector> Finder(const double& W, const double& Q2) return mult; } -vector> Helicity_amplitudes(const double& W, const double& Q2, const double& theta) +std::vector> Helicity_amplitudes(const double& W, const double& Q2, const double& theta) { - vector> AMP, H; double l; - complex buff; buff = 0; bool extra_check(false); + std::vector> AMP, H; double l; + std::complex buff; buff = 0; bool extra_check(false); double factor_all(1), factor_M1(1), Q2_1(0); Q2_1 = Q2; @@ -427,7 +425,7 @@ vector> 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> H; + std::vector> H; nu = (W*W + Q2 - m_p*m_p)/(2*m_p); eps = 1/(1 + 2*(nu*nu + Q2)/(4*(E0 - nu)*E0 - Q2)); @@ -448,9 +446,9 @@ double Section(const double& W, const double& Q2, const double& theta, const d return S; } -vector Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2) +std::vector Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2) { - vector Abc(2); + std::vector Abc(2); Abc[0] = (y1 - y2)/(x1 - x2); Abc[1] = (x1*y2 - y1*x2)/(x1 - x2); @@ -460,7 +458,7 @@ vector 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 ab; + double S; std::vector 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; @@ -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 ab; + double S; std::vector 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; @@ -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); diff --git a/header.h b/header.h index 465d312..729a08c 100644 --- a/header.h +++ b/header.h @@ -12,16 +12,14 @@ #include "TH2F.h" #include -using namespace std; - extern double E0, R, L, W_min, W_max, Q2_min, Q2_max, delta; -extern vector> data, data_interp, LUND_OUTPUT; +extern std::vector> 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 values_rad; +extern std::vector values_rad; extern double Q2_line, W_line; @@ -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>&V); /* Reads the data from csv */ +void Reading(std::string Path,std::vector>&V); /* Reads the data from csv */ double P(const int& der,const int& n, const double& theta); /* Legendre polynomials */ -vector> Finder(const double& W, const double& Q2); /* it searches for multipoles EMS */ +std::vector> 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] @@ -62,13 +60,13 @@ void generate_particle(const int& k); /* * Writes the event in output file */ -vector Coefficients_lin(const double& x1, const double& y1, const double& x2, const double& y2); /* Coef. for linear interp. */ +std::vector 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> Helicity_amplitudes(const double& W, const double& Q2, const double& theta);/* Helicity ampl. eval. */ +std::vector> 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 */ diff --git a/kinematics.cpp b/kinematics.cpp index 9d294b5..b769206 100644 --- a/kinematics.cpp +++ b/kinematics.cpp @@ -1,8 +1,5 @@ -#include #include "header.h" -using namespace std; - double fRand(const double& fMin, const double& fMax) { double f = (double)rand() / RAND_MAX; diff --git a/main.cpp b/main.cpp index 83b35b7..2d29d4a 100644 --- a/main.cpp +++ b/main.cpp @@ -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); @@ -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) @@ -39,7 +37,7 @@ 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) { @@ -47,12 +45,12 @@ int main(int argc, char* argv[]) { 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(); } @@ -65,12 +63,12 @@ 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++) { @@ -78,15 +76,15 @@ int main(int argc, char* argv[]) { 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) { @@ -135,7 +133,7 @@ int main(int argc, char* argv[]) auto finish = std::chrono::high_resolution_clock::now(); std::chrono::duration 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();