diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 69b7e92f0f8de8728607bce95ecd2770729b76a7..35115b31b03e6053f7d04d3236078412c003f110 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -18,8 +18,8 @@ #include <corsika/random/RNGManager.h> -#include <corsika/cascade/sibyll2.3c.h> #include <corsika/cascade/SibStack.h> +#include <corsika/cascade/sibyll2.3c.h> #include <corsika/process/sibyll/ParticleConversion.h> #include <corsika/units/PhysicalUnits.h> @@ -29,8 +29,8 @@ using namespace corsika::units; using namespace corsika::particles; using namespace corsika::random; -#include <typeinfo> #include <iostream> +#include <typeinfo> using namespace std; static int fCount = 0; @@ -44,44 +44,50 @@ public: // beam particles for sibyll : 1, 2, 3 for p, pi, k // read from cross section code table - int kBeam = 1; - - /* + int kBeam = 1; + + /* the target should be defined by the Environment, - ideally as full particle object so that the four momenta + ideally as full particle object so that the four momenta and the boosts can be defined.. */ // target nuclei: A < 18 // FOR NOW: assume target is oxygen int kTarget = 1; - double beamEnergy = p.GetEnergy() / 1_GeV; - std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl; - double prodCrossSection,dummy,dum1,dum2,dum3,dum4; + double beamEnergy = p.GetEnergy() / 1_GeV; + std::cout << "ProcessSplit: " + << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl; + double prodCrossSection, dummy, dum1, dum2, dum3, dum4; double dumdif[3]; - if(kTarget==1) - sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 ); + if (kTarget == 1) + sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4); else - sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy ); - - std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl; + sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy); + + std::cout << "ProcessSplit: " + << "MinStep: sibyll return: " << prodCrossSection << std::endl; CrossSectionType sig = prodCrossSection * 1_mbarn; - std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl; + std::cout << "ProcessSplit: " + << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl; const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared; - std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl; + std::cout << "ProcessSplit: " + << "nucleon mass " << nucleon_mass << std::endl; // calculate interaction length in medium - double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter ); + double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter); // pick random step lenth - std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl; + std::cout << "ProcessSplit: " + << "interaction length (g/cm2): " << int_length << std::endl; // add exponential sampling int a = 0; const double next_step = -int_length * log(s_rndm_(a)); /* what are the units of the output? slant depth or 3space length? - + */ - std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl; + std::cout << "ProcessSplit: " + << "next step (g/cm2): " << next_step << std::endl; return next_step; } @@ -93,154 +99,157 @@ public: template <typename Particle, typename Stack> void DoDiscrete(Particle& p, Stack& s) const { - cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl; - if( process::sibyll::CanInteract( p.GetPID() ) ){ - - // get energy of particle from stack - /* - stack is in GeV in lab. frame - convert to GeV in cm. frame - (assuming proton at rest as target AND - assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z) - */ - EnergyType E = p.GetEnergy(); - EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV ); + cout << "DoDiscrete: " << p.GetPID() << " interaction? " + << process::sibyll::CanInteract(p.GetPID()) << endl; + if (process::sibyll::CanInteract(p.GetPID())) { - int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() ); - - /* - the target should be defined by the Environment, - ideally as full particle object so that the four momenta - and the boosts can be defined.. - */ - // FOR NOW: set target to proton - int kTarget = 1; //p.GetPID(); - - std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl; - if (E < 8.5_GeV || Ecm < 10_GeV ) { - std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl; - p.Delete(); - fCount++; - } else { - // Sibyll does not know about units.. - double sqs = Ecm / 1_GeV; - // running sibyll, filling stack - sibyll_( kBeam, kTarget, sqs); - // running decays - //decsib_(); - // print final state - int print_unit = 6; - sib_list_( print_unit ); - - // delete current particle - p.Delete(); - - // add particles from sibyll to stack - // link to sibyll stack - SibStack ss; + // get energy of particle from stack /* - get transformation between Stack-frame and SibStack-frame - for EAS Stack-frame is lab. frame, could be different for CRMC-mode - the transformation should be derived from the input momenta - in general transformation is rotation + boost + stack is in GeV in lab. frame + convert to GeV in cm. frame + (assuming proton at rest as target AND + assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z) */ - const EnergyType proton_mass = 0.93827_GeV; - const double gamma = ( E + proton_mass ) / ( Ecm ); - const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm; - - // SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll - int i = -1; - for (auto &p: ss){ - ++i; - //transform to lab. frame, primitve - const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy(); - // add to corsika stack - auto pnew = s.NewParticle(); - pnew.SetEnergy( en_lab * 1_GeV ); - pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) ); + EnergyType E = p.GetEnergy(); + EnergyType Ecm = sqrt(2. * E * 0.93827_GeV); + + int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID()); + + /* + the target should be defined by the Environment, + ideally as full particle object so that the four momenta + and the boosts can be defined.. + */ + // FOR NOW: set target to proton + int kTarget = 1; // p.GetPID(); + + std::cout << "ProcessSplit: " + << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV + << std::endl; + if (E < 8.5_GeV || Ecm < 10_GeV) { + std::cout << "ProcessSplit: " + << " DoDiscrete: dropping particle.." << std::endl; + p.Delete(); + fCount++; + } else { + // Sibyll does not know about units.. + double sqs = Ecm / 1_GeV; + // running sibyll, filling stack + sibyll_(kBeam, kTarget, sqs); + // running decays + // decsib_(); + // print final state + int print_unit = 6; + sib_list_(print_unit); + + // delete current particle + p.Delete(); + + // add particles from sibyll to stack + // link to sibyll stack + SibStack ss; + /* + get transformation between Stack-frame and SibStack-frame + for EAS Stack-frame is lab. frame, could be different for CRMC-mode + the transformation should be derived from the input momenta + in general transformation is rotation + boost + */ + const EnergyType proton_mass = 0.93827_GeV; + const double gamma = (E + proton_mass) / (Ecm); + const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm; + + // SibStack does not know about momentum yet so we need counter to access momentum + // array in Sibyll + int i = -1; + for (auto& p : ss) { + ++i; + // transform to lab. frame, primitve + const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy(); + // add to corsika stack + auto pnew = s.NewParticle(); + pnew.SetEnergy(en_lab * 1_GeV); + pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID())); + } } - } - }else + } else p.Delete(); } - - - void Init() - { + + void Init() { fCount = 0; // define reference frame? --> defines boosts between corsika stack and model stack. - + // initialize random numbers for sibyll // FOR NOW USE SIBYLL INTERNAL !!! // rnd_ini_(); - - corsika::random::RNGManager & rmng = corsika::random::RNGManager::GetInstance();; + + corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance(); + ; const std::string str_name = "s_rndm"; rmng.RegisterRandomStream(str_name); - + // // corsika::random::RNG srng; // auto srng = rmng.GetRandomStream("s_rndm"); // test random number generator - std::cout << "ProcessSplit: " << " test sequence of random numbers." << std::endl; + std::cout << "ProcessSplit: " + << " test sequence of random numbers." << std::endl; int a = 0; - for(int i=0; i<8; ++i) - std::cout << i << " " << s_rndm_(a) << std::endl; - - //initialize Sibyll + for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl; + + // initialize Sibyll sibyll_ini_(); // set particles stable / unstable // use stack to loop over particles setup::Stack ds; - ds.NewParticle().SetPID( Code::Proton ); - ds.NewParticle().SetPID( Code::Neutron ); - ds.NewParticle().SetPID( Code::PiPlus ); - ds.NewParticle().SetPID( Code::PiMinus ); - ds.NewParticle().SetPID( Code::KPlus ); - ds.NewParticle().SetPID( Code::KMinus ); - ds.NewParticle().SetPID( Code::K0Long ); - ds.NewParticle().SetPID( Code::K0Short ); - - for( auto &p: ds){ - int s_id = process::sibyll::ConvertToSibyllRaw( p.GetPID() ); + ds.NewParticle().SetPID(Code::Proton); + ds.NewParticle().SetPID(Code::Neutron); + ds.NewParticle().SetPID(Code::PiPlus); + ds.NewParticle().SetPID(Code::PiMinus); + ds.NewParticle().SetPID(Code::KPlus); + ds.NewParticle().SetPID(Code::KMinus); + ds.NewParticle().SetPID(Code::K0Long); + ds.NewParticle().SetPID(Code::K0Short); + + for (auto& p : ds) { + int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID()); // set particle stable by setting table value negative - cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" << " stable in Sibyll .." << endl; - s_csydec_.idb[ s_id ] = -s_csydec_.idb[ s_id-1 ]; + cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" + << " stable in Sibyll .." << endl; + s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1]; p.Delete(); } } - + int GetCount() { return fCount; } private: }; -double s_rndm_(int &) -{ - static corsika::random::RNG& rmng = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");; - return rmng()/(double)rmng.max(); +double s_rndm_(int&) { + static corsika::random::RNG& rmng = + corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); + ; + return rmng() / (double)rmng.max(); } - -int main(){ +int main() { stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true); ProcessSplit p1; const auto sequence = p0 + p1; setup::Stack stack; - corsika::cascade::Cascade EAS(sequence, stack); stack.Clear(); auto particle = stack.NewParticle(); EnergyType E0 = 100_GeV; particle.SetEnergy(E0); - particle.SetPID( Code::Proton ); + particle.SetPID(Code::Proton); EAS.Init(); EAS.Run(); cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl; - } diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc index ebf824c0c2d7f53b2eaa797ab7cd8fb25b9f13ca..5bacd1489ab9bff7add254ad8edc70ccada9ea43 100644 --- a/Documentation/Examples/geometry_example.cc +++ b/Documentation/Examples/geometry_example.cc @@ -9,8 +9,8 @@ * the license. */ -#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Sphere.h> #include <corsika/geometry/Vector.h> #include <corsika/units/PhysicalUnits.h> @@ -25,7 +25,8 @@ using namespace corsika::units::si; int main() { // define the root coordinate system - geometry::CoordinateSystem& root = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); + geometry::CoordinateSystem& root = + geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // another CS defined by a translation relative to the root CS CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m}); diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc index 2c22325704dcbb2e4e1960217edf5dc8ded16edb..c51afce291181a4629b4fceaf776b0bfa9fc40cd 100644 --- a/Documentation/Examples/helix_example.cc +++ b/Documentation/Examples/helix_example.cc @@ -9,9 +9,9 @@ * the license. */ -#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Helix.h> #include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Vector.h> #include <corsika/units/PhysicalUnits.h> #include <array> @@ -23,7 +23,8 @@ using namespace corsika::geometry; using namespace corsika::units::si; int main() { - geometry::CoordinateSystem& root = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); + geometry::CoordinateSystem& root = + geometry::RootCoordinateSystem::GetInstance().GetRootCS(); Point const r0(root, {0_m, 0_m, 0_m}); auto const omegaC = 2 * M_PI * 1_Hz; diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index c97e33767dee08e4098c1e1167476118f4c16ba6..e46960dc9c3d8e8c33254c2751a601ab72843a18 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -67,14 +67,14 @@ namespace corsika::cascade { std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step); corsika::process::EProcessReturn status = - fProcesseList.DoContinuous(particle, step, fStack); + fProcesseList.DoContinuous(particle, step, fStack); if (status == corsika::process::EProcessReturn::eParticleAbsorbed) { fStack.Delete(particle); // TODO: check if this is really needed } else { fProcesseList.DoDiscrete(particle, fStack); } } - + private: Tracking& fTracking; ProcessList& fProcesseList; diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h index d38bf0da9a6c021cfd63e00ab53599a360c76c05..7e9bd3bac3bc743bea8fcb0c2b0d1d401e16efae 100644 --- a/Framework/Cascade/SibStack.h +++ b/Framework/Cascade/SibStack.h @@ -1,57 +1,59 @@ #ifndef _include_sibstack_h_ #define _include_sibstack_h_ -#include <vector> #include <string> +#include <vector> -#include <corsika/stack/Stack.h> #include <corsika/cascade/sibyll2.3c.h> #include <corsika/process/sibyll/ParticleConversion.h> +#include <corsika/stack/Stack.h> using namespace std; using namespace corsika::stack; class SibStackData { - - public: + +public: void Init(); - + void Clear() { s_plist_.np = 0; } - - int GetSize() const { return s_plist_.np; } + + int GetSize() const { return s_plist_.np; } int GetCapacity() const { return 8000; } - - void SetId(const int i, const int v) { s_plist_.llist[i]=v; } - void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v; } + void SetId(const int i, const int v) { s_plist_.llist[i] = v; } + void SetEnergy(const int i, const double v) { s_plist_.p[3][i] = v; } int GetId(const int i) const { return s_plist_.llist[i]; } double GetEnergy(const int i) const { return s_plist_.p[3][i]; } - + void Copy(const int i1, const int i2) { s_plist_.llist[i1] = s_plist_.llist[i2]; s_plist_.p[3][i1] = s_plist_.p[3][i2]; } - - protected: + +protected: void IncrementSize() { s_plist_.np++; } - void DecrementSize() { if ( s_plist_.np>0) { s_plist_.np--; } } + void DecrementSize() { + if (s_plist_.np > 0) { s_plist_.np--; } + } }; - -template<typename StackIteratorInterface> +template <typename StackIteratorInterface> class ParticleInterface : public ParticleBase<StackIteratorInterface> { using ParticleBase<StackIteratorInterface>::GetStackData; using ParticleBase<StackIteratorInterface>::GetIndex; - public: + +public: void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); } double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); } void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); } - corsika::process::sibyll::SibyllCode GetPID() const { return static_cast<corsika::process::sibyll::SibyllCode> (GetStackData().GetId(GetIndex())); } - + corsika::process::sibyll::SibyllCode GetPID() const { + return static_cast<corsika::process::sibyll::SibyllCode>( + GetStackData().GetId(GetIndex())); + } }; - typedef Stack<SibStackData, ParticleInterface> SibStack; #endif diff --git a/Framework/Cascade/sibyll2.3c.h b/Framework/Cascade/sibyll2.3c.h index e66f1f74426ede7caf3bbf19ff60097ae6c2aa3c..5d2d1e69644e812246b731aa4bf276d090d79910 100644 --- a/Framework/Cascade/sibyll2.3c.h +++ b/Framework/Cascade/sibyll2.3c.h @@ -3,94 +3,88 @@ //---------------------------------------------- // C++ interface for the SIBYLL event generator //---------------------------------------------- -//wrapper +// wrapper +extern "C" { +typedef char s_name[6]; +// SIBYLL particle stack (FORTRAN COMMON) +// variables are: np : numer of particles on stack +// p : 4momentum + mass of particles on stack +// llist : id of particles on stack +extern struct { + double p[5][8000]; + int llist[8000]; + int np; +} s_plist_; -extern"C"{ +extern struct { + double cbr[223 + 16 + 12 + 8]; + int kdec[1338 + 6 * (16 + 12 + 8)]; + int lbarp[99]; + int idb[99]; +} s_csydec_; - typedef char s_name[6]; +// additional particle stack for the mother particles of unstable particles +// stable particles have entry zero +extern struct { int llist1[8000]; } s_plist1_; - // SIBYLL particle stack (FORTRAN COMMON) - // variables are: np : numer of particles on stack - // p : 4momentum + mass of particles on stack - // llist : id of particles on stack - extern struct { - double p[5][8000]; - int llist[8000]; - int np; - }s_plist_; +// tables with particle properties +// charge, strangeness and baryon number +extern struct { + int ichp[99]; + int istr[99]; + int ibar[99]; +} s_chp_; +// tables with particle properties +// mass and mass squared +extern struct { + double am[99]; + double am2[99]; +} s_mass1_; - extern struct { - double cbr[223+16+12+8]; - int kdec[1338+6*(16+12+8)]; - int lbarp[99]; - int idb[99]; - }s_csydec_; +// table with particle names +extern struct { char namp[6][99]; } s_cnam_; - - // additional particle stack for the mother particles of unstable particles - // stable particles have entry zero - extern struct { - int llist1[8000]; - }s_plist1_; +// debug info +extern struct { + int ncall; + int ndebug; + int lun; +} s_debug_; +// lund random generator setup +// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_; - // tables with particle properties - // charge, strangeness and baryon number - extern struct { - int ichp[99]; - int istr[99]; - int ibar[99]; - }s_chp_; +// sibyll main subroutine +void sibyll_(int&, int&, double&); - // tables with particle properties - // mass and mass squared - extern struct { - double am[99]; - double am2[99]; - }s_mass1_; +// subroutine to initiate sibyll +void sibyll_ini_(); - // table with particle names - extern struct {char namp[6][99];}s_cnam_; +// subroutine to SET DECAYS +void dec_ini_(); - // debug info - extern struct {int ncall; int ndebug; int lun; }s_debug_; +// subroutine to initiate random number generator +// void rnd_ini_(); - // lund random generator setup - //extern struct {int mrlu[6]; float rrlu[100]; }ludatr_; +// print event +void sib_list_(int&); +// decay routine +void decsib_(); - // sibyll main subroutine - void sibyll_(int&,int&,double&); - - // subroutine to initiate sibyll - void sibyll_ini_(); +// interaction length +// double fpni_(double&, int&); - // subroutine to SET DECAYS - void dec_ini_(); +void sib_sigma_hnuc_(int&, int&, double&, double&, double&); +void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&); - // subroutine to initiate random number generator - //void rnd_ini_(); - - // print event - void sib_list_(int&); +double s_rndm_(int&); - // decay routine - void decsib_(); - - // interaction length - //double fpni_(double&, int&); - - void sib_sigma_hnuc_(int&,int&,double&,double&,double&); - void sib_sigma_hp_(int&,double&,double&,double&,double&,double*,double&,double&); - - double s_rndm_(int&); - - // phojet random generator setup - void pho_rndin_(int&,int&,int&,int&); +// phojet random generator setup +void pho_rndin_(int&, int&, int&, int&); } #endif - diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 0cb36f524e875bdd5055193d97443c592b4d38ec..254e08a8ea2a3aab7a10af76100278e526b67e14 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -79,8 +79,7 @@ class NewtonTracking { // naja.. not yet public: void Init() {} corsika::setup::Trajectory GetTrack(Particle& p) { - corsika::geometry::Vector<SpeedType::dimension_type> v = - p.GetDirection(); + corsika::geometry::Vector<SpeedType::dimension_type> v = p.GetDirection(); corsika::geometry::Line traj(p.GetPosition(), v); { CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS(); diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h index 50cf1b07944d67d3523df8413449d9e7cd92f565..2cd1c3e490be90485164a060cefc708dbd64806e 100644 --- a/Framework/Geometry/BaseVector.h +++ b/Framework/Geometry/BaseVector.h @@ -31,8 +31,6 @@ namespace corsika::geometry { BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) : qVector(pQVector) , cs(&pCS) {} - - }; } // namespace corsika::geometry diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h index a75ff439fad6d7e24cf9e569472d031a7c4733a6..3a06dfbd8bc65b06a83a58184088d5bec42e8b03 100644 --- a/Framework/Geometry/CoordinateSystem.h +++ b/Framework/Geometry/CoordinateSystem.h @@ -22,7 +22,7 @@ typedef Eigen::Translation<double, 3> EigenTranslation; namespace corsika::geometry { class RootCoordinateSystem; - + using corsika::units::si::length_d; class CoordinateSystem { @@ -32,17 +32,17 @@ namespace corsika::geometry { CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf) : reference(&reference) , transf(transf) {} - + CoordinateSystem() : // for creating the root CS transf(EigenTransform::Identity()) {} protected: - static auto CreateCS() { return CoordinateSystem(); } - friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can creat ONE unique root CS + static auto CreateCS() { return CoordinateSystem(); } + friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can + /// creat ONE unique root CS public: - static EigenTransform GetTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2); @@ -60,9 +60,9 @@ namespace corsika::geometry { auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const { if (axis.eVector.isZero()) { - throw std::string("null-vector given as axis parameter"); + throw std::string("null-vector given as axis parameter"); } - + EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())}; return CoordinateSystem(*this, rotation); @@ -71,9 +71,9 @@ namespace corsika::geometry { auto translateAndRotate(QuantityVector<phys::units::length_d> translation, QuantityVector<phys::units::length_d> axis, double angle) { if (axis.eVector.isZero()) { - throw std::string("null-vector given as axis parameter"); + throw std::string("null-vector given as axis parameter"); } - + EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) * EigenTranslation(translation.eVector)}; diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h index 5de0ca9cf2befca18250edcf8e8c9d3b922d0001..e093a0055fb82ceaeab2c656227b2c005f15958e 100644 --- a/Framework/Geometry/Line.h +++ b/Framework/Geometry/Line.h @@ -19,7 +19,7 @@ namespace corsika::geometry { class Line { - + using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>; Point const r0; @@ -30,12 +30,10 @@ namespace corsika::geometry { : r0(pR0) , v0(pV0) {} - Point GetPosition(corsika::units::si::TimeType t) const { - return r0 + v0 * t; - } - + Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; } + LengthType GetDistanceBetween(corsika::units::si::TimeType t1, - corsika::units::si::TimeType t2) const { + corsika::units::si::TimeType t2) const { // assert(t2 >= t1); return v0.norm() * (t2 - t1); } diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h index 9570c8808008de9a4490c2cc4c047bb649d7d06c..c9c88d7bf599e3f41d33bd92dfe746706f5e5339 100644 --- a/Framework/Geometry/Point.h +++ b/Framework/Geometry/Point.h @@ -34,7 +34,8 @@ namespace corsika::geometry { Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z) : BaseVector<phys::units::length_d>(cs, {x, y, z}) {} - // TODO: this should be private or protected, we don NOT want to expose numbers without reference to outside: + // TODO: this should be private or protected, we don NOT want to expose numbers + // without reference to outside: auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; } /// this always returns a QuantityVector as triple diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h index 58cc20c61ccdbb48ec0441017919b76b40241c67..4c3133becd9c59e4b50a9730b893a933657c5518 100644 --- a/Framework/Geometry/RootCoordinateSystem.h +++ b/Framework/Geometry/RootCoordinateSystem.h @@ -14,23 +14,20 @@ namespace corsika::geometry { class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> { - + friend class corsika::utl::Singleton<RootCoordinateSystem>; - + protected: - RootCoordinateSystem() {} - - public: + public: corsika::geometry::CoordinateSystem& GetRootCS() { return fRootCS; } const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; } private: corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT - }; - -} + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h index 99b02b377c65cb288c4a345925e9e23e775ead4f..686799756504138e349aacfa428b325166ef1e23 100644 --- a/Framework/Geometry/Trajectory.h +++ b/Framework/Geometry/Trajectory.h @@ -22,33 +22,29 @@ namespace corsika::geometry { template <typename T> class Trajectory : public T { - corsika::units::si::TimeType fTimeLength; + corsika::units::si::TimeType fTimeLength; public: - using T::GetPosition; using T::GetDistanceBetween; + using T::GetPosition; - Trajectory(T const& theT, - corsika::units::si::TimeType timeLength) + Trajectory(T const& theT, corsika::units::si::TimeType timeLength) : T(theT) , fTimeLength(timeLength) {} /*Point GetPosition(corsika::units::si::TimeType t) const { return fTraj.GetPosition(t + fTStart); }*/ - - Point GetPosition(const double u) const { - return T::GetPosition(fTimeLength * u); - } + + Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); } TimeType GetDuration() const { return fTimeLength; } - + LengthType GetDistance(const corsika::units::si::TimeType t) const { - assert(t>fTimeLength); - assert(t>=0*corsika::units::si::second); + assert(t > fTimeLength); + assert(t >= 0 * corsika::units::si::second); return T::DistanceBetween(0, t); } - }; } // namespace corsika::geometry diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc index c582b2d6bd94b66e9e873a294aac7435492dcf5b..9e3be3310bbaf3ed6c9ae22bb6faa67701f20687 100644 --- a/Framework/Geometry/testGeometry.cc +++ b/Framework/Geometry/testGeometry.cc @@ -14,10 +14,10 @@ #include <catch2/catch.hpp> #include <corsika/geometry/CoordinateSystem.h> -#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Helix.h> #include <corsika/geometry/Line.h> #include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Sphere.h> #include <corsika/geometry/Trajectory.h> #include <corsika/units/PhysicalUnits.h> @@ -147,10 +147,10 @@ TEST_CASE("Trajectories") { {1_m / second, 0_m / second, 0_m / second}); Line const line(r0, v0); - CHECK((line.GetPosition(2_s).GetCoordinates() - - QuantityVector<length_d>(2_m, 0_m, 0_m)) - .norm() - .magnitude() == Approx(0).margin(absMargin)); + CHECK( + (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m)) + .norm() + .magnitude() == Approx(0).margin(absMargin)); Trajectory<Line> base(line, 1_s); CHECK(line.GetPosition(2_s).GetCoordinates() == diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index f4e4eb80536410707c95f0715c65b21711e5d359..a07e3ef84c91e6e4e54743c4a4e87c09d933bf2f 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -162,7 +162,7 @@ namespace corsika::process { } template <typename Particle> - inline void MinStepLength(Particle& p, Trajectory& step) const { + inline void MinStepLength(Particle& p, Trajectory& step) const { A.MinStepLength(p, step); B.MinStepLength(p, step); } diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index 5e76808331194a25f5df18fcb29257898d6322aa..9757298cdc2269a90fd875cdedfa68a16f19ab5a 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -27,10 +27,8 @@ using namespace corsika::units::si; using namespace std; using namespace corsika::process; - static const int nData = 10; - class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> { public: ContinuousProcess1() {} @@ -151,16 +149,17 @@ TEST_CASE("Cascade", "[Cascade]") { DummyData p; DummyStack s; - + cout << "-->init" << endl; sequence2.Init(); cout << "-->docont" << endl; - //auto const root = corsika::geometry::CoordinateSystem::CreateRootCS(); - //corsika::geometry::Point pos(root, {0_m, 0_m, 0_m}); - //corsika::geometry::Vector<SpeedType::dimension_type> vec(root, {1_m/1_s,0_m/1_s,0_m/1_s}); - //corsika::geometry::Line traj(pos, vec); - Trajectory t;//(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns)); + // auto const root = corsika::geometry::CoordinateSystem::CreateRootCS(); + // corsika::geometry::Point pos(root, {0_m, 0_m, 0_m}); + // corsika::geometry::Vector<SpeedType::dimension_type> vec(root, + // {1_m/1_s,0_m/1_s,0_m/1_s}); corsika::geometry::Line traj(pos, vec); + Trajectory + t; //(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns)); sequence2.DoContinuous(p, t, s); cout << "-->dodisc" << endl; diff --git a/Framework/Random/RNGManager.h b/Framework/Random/RNGManager.h index c9fab9ef70e45dfc737c5cd2e4ed5dcad3146212..c440e7aaae0beb3ae745562fa4949eff5ff504ee 100644 --- a/Framework/Random/RNGManager.h +++ b/Framework/Random/RNGManager.h @@ -30,13 +30,13 @@ namespace corsika::random { class RNGManager : public corsika::utl::Singleton<RNGManager> { friend class corsika::utl::Singleton<RNGManager>; - + std::map<std::string, RNG> rngs; - std::map<std::string, std::seed_seq> seeds; + std::map<std::string, std::seed_seq> seeds; protected: RNGManager() {} - + public: /*! * This function is to be called by a module requiring a random-number diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h index 65a98e482f7727af4176ded2c2b84339e6fcdf6e..94f98ac2ccf6a0f9c96f06e7b4429b9c022e1e19 100644 --- a/Framework/Units/PhysicalConstants.h +++ b/Framework/Units/PhysicalConstants.h @@ -56,8 +56,7 @@ namespace corsika::units::si::constants { // barn moved to PhysicalUnits // constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter}; - - + // etc. } // namespace corsika::units::si::constants diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h index f1a1b845ab854284ed6c348f0c0df40b55fea4b7..e97003eee0a0e649021e20c420776afd6e586e59 100644 --- a/Framework/Units/PhysicalUnits.h +++ b/Framework/Units/PhysicalUnits.h @@ -16,7 +16,7 @@ namespace corsika::units::si { using namespace phys::units; using namespace phys::units::literals; // namespace literals = phys::units::literals; - + /// defining momentum you suckers /// dimensions, i.e. composition in base SI dimensions using momentum_d = phys::units::dimensions<1, 1, -1>; @@ -35,14 +35,12 @@ namespace corsika::units::si { using ElectricChargeType = phys::units::quantity<phys::units::electric_charge_d, double>; using EnergyType = phys::units::quantity<phys::units::energy_d, double>; - using MassType = phys::units::quantity<phys::units::mass_d, double>; + using MassType = phys::units::quantity<phys::units::mass_d, double>; using MomentumType = phys::units::quantity<momentum_d, double>; using CrossSectionType = phys::units::quantity<sigma_d, double>; } // end namespace corsika::units::si - - /** * @file PhysicalUnits * diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc index 712146c6b657fff07000ac891c71d29f632990ad..5785cc4aef1b285c3ac0676e692bbeaecdfa8d17 100644 --- a/Framework/Units/testUnits.cc +++ b/Framework/Units/testUnits.cc @@ -32,7 +32,7 @@ TEST_CASE("PhysicalUnits", "[Units]") { LengthType l1 = 10_nm; l1 = l1; - + LengthType arr0[5]; arr0[0] = 5_m; @@ -58,7 +58,7 @@ TEST_CASE("PhysicalUnits", "[Units]") { REQUIRE(1_K / 1_zK == Approx(1e21)); REQUIRE(1_K / 1_yK == Approx(1e24)); REQUIRE(1_barn / 1_mbarn == Approx(1e3)); - + REQUIRE(1_A / 1_hA == Approx(1e-2)); REQUIRE(1_m / 1_km == Approx(1e-3)); REQUIRE(1_m / 1_Mm == Approx(1e-6)); @@ -82,8 +82,8 @@ TEST_CASE("PhysicalUnits", "[Units]") { const MassType m = 1_kg; const SpeedType v = 1_m / 1_s; - REQUIRE( m*v == 1_newton_second); - + REQUIRE(m * v == 1_newton_second); + const double lgE = log10(E2 / 1_GeV); REQUIRE(lgE == Approx(log10(40.))); @@ -95,7 +95,6 @@ TEST_CASE("PhysicalUnits", "[Units]") { const LengthType farAway = std::numeric_limits<double>::infinity() * meter; REQUIRE(farAway > 100000_m); - REQUIRE_FALSE(farAway < 1e19*meter); - + REQUIRE_FALSE(farAway < 1e19 * meter); } } diff --git a/Framework/Utilities/Bit.h b/Framework/Utilities/Bit.h index 4c21449d710d04aa21cbb6082c4a1c7615796e97..67473b37a131e2baff58cf9ca5f4b0d744a84a18 100644 --- a/Framework/Utilities/Bit.h +++ b/Framework/Utilities/Bit.h @@ -14,27 +14,27 @@ // #include <utl/AugerException.h> - namespace corsika::utl { namespace Bit { - template<typename T> + template <typename T> class Array { public: - Array(T& target) : fTarget(target) { } + Array(T& target) + : fTarget(target) {} class Bit { public: - Bit(T& target, T mask) : fTarget(target), fMask(mask) { } + Bit(T& target, T mask) + : fTarget(target) + , fMask(mask) {} operator bool() const { return fTarget & fMask; } bool operator~() const { return !bool(*this); } - Bit& - operator=(const bool value) - { + Bit& operator=(const bool value) { if (value) fTarget |= fMask; else @@ -49,41 +49,38 @@ namespace corsika::utl { T fMask; }; - Bit operator[](unsigned int position) - { return Bit(fTarget, T(1) << position); } + Bit operator[](unsigned int position) { return Bit(fTarget, T(1) << position); } - Bit - At(unsigned int position) - { - if (position >= 8*sizeof(T)) - //throw std::exceptionOutOfBoundException("Running out of bits."); - throw std::exception("Running out of bits."); + Bit At(unsigned int position) { + if (position >= 8 * sizeof(T)) + // throw std::exceptionOutOfBoundException("Running out of bits."); + throw std::exception("Running out of bits."); return (*this)[position]; } - template<typename M> - Array& Mask(const M mask, const bool value) - { Bit(fTarget, mask) = value; return *this; } + template <typename M> + Array& Mask(const M mask, const bool value) { + Bit(fTarget, mask) = value; + return *this; + } - template<typename M> - T Get(const M mask) { return fTarget & T(mask); } + template <typename M> + T Get(const M mask) { + return fTarget & T(mask); + } private: T& fTarget; }; - } + } // namespace Bit // helper - template<typename T> - inline - Bit::Array<T> - AsBitArray(T& target) - { + template <typename T> + inline Bit::Array<T> AsBitArray(T& target) { return Bit::Array<T>(target); } -} - +} // namespace corsika::utl #endif diff --git a/Framework/Utilities/Dummy.h b/Framework/Utilities/Dummy.h index 09260f7eecc293d334aed170b614bad360e15056..dac2021b49df401d221c5160f825909440468e22 100644 --- a/Framework/Utilities/Dummy.h +++ b/Framework/Utilities/Dummy.h @@ -4,7 +4,6 @@ namespace corsika::utl { // void.... - } #endif diff --git a/Framework/Utilities/Singleton.h b/Framework/Utilities/Singleton.h index c718065a19d96ef9fb17320139ec16dba17a4cc1..e06d974cc3c7de8003991227ded2cd6df0b0ab53 100644 --- a/Framework/Utilities/Singleton.h +++ b/Framework/Utilities/Singleton.h @@ -34,25 +34,23 @@ namespace corsika::utl { * \ingroup stl */ - template<typename T> + template <typename T> class Singleton { public: - static - T& - GetInstance() -# ifdef __MAKECINT__ - ; -# else + static T& GetInstance() +#ifdef __MAKECINT__ + ; +#else { static T instance; return instance; } -# endif +#endif protected: // derived class can call ctor and dtor - Singleton() { } - ~Singleton() { } + Singleton() {} + ~Singleton() {} private: // no one should do copies @@ -63,22 +61,18 @@ namespace corsika::utl { #else /// classical Gamma singleton - template<typename T> + template <typename T> class Singleton { public: - static - T& - GetInstance() - { - if (!fgInstance) - fgInstance = new T; + static T& GetInstance() { + if (!fgInstance) fgInstance = new T; return *fgInstance; } protected: // derived class can call ctor and dtor - Singleton() { } - ~Singleton() { } + Singleton() {} + ~Singleton() {} private: // no one should do copies @@ -90,7 +84,6 @@ namespace corsika::utl { #endif - /** * \class LeakingSingleton Singleton.h utl/Singleton.h * @@ -123,21 +116,18 @@ namespace corsika::utl { * \ingroup stl */ - template<class T> + template <class T> class LeakingSingleton { public: - static - T& - GetInstance() - { + static T& GetInstance() { static T* const instance = new T; return *instance; } protected: // derived class can call ctor and dtor - LeakingSingleton() { } - ~LeakingSingleton() { } + LeakingSingleton() {} + ~LeakingSingleton() {} private: // no one should do copies @@ -145,7 +135,6 @@ namespace corsika::utl { LeakingSingleton& operator=(const LeakingSingleton&); }; -} - +} // namespace corsika::utl #endif diff --git a/Framework/Utilities/Test.h b/Framework/Utilities/Test.h index 8135219a0d71e5913721e4e904e367aac3970771..41c7cecd7a9856e2c5b68f845d22a25f64a3cf0a 100644 --- a/Framework/Utilities/Test.h +++ b/Framework/Utilities/Test.h @@ -13,73 +13,73 @@ \ingroup testing */ -#include <cmath> +#include <utl/Triple.h> #include <boost/format.hpp> #include <boost/tuple/tuple.hpp> #include <boost/tuple/tuple_comparison.hpp> #include <boost/tuple/tuple_io.hpp> -#include <utl/Triple.h> - +#include <cmath> namespace utl { //! Predicate used in STL for searching for whitespace struct IsSpace { - bool operator()(const char x) const - { return x == ' ' || x == '\r' || x == '\n' || x == '\t'; } + bool operator()(const char x) const { + return x == ' ' || x == '\r' || x == '\n' || x == '\t'; + } }; - /// Predicate for equality class Equal { public: - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return lhs == rhs; } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return lhs == rhs; + } static const char* Name() { return "equal"; } }; - /// Predicate for less class Less { public: - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return lhs < rhs; } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return lhs < rhs; + } static const char* Name() { return "less"; } }; - /// Predicate for less or equal class LessOrEqual { public: - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return lhs <= rhs; } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return lhs <= rhs; + } static const char* Name() { return "less or equal"; } }; - /// Predicate for greater class Greater { public: - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return lhs > rhs; } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return lhs > rhs; + } static const char* Name() { return "greater"; } }; - /// Predicate for greater or equal class GreaterOrEqual { public: - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return lhs >= rhs; } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return lhs >= rhs; + } static const char* Name() { return "greater or equal"; } }; @@ -90,41 +90,41 @@ namespace utl { */ class CloseTo { public: - CloseTo(const double eps = 1e-6) : fEpsilon(eps) { } + CloseTo(const double eps = 1e-6) + : fEpsilon(eps) {} - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return IsCloseTo(lhs, rhs); } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return IsCloseTo(lhs, rhs); + } - boost::format Name() const - { return boost::format("close (@%g) to") % fEpsilon; } + boost::format Name() const { return boost::format("close (@%g) to") % fEpsilon; } protected: - template<typename T> - bool IsCloseAbs(const T& lhs, const T& rhs) const - { return std::abs(double(lhs) - double(rhs)) < fEpsilon; } - - bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const - { return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon; } - - template<typename T> - bool IsCloseRel(const T& lhs, const T& rhs) const - { return 2*std::abs(double(lhs) - double(rhs)) / (std::abs(double(lhs)) + std::abs(double(rhs))) < fEpsilon; } - - bool - IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs) - const - { - return (2*sqrt(TupleDist2(lhs, rhs)) / - (sqrt(TupleDist2(lhs, utl::Triple(0,0,0))) + - sqrt(TupleDist2(rhs, utl::Triple(0,0,0))))) < fEpsilon; - } - - template<typename T> - bool - IsCloseTo(const T& lhs, const T& rhs) - const - { + template <typename T> + bool IsCloseAbs(const T& lhs, const T& rhs) const { + return std::abs(double(lhs) - double(rhs)) < fEpsilon; + } + + bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const { + return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon; + } + + template <typename T> + bool IsCloseRel(const T& lhs, const T& rhs) const { + return 2 * std::abs(double(lhs) - double(rhs)) / + (std::abs(double(lhs)) + std::abs(double(rhs))) < + fEpsilon; + } + + bool IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs) const { + return (2 * sqrt(TupleDist2(lhs, rhs)) / + (sqrt(TupleDist2(lhs, utl::Triple(0, 0, 0))) + + sqrt(TupleDist2(rhs, utl::Triple(0, 0, 0))))) < fEpsilon; + } + + template <typename T> + bool IsCloseTo(const T& lhs, const T& rhs) const { if (IsCloseAbs(lhs, rhs)) return true; else @@ -132,107 +132,106 @@ namespace utl { } // tuple distance - template<typename Head, typename Tail> - static - double - TupleDist2(const boost::tuples::cons<Head, Tail>& lhs, - const boost::tuples::cons<Head, Tail>& rhs) - { + template <typename Head, typename Tail> + static double TupleDist2(const boost::tuples::cons<Head, Tail>& lhs, + const boost::tuples::cons<Head, Tail>& rhs) { const double t = lhs.get_head() - rhs.get_head(); - return t*t + TupleDist2(lhs.get_tail(), rhs.get_tail()); + return t * t + TupleDist2(lhs.get_tail(), rhs.get_tail()); } static double TupleDist2(const boost::tuples::null_type& /*lhs*/, - const boost::tuples::null_type& /*rhs*/) - { return 0; } + const boost::tuples::null_type& /*rhs*/) { + return 0; + } double fEpsilon; }; - class CloseAbs : public CloseTo { public: - CloseAbs(const double eps = 1e-6) : CloseTo(eps) { } + CloseAbs(const double eps = 1e-6) + : CloseTo(eps) {} - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return IsCloseAbs(lhs, rhs); } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return IsCloseAbs(lhs, rhs); + } - boost::format Name() const - { return boost::format("absolutely close (@%g) to") % fEpsilon; } + boost::format Name() const { + return boost::format("absolutely close (@%g) to") % fEpsilon; + } }; - class CloseRel : public CloseTo { public: - CloseRel(const double eps = 1e-6) : CloseTo(eps) { } + CloseRel(const double eps = 1e-6) + : CloseTo(eps) {} - template<typename T> - bool operator()(const T& lhs, const T& rhs) const - { return IsCloseRel(lhs, rhs); } + template <typename T> + bool operator()(const T& lhs, const T& rhs) const { + return IsCloseRel(lhs, rhs); + } - boost::format Name() const - { return boost::format("relatively close (@%g) to") % fEpsilon; } + boost::format Name() const { + return boost::format("relatively close (@%g) to") % fEpsilon; + } }; - - template<typename Predicate> + template <typename Predicate> class Not : public Predicate { public: - Not() : Predicate() { } - - Not(const double eps) : Predicate(eps) { } - - template<typename T> - bool operator()(const T& x) const - { return !Predicate::operator()(x); } - - template<typename T, typename U> - bool operator()(const T& x, const U& y) const - { return !Predicate::operator()(x, y); } - - template<typename T, typename U, typename W> - bool operator()(const T& x, const U& y, const W& z) const - { return !Predicate::operator()(x, y, z); } - - static boost::format Name() - { return boost::format("not-%s") % Predicate().Name(); } - }; + Not() + : Predicate() {} + + Not(const double eps) + : Predicate(eps) {} + + template <typename T> + bool operator()(const T& x) const { + return !Predicate::operator()(x); + } + template <typename T, typename U> + bool operator()(const T& x, const U& y) const { + return !Predicate::operator()(x, y); + } + + template <typename T, typename U, typename W> + bool operator()(const T& x, const U& y, const W& z) const { + return !Predicate::operator()(x, y, z); + } - inline - utl::Triple - Diff(const utl::Triple& lhs, const utl::Triple& rhs) - { - return utl::Triple(lhs.get<0>() - rhs.get<0>(), - lhs.get<1>() - rhs.get<1>(), + static boost::format Name() { return boost::format("not-%s") % Predicate().Name(); } + }; + + inline utl::Triple Diff(const utl::Triple& lhs, const utl::Triple& rhs) { + return utl::Triple(lhs.get<0>() - rhs.get<0>(), lhs.get<1>() - rhs.get<1>(), lhs.get<2>() - rhs.get<2>()); } - /// Test condition by evaluating a predicate /** If the predicate evaluates to false, we print a failure message with the values of the left- and right-hand side and the name of the predicate. This information is normally not available when using the CPPUNIT_ASSERT macro. */ - template<class Predicate, typename T> - inline bool Test(const Predicate& pred, const T& lhs, const T& rhs) - { return pred(lhs, rhs); } - + template <class Predicate, typename T> + inline bool Test(const Predicate& pred, const T& lhs, const T& rhs) { + return pred(lhs, rhs); + } /// Main test function - template<class Predicate, typename T> - inline bool Test(const T& lhs, const T& rhs) - { return Test(Predicate(), lhs, rhs); } - + template <class Predicate, typename T> + inline bool Test(const T& lhs, const T& rhs) { + return Test(Predicate(), lhs, rhs); + } /// Test function for predicates that take an option - template<class Predicate, typename T, typename U> - inline bool Test(const T& lhs, const T& rhs, const U& eps) - { return Test(Predicate(eps), lhs, rhs); } - -} + template <class Predicate, typename T, typename U> + inline bool Test(const T& lhs, const T& rhs, const U& eps) { + return Test(Predicate(eps), lhs, rhs); + } +} // namespace utl #endif diff --git a/Framework/Utilities/testBit.cc b/Framework/Utilities/testBit.cc index dc4f454368d80f4388baad7b5b9e772e2b319bf5..a4d36d1bcae35a23c3baaa0fa96a0c05821c8fa5 100644 --- a/Framework/Utilities/testBit.cc +++ b/Framework/Utilities/testBit.cc @@ -10,18 +10,17 @@ */ #include <corsika/utl/Test.h> +#include <cppunit/extensions/HelperMacros.h> #include <tst/Verify.h> #include <utl/Bit.h> -#include <cppunit/extensions/HelperMacros.h> +#include <bitset> #include <cstdio> #include <iostream> -#include <bitset> using namespace tst; using namespace utl; using namespace std; - /** \ingroup testing */ @@ -34,14 +33,12 @@ class TestBit : public CppUnit::TestFixture { CPPUNIT_TEST_SUITE_END(); public: - void setUp() { } + void setUp() {} - void tearDown() { } + void tearDown() {} - void - TestGet() - { - const int size = sizeof(int)*8; + void TestGet() { + const int size = sizeof(int) * 8; const int bc2 = 12345; int b2 = bc2; bitset<size> b1(bc2); @@ -59,23 +56,18 @@ public: CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out3.str())); } - void - TestSet() - { - const int size = sizeof(int)*8; + void TestSet() { + const int size = sizeof(int) * 8; const int number = 12345; bitset<size> b1(number); int b2 = 11111; - for (int i = 0; i < size; ++i) - AsBitArray(b2)[i] = b1[i]; + for (int i = 0; i < size; ++i) AsBitArray(b2)[i] = b1[i]; CPPUNIT_ASSERT(Verify<Equal>(b2, number)); } - void - TestMask() - { + void TestMask() { const int n = (1 << 18) | (1 << 5); int m = 0; @@ -83,19 +75,16 @@ public: AsBitArray(m)[5] = true; CPPUNIT_ASSERT(Verify<Equal>(n, m)); - for (unsigned int i = 0; i < 8*sizeof(int); ++i) - AsBitArray(m)[i] = 0; + for (unsigned int i = 0; i < 8 * sizeof(int); ++i) AsBitArray(m)[i] = 0; CPPUNIT_ASSERT(Verify<Equal>(m, 0)); m = 1; AsBitArray(m).Mask(n, true); - CPPUNIT_ASSERT(Verify<Equal>(m, n+1)); + CPPUNIT_ASSERT(Verify<Equal>(m, n + 1)); AsBitArray(m).Mask(n, false); CPPUNIT_ASSERT(Verify<Equal>(m, 1)); } - }; - CPPUNIT_TEST_SUITE_REGISTRATION(TestBit); diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h index 9d7227ce97b324ba33977ce351c476cb1a2fe0c9..35bfbc617df904ba0e859fbf451ab1356092f207 100644 --- a/Processes/Sibyll/ParticleConversion.h +++ b/Processes/Sibyll/ParticleConversion.h @@ -35,21 +35,23 @@ namespace corsika::process::sibyll { SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) { //~ assert(handledBySibyll(pCode)); - return static_cast<SibyllCode>(corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]); + return static_cast<SibyllCode>( + corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]); } - + corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) { return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll]; } - int ConvertToSibyllRaw(corsika::particles::Code pCode){ - return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>( corsika::process::sibyll::ConvertToSibyll( pCode ) ); + int ConvertToSibyllRaw(corsika::particles::Code pCode) { + return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>( + corsika::process::sibyll::ConvertToSibyll(pCode)); } int GetSibyllXSCode(corsika::particles::Code pCode) { return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)]; } - + } // namespace corsika::process::sibyll #endif diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index 12917bfd91cf85d9f3f310256386f70bcfce483d..67c950e5585f831f8a95ae4f0e00f59df95c8919 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -50,7 +50,7 @@ TEST_CASE("Sibyll", "[processes]") { REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::Electron::GetCode())); REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::SigmaC0::GetCode())); } - + SECTION("cross-section type") { REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::Electron) == 0); diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index 807386f8c2c79ab72811d4bcbb9661858f69aebf..b0acefc4a86c72acba373feafc82b7434aa2b9a0 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -9,9 +9,9 @@ * the license. */ +#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/units/PhysicalUnits.h> -#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/logging/Logger.h> @@ -32,8 +32,8 @@ template <typename Stack> StackInspector<Stack>::~StackInspector() {} template <typename Stack> -process::EProcessReturn StackInspector<Stack>::DoContinuous( - Particle&, setup::Trajectory&, Stack& s) const { +process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&, + Stack& s) const { static int countStep = 0; if (!fReport) return EProcessReturn::eOk; [[maybe_unused]] int i = 0; @@ -42,22 +42,21 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous( for (auto& iterP : s) { EnergyType E = iterP.GetEnergy(); Etot += E; - geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // for printout + geometry::CoordinateSystem& rootCS = + geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // for printout auto pos = iterP.GetPosition().GetCoordinates(rootCS); - cout << "StackInspector: i=" << setw(5) << fixed << (i++) - << ", id=" << setw(30) << iterP.GetPID() - << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, " - << " pos=" << pos - << endl; + cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30) + << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, " + << " pos=" << pos << endl; } countStep++; - cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize() << " Estack=" << Etot / 1_GeV << " GeV" << endl; + cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize() + << " Estack=" << Etot / 1_GeV << " GeV" << endl; return EProcessReturn::eOk; } template <typename Stack> -void StackInspector<Stack>::MinStepLength(Particle&, - setup::Trajectory&) const { +void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const { // return 0; } diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h index 9d8d14f65d980d359f905aa3f1c3bc570c0aa52f..2469f83387c0b08044890ca1babfc4b2c454f8b6 100644 --- a/Processes/StackInspector/StackInspector.h +++ b/Processes/StackInspector/StackInspector.h @@ -16,8 +16,8 @@ #include <corsika/setup/SetupTrajectory.h> -//namespace corsika::setup { -//class BaseTrajectory; +// namespace corsika::setup { +// class BaseTrajectory; //} namespace corsika::process { diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h index bad174c4ddff353592ce866ea6a14d10d5e2ee07..f1dfd497a6bec31eeb9e4d84967a8b80f380ebb0 100644 --- a/Setup/SetupTrajectory.h +++ b/Setup/SetupTrajectory.h @@ -33,7 +33,7 @@ namespace corsika::setup { /// helper visitor to modify Particle by moving along Trajectory template <typename Particle> class ParticleUpdate { - + Particle& fP; public: @@ -48,11 +48,13 @@ namespace corsika::setup { }; /// helper visitor to modify Particle by moving along Trajectory - class GetDuration { + class GetDuration { public: - corsika::units::si::TimeType operator()(std::monostate const&) { return 0*corsika::units::si::second; } + corsika::units::si::TimeType operator()(std::monostate const&) { + return 0 * corsika::units::si::second; + } template <typename T> - corsika::units::si::TimeType operator()(T const &trajectory) { + corsika::units::si::TimeType operator()(T const& trajectory) { return trajectory.GetDuration(); } }; diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h index e560bdc08552b53775f97ba5c94c4d7c433dc868..064edffa2c4497743cfd4a5e5ef2caca6c86860a 100644 --- a/Stack/SuperStupidStack/SuperStupidStack.h +++ b/Stack/SuperStupidStack/SuperStupidStack.h @@ -16,12 +16,12 @@ #include <corsika/stack/Stack.h> #include <corsika/units/PhysicalUnits.h> -#include <corsika/geometry/RootCoordinateSystem.h> // remove #include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> // remove #include <corsika/geometry/Vector.h> -#include <vector> #include <algorithm> +#include <vector> using namespace corsika; @@ -29,49 +29,53 @@ namespace corsika::stack { namespace super_stupid { + using corsika::geometry::Point; + using corsika::geometry::Vector; using corsika::particles::Code; + using corsika::units::si::energy_d; using corsika::units::si::EnergyType; - using corsika::units::si::TimeType; - using corsika::units::si::SpeedType; - using corsika::units::si::second; - using corsika::units::si::meter; using corsika::units::si::joule; - using corsika::units::si::newton_second; - using corsika::units::si::energy_d; + using corsika::units::si::meter; using corsika::units::si::momentum_d; - using corsika::geometry::Point; - using corsika::geometry::Vector; + using corsika::units::si::newton_second; + using corsika::units::si::second; + using corsika::units::si::SpeedType; + using corsika::units::si::TimeType; - typedef Vector<momentum_d> MomentumVector; + typedef Vector<momentum_d> MomentumVector; /** * Example of a particle object on the stack. */ template <typename StackIteratorInterface> - class ParticleInterface : public ParticleBase<StackIteratorInterface> { - + class ParticleInterface : public ParticleBase<StackIteratorInterface> { + using ParticleBase<StackIteratorInterface>::GetStackData; using ParticleBase<StackIteratorInterface>::GetIndex; public: void SetPID(const Code id) { GetStackData().SetPID(GetIndex(), id); } void SetEnergy(const EnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); } - void SetMomentum(const MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); } + void SetMomentum(const MomentumVector& v) { + GetStackData().SetMomentum(GetIndex(), v); + } void SetPosition(const Point& v) { GetStackData().SetPosition(GetIndex(), v); } void SetTime(const TimeType& v) { GetStackData().SetTime(GetIndex(), v); } Code GetPID() const { return GetStackData().GetPID(GetIndex()); } EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); } - MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); } + MomentumVector GetMomentum() const { + return GetStackData().GetMomentum(GetIndex()); + } Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); } TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); } - + #warning this does not really work, nor make sense: Vector<SpeedType::dimension_type> GetDirection() const { - auto P = GetMomentum(); - return P/P.norm() * 1e10 * (units::si::meter/units::si::second); } - + auto P = GetMomentum(); + return P / P.norm() * 1e10 * (units::si::meter / units::si::second); + } }; /** @@ -91,13 +95,13 @@ namespace corsika::stack { int GetSize() const { return fDataPID.size(); } int GetCapacity() const { return fDataPID.size(); } - + void SetPID(const int i, const Code id) { fDataPID[i] = id; } void SetEnergy(const int i, const EnergyType e) { fDataE[i] = e; } void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; } void SetPosition(const int i, const Point& v) { fPosition[i] = v; } void SetTime(const int i, const TimeType& v) { fTime[i] = v; } - + Code GetPID(const int i) const { return fDataPID[i]; } EnergyType GetEnergy(const int i) const { return fDataE[i]; } MomentumVector GetMomentum(const int i) const { return fMomentum[i]; } @@ -110,20 +114,20 @@ namespace corsika::stack { void Copy(const int i1, const int i2) { fDataPID[i2] = fDataPID[i1]; fDataE[i2] = fDataE[i1]; - fMomentum[i2] = fMomentum[i1]; - fPosition[i2] = fPosition[i1]; - fTime[i2] = fTime[i1]; + fMomentum[i2] = fMomentum[i1]; + fPosition[i2] = fPosition[i1]; + fTime[i2] = fTime[i1]; } /** * Function to copy particle at location i2 in stack to i1 */ void Swap(const int i1, const int i2) { - std::swap(fDataPID[i2], fDataPID[i1]); - std::swap(fDataE[i2], fDataE[i1]); - std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!! - std::swap(fPosition[i2], fPosition[i1]); - std::swap(fTime[i2], fTime[i1]); + std::swap(fDataPID[i2], fDataPID[i1]); + std::swap(fDataE[i2], fDataE[i1]); + std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!! + std::swap(fPosition[i2], fPosition[i1]); + std::swap(fTime[i2], fTime[i1]); } protected: @@ -131,20 +135,20 @@ namespace corsika::stack { fDataPID.push_back(Code::Unknown); fDataE.push_back(0 * joule); #warning this here makes no sense: see issue #48 - geometry::CoordinateSystem& dummyCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); - fMomentum.push_back(MomentumVector(dummyCS, - {0 * newton_second, 0 * newton_second, 0 * newton_second})); - fPosition.push_back(Point(dummyCS, - {0 * meter, 0 * meter, 0 * meter})); - fTime.push_back(0 * second); + geometry::CoordinateSystem& dummyCS = + geometry::RootCoordinateSystem::GetInstance().GetRootCS(); + fMomentum.push_back(MomentumVector( + dummyCS, {0 * newton_second, 0 * newton_second, 0 * newton_second})); + fPosition.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter})); + fTime.push_back(0 * second); } void DecrementSize() { if (fDataE.size() > 0) { - fDataPID.pop_back(); + fDataPID.pop_back(); fDataE.pop_back(); - fMomentum.pop_back(); - fPosition.pop_back(); - fTime.pop_back(); + fMomentum.pop_back(); + fPosition.pop_back(); + fTime.pop_back(); } } @@ -153,7 +157,7 @@ namespace corsika::stack { std::vector<Code> fDataPID; std::vector<EnergyType> fDataE; - std::vector<MomentumVector> fMomentum; + std::vector<MomentumVector> fMomentum; std::vector<Point> fPosition; std::vector<TimeType> fTime; diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc index f40f2f1813224276b12606fea9ca49a03f9d005d..28cb891bdc8dfead5e3b04e8113d46f7bc127108 100644 --- a/Stack/SuperStupidStack/testSuperStupidStack.cc +++ b/Stack/SuperStupidStack/testSuperStupidStack.cc @@ -9,9 +9,9 @@ * the license. */ +#include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/stack/super_stupid/SuperStupidStack.h> #include <corsika/units/PhysicalUnits.h> -#include <corsika/geometry/RootCoordinateSystem.h> using namespace corsika::geometry; using namespace corsika::units::si; @@ -23,11 +23,9 @@ using namespace corsika::units::si; using namespace corsika; using namespace corsika::stack::super_stupid; - #include <iostream> using namespace std; - TEST_CASE("SuperStupidStack", "[stack]") { SECTION("read+write") { @@ -36,33 +34,33 @@ TEST_CASE("SuperStupidStack", "[stack]") { auto p = s.NewParticle(); p.SetPID(particles::Code::Electron); p.SetEnergy(1.5_GeV); - geometry::CoordinateSystem& dummyCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); - p.SetMomentum(MomentumVector(dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second})); + geometry::CoordinateSystem& dummyCS = + geometry::RootCoordinateSystem::GetInstance().GetRootCS(); + p.SetMomentum(MomentumVector( + dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second})); p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter})); p.SetTime(100_s); - + // read REQUIRE(s.GetSize() == 1); auto pout = s.GetNextParticle(); REQUIRE(pout.GetPID() == particles::Code::Electron); REQUIRE(pout.GetEnergy() == 1.5_GeV); #warning Fix the next two lines: - //REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule})); - //REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 * meter})); + // REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * + // joule})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 * + // meter})); REQUIRE(pout.GetTime() == 100_s); } - + SECTION("write+delete") { SuperStupidStack s; - for (int i=0; i<99; ++i) - s.NewParticle(); + for (int i = 0; i < 99; ++i) s.NewParticle(); REQUIRE(s.GetSize() == 99); - for (int i=0; i<99; ++i) - s.GetNextParticle().Delete(); - + for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete(); REQUIRE(s.GetSize() == 0); }