From dc7954aab9f854b36b62f1ed3c5dcfea4de0bd58 Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:26:46 -0300 Subject: [PATCH 01/17] Add VChamberShape type, add xy loss option and fix identations --- include/trackcpp/auxiliary.h | 130 ++++++++++++++++++----------------- 1 file changed, 68 insertions(+), 62 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 5496446..077d7d3 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -26,33 +26,33 @@ class PassMethodsClass { public: - static const int pm_identity_pass = 0; - static const int pm_drift_pass = 1; - static const int pm_str_mpole_symplectic4_pass = 2; - static const int pm_bnd_mpole_symplectic4_pass = 3; - static const int pm_corrector_pass = 4; - static const int pm_cavity_pass = 5; - static const int pm_thinquad_pass = 6; - static const int pm_thinsext_pass = 7; - static const int pm_kickmap_pass = 8; - static const int pm_matrix_pass = 9; - static const int pm_nr_pms = 10; // counter for number of passmethods - PassMethodsClass() { - passmethods.push_back("identity_pass"); - passmethods.push_back("drift_pass"); - passmethods.push_back("str_mpole_symplectic4_pass"); - passmethods.push_back("bnd_mpole_symplectic4_pass"); - passmethods.push_back("corrector_pass"); - passmethods.push_back("cavity_pass"); - passmethods.push_back("thinquad_pass"); - passmethods.push_back("thinsext_pass"); - passmethods.push_back("kicktable_pass"); - passmethods.push_back("matrix_pass"); - } - int size() const { return passmethods.size(); } - std::string operator[](const int i) const { return passmethods[i]; } + static const int pm_identity_pass = 0; + static const int pm_drift_pass = 1; + static const int pm_str_mpole_symplectic4_pass = 2; + static const int pm_bnd_mpole_symplectic4_pass = 3; + static const int pm_corrector_pass = 4; + static const int pm_cavity_pass = 5; + static const int pm_thinquad_pass = 6; + static const int pm_thinsext_pass = 7; + static const int pm_kickmap_pass = 8; + static const int pm_matrix_pass = 9; + static const int pm_nr_pms = 10; // counter for number of passmethods + PassMethodsClass() { + passmethods.push_back("identity_pass"); + passmethods.push_back("drift_pass"); + passmethods.push_back("str_mpole_symplectic4_pass"); + passmethods.push_back("bnd_mpole_symplectic4_pass"); + passmethods.push_back("corrector_pass"); + passmethods.push_back("cavity_pass"); + passmethods.push_back("thinquad_pass"); + passmethods.push_back("thinsext_pass"); + passmethods.push_back("kicktable_pass"); + passmethods.push_back("matrix_pass"); + } + int size() const { return passmethods.size(); } + std::string operator[](const int i) const { return passmethods[i]; } private: - std::vector passmethods; + std::vector passmethods; }; // this is to be superseede by PassMethodClass @@ -74,17 +74,17 @@ struct PassMethod { // this is to be superseede by PassMethodClass const std::vector pm_dict = { - // this vector has to have the same number of entries as enum above - "identity_pass", - "drift_pass", - "str_mpole_symplectic4_pass", - "bnd_mpole_symplectic4_pass", - "corrector_pass", - "cavity_pass", - "thinquad_pass", - "thinsext_pass", - "kicktable_pass", - "matrix_pass" + // this vector has to have the same number of entries as enum above + "identity_pass", + "drift_pass", + "str_mpole_symplectic4_pass", + "bnd_mpole_symplectic4_pass", + "corrector_pass", + "cavity_pass", + "thinquad_pass", + "thinsext_pass", + "kicktable_pass", + "matrix_pass" }; struct Status { @@ -108,21 +108,21 @@ struct Status { }; const std::vector string_error_messages = { - "success", - "passmethod_not_defined", - "passmethod_not_implemented", - "particle_lost", - "inconsistent_dimensions", - "uninitialized_memory", - "findorbit_not_converged", - "findorbit_one_turn_matrix_problem", - "file_not_found", - "file_not_opened", - "kicktable_not_defined", - "kicktable_out_of_range", - "flat_file_error", - "newton_not_converged", - "not_implemented", + "success", + "passmethod_not_defined", + "passmethod_not_implemented", + "particle_lost", + "inconsistent_dimensions", + "uninitialized_memory", + "findorbit_not_converged", + "findorbit_one_turn_matrix_problem", + "file_not_found", + "file_not_opened", + "kicktable_not_defined", + "kicktable_out_of_range", + "flat_file_error", + "newton_not_converged", + "not_implemented", }; #define STR_HELPER(x) #x @@ -132,16 +132,22 @@ const std::vector string_error_messages = { const std::string string_version = "TRACKCPP version " + std::string(VER_STR) + " (" + std::string(__DATE__) + " " + std::string(__TIME__) + ")"; struct Plane { - enum type { - no_plane = 0, - x = 1, - y = 2, - z = 3 - }; + enum type { + no_plane = 0, + x = 1, + y = 2, + z = 3, + xy = 4 + }; }; -// template class Pos; -// class Element; +struct VChamberShape { + enum type { + rectangle = 0, + kite = 1, + ellipse = 2, + }; +}; const double light_speed = 299792458; // [m/s] - definition const double electron_charge = 1.602176634e-19; // [C] - definition @@ -155,12 +161,12 @@ const double electron_radius = pow(electron_charge,2)/(4*M_PI*vaccum_pe template int sgn(T val) { - if (val >= 0) return 1; else return -1; + if (val >= 0) return 1; else return -1; } #if __GNUC__ < 6 - bool isfinite(const double& v); + bool isfinite(const double& v); #endif std::string get_timestamp(); From bc1e48d46249ba70e62953e380fdf375e2887ecb Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:36:34 -0300 Subject: [PATCH 02/17] Add logical comparisons of Tpsa objects --- include/trackcpp/tpsa.h | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/include/trackcpp/tpsa.h b/include/trackcpp/tpsa.h index f9480f4..0ae3748 100644 --- a/include/trackcpp/tpsa.h +++ b/include/trackcpp/tpsa.h @@ -140,10 +140,12 @@ class Tpsa { bool operator > (const TYPE& o_) const; bool operator >= (const TYPE& o_) const; - bool operator > (const Tpsa& o_) const; - bool operator >= (const Tpsa& o_) const; bool operator == (const Tpsa& o_) const; bool operator != (const Tpsa& o_) const; + bool operator < (const Tpsa& o_) const; + bool operator <= (const Tpsa& o_) const; + bool operator > (const Tpsa& o_) const; + bool operator >= (const Tpsa& o_) const; explicit operator int() const { return int(c[0]); } explicit operator double() const { return double(c[0]); } @@ -349,6 +351,7 @@ bool Tpsa::operator < (const TYPE& o_) const { return false; } + template bool Tpsa::operator <= (const TYPE& o_) const { for(unsigned int i=0; i::operator >= (const TYPE& o_) const { return false; } -template -bool Tpsa::operator >= (const Tpsa& o_) const { - return (*this - o_) >= (TYPE) 0; -} - template bool Tpsa::operator == (const Tpsa& o_) const { return (*this - o_) == (TYPE) 0; @@ -382,11 +380,26 @@ bool Tpsa::operator != (const Tpsa& o_) const { return (*this - o_) != (TYPE) 0; } +template +bool Tpsa::operator < (const Tpsa& o_) const { + return (*this - o_) < (TYPE) 0; +} + +template +bool Tpsa::operator <= (const Tpsa& o_) const { + return (*this - o_) <= (TYPE) 0; +} + template bool Tpsa::operator > (const Tpsa& o_) const { return (*this - o_) > (TYPE) 0; } +template +bool Tpsa::operator >= (const Tpsa& o_) const { + return (*this - o_) >= (TYPE) 0; +} + // Implementation: AUXILIARY MEMBER FUNCTIONS // ------------------------------------------ From 13160aee0a44b38035ad5d8c0fe91e06459ffcc3 Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:38:28 -0300 Subject: [PATCH 03/17] Add 'vchamber' shape attribute to Element --- include/trackcpp/elements.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 357f6fe..2e99fa2 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -43,6 +43,7 @@ class Element { int pass_method = PassMethod::pm_drift_pass; double length = 0; // [m] int nr_steps = 1; // + int vchamber = VChamberShape::rectangle; // vchamber type double hmin = -DBL_MAX; // [m] double hmax = DBL_MAX; // [m] double vmin = -DBL_MAX; // [m] @@ -60,8 +61,8 @@ class Element { double frequency = 0; // [Hz] double voltage = 0; // [V] double phase_lag = 0; // [rad] - int kicktable_idx = -1; // index of kickmap object in kicktable_list - double rescale_kicks = 1.0; // for kickmaps + int kicktable_idx = -1; // index of kickmap object in kicktable_list + double rescale_kicks = 1.0; // for kickmaps std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; From 7873c98373c48ef429c3257812e0cbb297886b7b Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:39:18 -0300 Subject: [PATCH 04/17] Implement different vchamber shape collision checks --- include/trackcpp/tracking.h | 94 +++++++++++++++++++++++++++++++++---- 1 file changed, 84 insertions(+), 10 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index a734c56..9e28de6 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -163,19 +163,69 @@ Status::type track_linepass ( status = track_elementpass (element, orig_pos, accelerator); + const T& rx = orig_pos.rx; + const T& ry = orig_pos.ry; + // checks if particle is lost - if ((not isfinite(orig_pos.rx)) or - ((accelerator.vchamber_on) and - ((orig_pos.rx < element.hmin) or - (orig_pos.rx > element.hmax)))) { + + if (not isfinite(rx)) { lost_plane = Plane::x; status = Status::particle_lost; - }else if ((not isfinite(orig_pos.ry)) or - ((accelerator.vchamber_on) and - ((orig_pos.ry < element.vmin) or - (orig_pos.ry > element.vmax)))) { - lost_plane = Plane::y; - status = Status::particle_lost; + } + if (not isfinite(ry)) { + if (status != Status::particle_lost) { + lost_plane = Plane::y; + status = Status::particle_lost; + } else { + lost_plane = Plane::xy; + } + } + if ((status != Status::particle_lost) and accelerator.vchamber_on) { + if (element.vchamber == VChamberShape::rectangle) { + // rectangular vaccum chamber + if (((rx < element.hmin) or (rx > element.hmax))) { + lost_plane = Plane::x; + status = Status::particle_lost; + } + if (((ry < element.vmin) or (ry > element.vmax))) { + if (status != Status::particle_lost) { + lost_plane = Plane::y; + status = Status::particle_lost; + } else { + lost_plane = Plane::xy; + } + } + } else if (element.vchamber == VChamberShape::kite) { + // kite-shaped vaccum chamber + if (((rx < element.hmin) or (rx > element.hmax))) { + lost_plane = Plane::xy; + status = Status::particle_lost; + } + if ((ry > get_kite_ry(element, 0, rx)) or // upper right + (ry > get_kite_ry(element, 1, rx)) or // upper left + (ry < get_kite_ry(element, 2, rx)) or // lower left + (ry < get_kite_ry(element, 3, rx))) { // lower right + lost_plane = Plane::xy; + status = Status::particle_lost; + } + } else if (element.vchamber == VChamberShape::ellipse) { + // elliptical vaccum chamber + double lx = (element.hmax - element.hmin) / 2; + double ly = (element.vmax - element.vmin) / 2; + double xc = (element.hmax + element.hmin) / 2; + double yc = (element.vmax + element.vmin) / 2; + double xn = (rx - xc)/lx; + double yn = (ry - yc)/ly; + double amplitude = xn*xn + yn*yn; + if (amplitude > 1) { + lost_plane = Plane::xy; + status = Status::particle_lost; + } + } else { + // any other shape not implemented safely signals lost particle + lost_plane = Plane::xy; + status = Status::particle_lost; + } } if (status != Status::success) { @@ -197,6 +247,30 @@ Status::type track_linepass ( } +template +T get_kite_ry(const Element& elem, int side, const T& rx) { + double x1, y1, x2, y2; + if (side == 0) { + // upper-right + x1 = 0; y1 = elem.vmax; + x2 = elem.hmax; y2 = 0; + } else if (side == 1) { + // upper-left + x1 = elem.hmin; y1 = 0; + x2 = 0; y1 = elem.vmax; + } else if (side == 2) { + // lower-left + x1 = elem.hmin; y1 = 0; + x2 = 0; y1 = elem.vmin; + } else { + // lower-right + x1 = 0; y1 = elem.vmin; + x2 = elem.hmax; y1 = 0; + } + return y1 + (rx - x1) * (y2 - y1) / (x2 - x1); +} + + template Status::type track_linepass ( const Accelerator& accelerator, From 22edf9493e8ed024a317e0744967ceb3b9b569cd Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:40:10 -0300 Subject: [PATCH 05/17] Update VERSION --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 1d068c6..a84947d 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -4.4.2 +4.5.0 From 0de9a811c5765a57a9e4b36f9194ba318e24fe48 Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 13:54:01 -0300 Subject: [PATCH 06/17] Fix elliptical vchamber check --- include/trackcpp/tracking.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 9e28de6..1dd146f 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -214,9 +214,9 @@ Status::type track_linepass ( double ly = (element.vmax - element.vmin) / 2; double xc = (element.hmax + element.hmin) / 2; double yc = (element.vmax + element.vmin) / 2; - double xn = (rx - xc)/lx; - double yn = (ry - yc)/ly; - double amplitude = xn*xn + yn*yn; + T xn = (rx - xc)/lx; + T yn = (ry - yc)/ly; + T amplitude = xn*xn + yn*yn; if (amplitude > 1) { lost_plane = Plane::xy; status = Status::particle_lost; From de8f65aff46192b18573874c958fc05a9690c8ae Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 15:27:23 -0300 Subject: [PATCH 07/17] Add 'vchamber' attribute o flat_files functions --- src/flat_file.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 8031f94..9926673 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -110,6 +110,7 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) } if (has_polynom(e.polynom_a)) write_polynom(fp, "polynom_a", e.polynom_a); if (has_polynom(e.polynom_b)) write_polynom(fp, "polynom_b", e.polynom_b); + if (e.vchamber != 0) { fp << std::setw(pw) << "vchamber" << e.vchamber << '\n'; } if (e.hmin != 0) { fp << std::setw(pw) << "hmin" << e.hmin << '\n'; } if (e.hmax != 0) { fp << std::setw(pw) << "hmax" << e.hmax << '\n'; } if (e.vmin != 0) { fp << std::setw(pw) << "vmin" << e.vmin << '\n'; } @@ -195,6 +196,7 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } if (cmd.compare("length") == 0) { ss >> e.length; continue; } + if (cmd.compare("vchamber") == 0) { ss >> e.vchamber; continue; } if (cmd.compare("hmin") == 0) { ss >> e.hmin; found_hmin = true; continue; } if (cmd.compare("hmax") == 0) { ss >> e.hmax; From 6991383ae610eb7eee9d5fc8e022521a3b653a6a Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 15:28:20 -0300 Subject: [PATCH 08/17] Add 'vchamber' and 'rescale_kicks' to element comparison --- src/elements.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/elements.cpp b/src/elements.cpp index 075e285..3f58b93 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -150,6 +150,7 @@ bool Element::operator==(const Element& o) const { if (this->fam_name != o.fam_name) return false; if (this->pass_method != o.pass_method) return false; if (this->length != o.length) return false; + if (this->vchamber != o.vchamber) return false; if (this->hmin != o.hmin) return false; if (this->hmax != o.hmax) return false; if (this->vmin != o.vmin) return false; @@ -191,6 +192,8 @@ bool Element::operator==(const Element& o) const { if (this->r_out[i] != o.r_out[i]) return false; } if (this->kicktable_idx != o.kicktable_idx) return false; + if (this->rescale_kicks != o.rescale_kicks) return false; + return true; } From a64954f2c1a1f3dc1eb55d660723fa8c3d7fc9bb Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 17:01:50 -0300 Subject: [PATCH 09/17] Update VERSION --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index a84947d..f6cdf40 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -4.5.0 +4.7.0 From c9b4435544d00a025d5ec60aafadefddab9e0ac5 Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 17:31:21 -0300 Subject: [PATCH 10/17] Fix bug in collision check for kyte shaped vchambers --- include/trackcpp/tracking.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 1dd146f..dcb4467 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -251,21 +251,21 @@ template T get_kite_ry(const Element& elem, int side, const T& rx) { double x1, y1, x2, y2; if (side == 0) { - // upper-right + // std::cout << "upper-right" << std::endl; x1 = 0; y1 = elem.vmax; x2 = elem.hmax; y2 = 0; } else if (side == 1) { - // upper-left + // std::cout << "upper-left" << std::endl; x1 = elem.hmin; y1 = 0; - x2 = 0; y1 = elem.vmax; + x2 = 0; y2 = elem.vmax; } else if (side == 2) { - // lower-left + // std::cout << "lower-left" << std::endl; x1 = elem.hmin; y1 = 0; - x2 = 0; y1 = elem.vmin; + x2 = 0; y2 = elem.vmin; } else { - // lower-right + // std::cout << "lower-right" << std::endl; x1 = 0; y1 = elem.vmin; - x2 = elem.hmax; y1 = 0; + x2 = elem.hmax; y2 = 0; } return y1 + (rx - x1) * (y2 - y1) / (x2 - x1); } From deb0601bd225157183cc8b4414e95ae6b8c2756d Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 17:46:43 -0300 Subject: [PATCH 11/17] Fix vacuum spelling --- include/trackcpp/auxiliary.h | 4 ++-- include/trackcpp/tracking.h | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 077d7d3..75be914 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -155,9 +155,9 @@ const double reduced_planck_constant = 1.054571817e-34; // [J.s] - definiti const double electron_mass = 9.1093837015e-31; // [Kg] - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=electron+mass const double vacuum_permeability = 1.25663706212e-6; // [T.m/A] - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?mu0|search_for=vacuum+permeability const double electron_rest_energy = electron_mass * pow(light_speed,2); // [Kg.m^2/s^2] - derived -const double vaccum_permitticity = 1/(vacuum_permeability * pow(light_speed,2)); // [V.s/(A.m)] - derived +const double vacuum_permitticity = 1/(vacuum_permeability * pow(light_speed,2)); // [V.s/(A.m)] - derived const double electron_rest_energy_MeV = (electron_rest_energy / electron_charge) / 1e6; // [MeV] - derived -const double electron_radius = pow(electron_charge,2)/(4*M_PI*vaccum_permitticity*electron_rest_energy); // [m] - derived +const double electron_radius = pow(electron_charge,2)/(4*M_PI*vacuum_permitticity*electron_rest_energy); // [m] - derived template int sgn(T val) { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index dcb4467..dacc6b8 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -182,7 +182,7 @@ Status::type track_linepass ( } if ((status != Status::particle_lost) and accelerator.vchamber_on) { if (element.vchamber == VChamberShape::rectangle) { - // rectangular vaccum chamber + // rectangular vacuum chamber if (((rx < element.hmin) or (rx > element.hmax))) { lost_plane = Plane::x; status = Status::particle_lost; @@ -196,7 +196,7 @@ Status::type track_linepass ( } } } else if (element.vchamber == VChamberShape::kite) { - // kite-shaped vaccum chamber + // kite-shaped vacuum chamber if (((rx < element.hmin) or (rx > element.hmax))) { lost_plane = Plane::xy; status = Status::particle_lost; @@ -209,7 +209,7 @@ Status::type track_linepass ( status = Status::particle_lost; } } else if (element.vchamber == VChamberShape::ellipse) { - // elliptical vaccum chamber + // elliptical vacuum chamber double lx = (element.hmax - element.hmin) / 2; double ly = (element.vmax - element.vmin) / 2; double xc = (element.hmax + element.hmin) / 2; From c38f30ce133a843204aee1a11b2ec8492e72ce03 Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 19:02:27 -0300 Subject: [PATCH 12/17] check physical constants --- include/trackcpp/auxiliary.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 75be914..e884e3d 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -152,8 +152,8 @@ struct VChamberShape { const double light_speed = 299792458; // [m/s] - definition const double electron_charge = 1.602176634e-19; // [C] - definition const double reduced_planck_constant = 1.054571817e-34; // [J.s] - definition -const double electron_mass = 9.1093837015e-31; // [Kg] - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=electron+mass -const double vacuum_permeability = 1.25663706212e-6; // [T.m/A] - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?mu0|search_for=vacuum+permeability +const double electron_mass = 9.1093837015e-31; // [Kg] - 2022-03-19 - https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=electron+mass +const double vacuum_permeability = 1.25663706212e-6; // [T.m/A] - 2022-03-19 - https://physics.nist.gov/cgi-bin/cuu/Value?mu0|search_for=vacuum+permeability const double electron_rest_energy = electron_mass * pow(light_speed,2); // [Kg.m^2/s^2] - derived const double vacuum_permitticity = 1/(vacuum_permeability * pow(light_speed,2)); // [V.s/(A.m)] - derived const double electron_rest_energy_MeV = (electron_rest_energy / electron_charge) / 1e6; // [MeV] - derived From f50d9b5f2e5a774656c0a3be26111d3da6a40bef Mon Sep 17 00:00:00 2001 From: ximenes Date: Sat, 19 Mar 2022 19:16:45 -0300 Subject: [PATCH 13/17] Add comment line with a URL about SI units --- include/trackcpp/auxiliary.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index e884e3d..aa8b3cc 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -149,6 +149,8 @@ struct VChamberShape { }; }; +// See https://en.wikipedia.org/wiki/International_System_of_Units + const double light_speed = 299792458; // [m/s] - definition const double electron_charge = 1.602176634e-19; // [C] - definition const double reduced_planck_constant = 1.054571817e-34; // [J.s] - definition From c49edd3fd82702df77d4cded27528bef4a0ee71a Mon Sep 17 00:00:00 2001 From: ximenes Date: Tue, 22 Mar 2022 12:52:27 -0300 Subject: [PATCH 14/17] Generalize vchamber --- include/trackcpp/auxiliary.h | 5 +-- include/trackcpp/tpsa.h | 17 +++++++++-- include/trackcpp/tracking.h | 59 +++++++++++------------------------- 3 files changed, 35 insertions(+), 46 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index aa8b3cc..98a5e78 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -142,10 +142,11 @@ struct Plane { }; struct VChamberShape { + // integer corresponds to p-norm number defining the normalized shape enum type { - rectangle = 0, - kite = 1, + rhombus = 1, ellipse = 2, + rectangle = 0, // (inf-norm) }; }; diff --git a/include/trackcpp/tpsa.h b/include/trackcpp/tpsa.h index 0ae3748..c0b7116 100644 --- a/include/trackcpp/tpsa.h +++ b/include/trackcpp/tpsa.h @@ -70,6 +70,7 @@ template struct et_osip { enum { val = et_binomia // -------------------- template class Tpsa; +template Tpsa pow (const Tpsa&, const int n); template Tpsa abs (const Tpsa&); template Tpsa sqrt (const Tpsa&); template Tpsa log (const Tpsa&); @@ -91,6 +92,7 @@ template Tpsa D (c template class Tpsa { + friend Tpsa pow<> (const Tpsa&, const int); friend Tpsa abs<> (const Tpsa&); friend Tpsa sqrt<> (const Tpsa&); friend Tpsa log<> (const Tpsa&); @@ -477,14 +479,23 @@ template Tpsa operator / (const T& o1, const Tpsa& o2) { return o2.inverse() * o1; } +template +Tpsa pow(const Tpsa& a_, const int n) { + Tpsa r; + Tpsa x(a_); x.c[0] = 0; x /= a_.c[0]; + int n_ = N < n ? N : n; + for(unsigned int i=0; i<=n_; i++) { + r *= 1 + x; + } + r *= std::pow(a_.c[0], n); + return r; +} + template Tpsa abs(const Tpsa& a_) { if (a_ >= 0) return a_; else return -a_; } -//template -//Tpsa fabs(const Tpsa& a_) { return abs(a_); } - #include template Tpsa sqrt(const Tpsa& a_) { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index dacc6b8..82a24dc 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -195,29 +195,10 @@ Status::type track_linepass ( lost_plane = Plane::xy; } } - } else if (element.vchamber == VChamberShape::kite) { - // kite-shaped vacuum chamber - if (((rx < element.hmin) or (rx > element.hmax))) { - lost_plane = Plane::xy; - status = Status::particle_lost; - } - if ((ry > get_kite_ry(element, 0, rx)) or // upper right - (ry > get_kite_ry(element, 1, rx)) or // upper left - (ry < get_kite_ry(element, 2, rx)) or // lower left - (ry < get_kite_ry(element, 3, rx))) { // lower right - lost_plane = Plane::xy; - status = Status::particle_lost; - } - } else if (element.vchamber == VChamberShape::ellipse) { - // elliptical vacuum chamber - double lx = (element.hmax - element.hmin) / 2; - double ly = (element.vmax - element.vmin) / 2; - double xc = (element.hmax + element.hmin) / 2; - double yc = (element.vmax + element.vmin) / 2; - T xn = (rx - xc)/lx; - T yn = (ry - yc)/ly; - T amplitude = xn*xn + yn*yn; - if (amplitude > 1) { + } else if ( + (element.vchamber == VChamberShape::rhombus) or + (element.vchamber == VChamberShape::ellipse)) { + if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { lost_plane = Plane::xy; status = Status::particle_lost; } @@ -248,26 +229,22 @@ Status::type track_linepass ( template -T get_kite_ry(const Element& elem, int side, const T& rx) { - double x1, y1, x2, y2; - if (side == 0) { - // std::cout << "upper-right" << std::endl; - x1 = 0; y1 = elem.vmax; - x2 = elem.hmax; y2 = 0; - } else if (side == 1) { - // std::cout << "upper-left" << std::endl; - x1 = elem.hmin; y1 = 0; - x2 = 0; y2 = elem.vmax; - } else if (side == 2) { - // std::cout << "lower-left" << std::endl; - x1 = elem.hmin; y1 = 0; - x2 = 0; y2 = elem.vmin; +T get_norm_amp_in_vchamber(const Element& elem, const T& rx, const T& ry) { + double lx = (elem.hmax - elem.hmin) / 2; + double ly = (elem.vmax - elem.vmin) / 2; + double xc = (elem.hmax + elem.hmin) / 2; + double yc = (elem.vmax + elem.vmin) / 2; + T xn = abs((rx - xc)/lx); + T yn = abs((ry - yc)/ly); + T amplitude; + if (elem.vchamber == VChamberShape::rhombus) { + amplitude = xn + yn; + } else if (elem.vchamber == VChamberShape::ellipse) { + amplitude = xn*xn + yn*yn; } else { - // std::cout << "lower-right" << std::endl; - x1 = 0; y1 = elem.vmin; - x2 = elem.hmax; y2 = 0; + amplitude = pow(xn, elem.vchamber) + pow(yn, elem.vchamber); } - return y1 + (rx - x1) * (y2 - y1) / (x2 - x1); + return amplitude; } From bca15966d1ef0d9d92d35a63cce4cf7c76e592dd Mon Sep 17 00:00:00 2001 From: ximenes Date: Thu, 31 Mar 2022 15:42:26 -0300 Subject: [PATCH 15/17] Improve handling of vchamber shapes --- VERSION | 2 +- include/trackcpp/tracking.h | 22 ++++++++++------------ 2 files changed, 11 insertions(+), 13 deletions(-) mode change 100644 => 120000 VERSION diff --git a/VERSION b/VERSION deleted file mode 100644 index f6cdf40..0000000 --- a/VERSION +++ /dev/null @@ -1 +0,0 @@ -4.7.0 diff --git a/VERSION b/VERSION new file mode 120000 index 0000000..a0036b8 --- /dev/null +++ b/VERSION @@ -0,0 +1 @@ +./python_package/trackcpp/VERSION \ No newline at end of file diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 82a24dc..bed5d21 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -156,7 +156,8 @@ Status::type track_linepass ( for(int i=0; i element.hmax))) { lost_plane = Plane::x; @@ -195,15 +200,8 @@ Status::type track_linepass ( lost_plane = Plane::xy; } } - } else if ( - (element.vchamber == VChamberShape::rhombus) or - (element.vchamber == VChamberShape::ellipse)) { - if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { - lost_plane = Plane::xy; - status = Status::particle_lost; - } - } else { - // any other shape not implemented safely signals lost particle + } else if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { + // lost in rhombus, elliptic and all finite p-norm shapes lost_plane = Plane::xy; status = Status::particle_lost; } From b271092e30dec85024bdadf21f32186f564ecb6a Mon Sep 17 00:00:00 2001 From: ximenes Date: Thu, 31 Mar 2022 15:44:55 -0300 Subject: [PATCH 16/17] Update VERSION --- python_package/trackcpp/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_package/trackcpp/VERSION b/python_package/trackcpp/VERSION index 6016e8a..f6cdf40 100644 --- a/python_package/trackcpp/VERSION +++ b/python_package/trackcpp/VERSION @@ -1 +1 @@ -4.6.0 +4.7.0 From d9bf41c547e957eaab5546f7e31e52211bc9ba76 Mon Sep 17 00:00:00 2001 From: ximenes Date: Thu, 31 Mar 2022 16:11:10 -0300 Subject: [PATCH 17/17] Turn off temporarily lint checking using clang-tidy (it is running into segmentation faults) --- .github/workflows/lint.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index 6e8a2e7..938e26d 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -18,6 +18,6 @@ jobs: make clean CONDA_PREFIX=/dummy-prefix PIP=pip3 PYTHON=python3 bear make -j$(nproc) - - name: Build using make and generate compile commands - run: | - bash ./tools/static-analyser.sh + # - name: Build using make and generate compile commands + # run: | + # bash ./tools/static-analyser.sh