Skip to content

Commit

Permalink
Update
Browse files Browse the repository at this point in the history
  • Loading branch information
@dr4kan
dr4kan committed Apr 6, 2023
1 parent 317d2a8 commit 9d7823c
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 35 deletions.
27 changes: 24 additions & 3 deletions EcoMug.h
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/////////////////////////////////////////////////////////////////////////////////////
// EcoMug: Efficient COsmic MUon Generator //
// Copyright (C) 2021 Davide Pagano <[email protected]> //
// Copyright (C) 2023 Davide Pagano <[email protected]> //
// EcoMug is based on the following work: //
// D. Pagano, G. Bonomi, A. Donzella, A. Zenoni, G. Zumerle, N. Zurlo, //
// "EcoMug: an Efficient COsmic MUon Generator for cosmic-ray muons applications", //
Expand Down Expand Up @@ -332,6 +332,7 @@ class EcoMug {
double mHSphereRadius;
double mMaxFuncSkyCylinder;
std::array<double, 3> mHSphereCenterPosition;
bool mCustomJ;
EMRandom mRandom;
std::default_random_engine mEngineC;
std::discrete_distribution<int> mDiscDistC;
Expand All @@ -351,7 +352,7 @@ class EcoMug {
mJPrime(0.), mN(0.), mRandAccRej(0.), mPhi0(0.), mTheta0(0.), mAccepted(false),
mSkySize({{0., 0.}}), mSkyCenterPosition({{0., 0., 0.}}), mCylinderHeight(0.),
mCylinderRadius(0.), mCylinderCenterPosition({{0., 0., 0.}}), mHSphereRadius(0.),
mMaxFuncSkyCylinder(5.3176), mHSphereCenterPosition({{0., 0., 0.}}),
mMaxFuncSkyCylinder(5.3176), mHSphereCenterPosition({{0., 0., 0.}}), mCustomJ(false),
mEngineC(std::random_device{}()) {
mDiscDistC = std::discrete_distribution<int>({128, 100});
mMaxJ = {-1., -1., -1.};
Expand Down Expand Up @@ -427,6 +428,7 @@ class EcoMug {
/// momentum has to be in GeV/c and theta in radians
void SetDifferentialFlux(std::function<double(double, double)> J) {
mJ = J;
mCustomJ = true;
};
/// Set the seed for the internal PRNG (if 0 a random seed is used)
void SetSeed(std::uint64_t seed) {
Expand Down Expand Up @@ -507,7 +509,13 @@ class EcoMug {
/// Get the average rate of cosmic muons as well the error on it
/// The optional parameter npoints defines the number
/// of points to be used in the MC integration of the J'
void GetAverageGenRate(double &rate, double &error, int npoints = 1e7) {
void GetAverageGenRateAndError(double &rate, double &error, int npoints = 1e7) {
// Only works for the default flux definition
if (mCustomJ) {
rate = 0.;
error = 0;
return;
}
// 129.0827 is the integral of the J'prime for the sky in
// the full range of theta, phi and up to 3 TeV in energy
double k = mHorizontalRate/129.0827;
Expand All @@ -521,6 +529,19 @@ class EcoMug {
rate *= k;
error *= k;
};
/// Get the average rate of cosmic muons
/// The optional parameter npoints defines the number
/// of points to be used in the MC integration of the J'
double GetAverageGenRate(int npoints = 1e7) {
double rate, error;
GetAverageGenRateAndError(rate, error, npoints);
return rate;
};
/// Get the estimated corresponding to the provided statistics
double GetEstimatedTime(int nmuons) {
if (mCustomJ) return 0.;
return (nmuons/(GetGenSurfaceArea()/EMUnits::m2))/(GetAverageGenRate()/EMUnits::hertz*EMUnits::m2);
};
///////////////////////////////////////////////////////////////


Expand Down
70 changes: 38 additions & 32 deletions TestSuite.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@
#include <chrono>

using namespace std;
using namespace EMUnits;

class ErrorsUtility {
public:
Expand Down Expand Up @@ -99,26 +98,22 @@ void SuiteNo1(int number_of_events) {

EcoMug genPlane;
genPlane.SetUseSky();
genPlane.SetSkySize({{200.*cm, 200.*cm}});
genPlane.SetSkyCenterPosition({0., 0., 1.*mm});
genPlane.SetSkySize({{200.*EMUnits::cm, 200.*EMUnits::cm}});
genPlane.SetSkyCenterPosition({0., 0., 1.*EMUnits::mm});

EcoMug genHSphere;
genHSphere.SetUseHSphere();
genHSphere.SetHSphereRadius(200*cm);
genHSphere.SetHSphereRadius(200*EMUnits::cm);
genHSphere.SetHSphereCenterPosition({0., 0., 0.});

double offsetX = -0.*cm;
double offsetY = -0.*cm;
double offsetX = -0.*EMUnits::cm;
double offsetY = -0.*EMUnits::cm;

TVector3 P1 = {-50.*cm + offsetX, -50.*cm + offsetY, 0.};
TVector3 P2 = { 50.*cm + offsetX, -50.*cm + offsetY, 0.};
TVector3 P3 = { 50.*cm + offsetX, 50.*cm + offsetY, 0.};
TVector3 P1 = {-50.*EMUnits::cm + offsetX, -50.*EMUnits::cm + offsetY, 0.};
TVector3 P2 = { 50.*EMUnits::cm + offsetX, -50.*EMUnits::cm + offsetY, 0.};
TVector3 P3 = { 50.*EMUnits::cm + offsetX, 50.*EMUnits::cm + offsetY, 0.};
PlaneDet detector(P1, P2, P3);

cout << "OFFSET (x,y): " << offsetX << ", " << offsetY << endl;

cout << "Area = " << detector.GetArea()/m2 << endl;

auto n_gen_events = 0;
auto n_good_events = 0;
while (n_good_events < number_of_events) {
Expand All @@ -142,11 +137,8 @@ void SuiteNo1(int number_of_events) {
cout << "\n--- Generation from horizontal plane ---" << endl;
cout << "number of generated muons = " << n_gen_events << endl;
cout << "number of muons through the detector = " << n_good_events << endl;
cout << "# gen muons/generation surface [m2] = " << n_gen_events/(genPlane.GetGenSurfaceArea()/m2) << endl;
cout << "Estimated time [s] = " << n_gen_events/(genPlane.GetGenSurfaceArea()/m2)/genPlane.GetHorizontalRate()/hertz*m2 << endl;
// cout << "horizonthal generation surface are [m2] = " << genPlane.GetHorGenSurfaceArea()/m2 << endl;
// cout << "generation rate [Hz/m2] = " << genPlane.GetAverageGenRate()/hertz*m2 << endl;
// cout << "estimated time of data taking [s] = " << genPlane.GetEstimatedTime(n_gen_events)/s << endl;
cout << "# gen muons/generation surface [m2] = " << n_gen_events/(genPlane.GetGenSurfaceArea()/EMUnits::m2) << endl;
cout << "Estimated time [s] = " << genPlane.GetEstimatedTime(n_gen_events) << endl;

////////////////////////////////////////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -178,45 +170,59 @@ void SuiteNo1(int number_of_events) {
cout << "\n--- Generation from half-sphere ---" << endl;
cout << "number of generated muons = " << n_gen_events << endl;
cout << "number of muons through the detector = " << n_good_events << endl;
cout << "# gen muons/generation surface [m2] = " << n_gen_events/(genHSphere.GetGenSurfaceArea()/m2) << endl;
cout << "# gen muons/generation surface [m2] = " << n_gen_events/(genHSphere.GetGenSurfaceArea()/EMUnits::m2) << endl;
cout << "horizonthal to half-spherical rate = " << (n_good_events/detector.GetArea())/(n_gen_events/genHSphere.GetGenSurfaceArea()) << endl;
cout << "Estimated time [s] = " << n_gen_events/(genHSphere.GetGenSurfaceArea()/m2)/genPlane.GetHorizontalRate()/hertz*m2 << endl;
// cout << "generation rate [Hz/m2] = " << genHSphere.GetAverageGenRate()/hertz*m2 << endl;
// cout << "estimated time of data taking [s] = " << genHSphere.GetEstimatedTime(n_gen_events)/s << endl;
cout << "Estimated time [s] = " << genHSphere.GetEstimatedTime(n_gen_events) << endl;

cout << endl;
gApplication->Terminate();
};

double J(double p, double theta) {
double A = 0.14*pow(p, -2.7);
double B = 1. / (1. + 1.1*p*cos(theta)/115.);
double C = 0.054 / (1. + 1.1*p*cos(theta)/850.);
return A*(B+C);
}

void SuiteNo2(int number_of_events) {
// ----------------------------------------
// Suite No. 2
// ----------------------------------------

EcoMug genPlane;
genPlane.SetUseSky();
genPlane.SetSkySize({{200.*cm, 200.*cm}});
genPlane.SetSkyCenterPosition({0., 0., 1.*mm});
genPlane.SetSkySize({{200.*EMUnits::cm, 200.*EMUnits::cm}});
genPlane.SetSkyCenterPosition({0., 0., 1.*EMUnits::mm});

EcoMug genCylinder;
genCylinder.SetUseCylinder();
genCylinder.SetCylinderRadius(100.*cm);
genCylinder.SetCylinderHeight(10.*m);
genPlane.SetCylinderCenterPosition({0., 0., 5.*m});
genCylinder.SetCylinderRadius(100.*EMUnits::cm);
genCylinder.SetCylinderHeight(10.*EMUnits::m);
genPlane.SetCylinderCenterPosition({0., 0., 5.*EMUnits::m});

EcoMug genHSphere;
genHSphere.SetUseHSphere();
genHSphere.SetHSphereRadius(300*cm);
genHSphere.SetHSphereRadius(300*EMUnits::cm);
genHSphere.SetHSphereCenterPosition({0., 0., 0.});

double rateSky, rateCyl, rateHS, errorSky, errorCyl, errorHS;
genPlane.GetAverageGenRate(rateSky, errorSky, 1e7);
genCylinder.GetAverageGenRate(rateCyl, errorCyl, 1e7);
genHSphere.GetAverageGenRate(rateHS, errorHS, 1e7);
EcoMug genCustom;
genCustom.SetUseSky();
genCustom.SetSkySize({{200.*EMUnits::cm, 200.*EMUnits::cm}});
genCustom.SetSkyCenterPosition({0., 0., 1.*EMUnits::mm});
genCustom.SetDifferentialFlux(&J);

double rateSky, rateCyl, rateHS, rateCustom, errorSky, errorCyl, errorHS, errorCustom;
genPlane.GetAverageGenRateAndError(rateSky, errorSky, 1e7);
genCylinder.GetAverageGenRateAndError(rateCyl, errorCyl, 1e7);
genHSphere.GetAverageGenRateAndError(rateHS, errorHS, 1e7);
genCustom.GetAverageGenRateAndError(rateCustom, errorCustom, 1e7);

cout << "rate sky = " << rateSky << " +- " << errorSky << endl;
cout << "rate cylinder = " << rateCyl << " +- " << errorCyl << endl;
cout << "rate half-sphere = " << rateHS << " +- " << errorHS << endl;
cout << "rate custom J = " << rateCustom << " +- " << errorCustom << endl;
cout << "time custom J = " << genCustom.GetEstimatedTime(10000) << endl;
cout << "ratio (sky-to-cyl) = " << rateSky/rateCyl << " +- " << ErrorsUtility::ErrorRatio(rateSky, rateCyl, errorSky, errorCyl) << endl;
cout << "ratio (sky-to-hs) = " << rateSky/rateHS << " +- " << ErrorsUtility::ErrorRatio(rateSky, rateHS, errorSky, errorHS) << endl;

Expand Down

0 comments on commit 9d7823c

Please sign in to comment.