154 #include "art/Framework/Core/EDProducer.h"
155 #include "art/Framework/Core/ModuleMacros.h"
156 #include "art/Framework/Principal/Event.h"
157 #include "art/Framework/Principal/Run.h"
158 #include "art/Framework/Services/Registry/ServiceHandle.h"
159 #include "art_root_io/TFileService.h"
160 #include "cetlib_except/exception.h"
161 #include "fhiclcpp/ParameterSet.h"
162 #include "messagefacility/MessageLogger/MessageLogger.h"
165 #include "nurandom/RandomUtils/NuRandomService.h"
168 #include "nusimdata/SimulationBase/MCParticle.h"
169 #include "nusimdata/SimulationBase/MCTruth.h"
175 #include "TDatabasePDG.h"
178 #include "CLHEP/Random/RandFlat.h"
179 #include "CLHEP/Random/RandGaussQ.h"
191 class MUSUN :
public art::EDProducer {
194 explicit MUSUN(fhicl::ParameterSet
const& pset);
200 void endRun(art::Run& run);
203 void SampleOne(
unsigned int i, simb::MCTruth& mct, CLHEP::HepRandomEngine& engine);
219 CLHEP::HepRandomEngine& engine);
324 : art::EDProducer{pset}
327 ,
fEngine(art::ServiceHandle<rndm::NuRandomService> {}->createEngine(*
this, pset,
"Seed"))
328 ,
fPDG{pset.get<
int>(
"PDG")}
330 ,
fInputDir{pset.get<std::string>(
"InputDir")}
336 ,
fEmin{pset.get<
double>(
"Emin")}
337 ,
fEmax{pset.get<
double>(
"Emax")}
338 ,
fThetamin{pset.get<
double>(
"Thetamin")}
339 ,
fThetamax{pset.get<
double>(
"Thetamax")}
349 ,
fT0{pset.get<
double>(
"T0")}
350 ,
fSigmaT{pset.get<
double>(
"SigmaT")}
351 ,
fTDist{pset.get<
int>(
"TDist")}
353 produces<std::vector<simb::MCTruth>>();
354 produces<sumdata::RunData, art::InRun>();
364 art::ServiceHandle<art::TFileService const>
tfs;
384 fTree = tfs->make<TTree>(
"Generator",
"analysis tree");
388 fTree->Branch(
"posX", &
x0,
"posX/D");
389 fTree->Branch(
"posY", &
y0,
"posY/D");
390 fTree->Branch(
"posZ", &
z0,
"posZ/D");
391 fTree->Branch(
"cosX", &
cx,
"cosX/D");
392 fTree->Branch(
"cosY", &
cy,
"cosY/D");
393 fTree->Branch(
"cosZ", &
cz,
"cosZ/D");
395 fTree->Branch(
"phi", &
phi,
"phi/D");
396 fTree->Branch(
"depth", &
dep,
"dep/D");
419 throw cet::exception(
"MUSUNGen")
420 <<
"\nThetamax has to be less than " << M_PI / 2 <<
", but was entered as " <<
fThetamax
421 <<
", this causes an error so leaving program now...\n\n";
423 throw cet::exception(
"MUSUNGen")
424 <<
"\nThetamin has to be more than 0, but was entered as " <<
fThetamin
425 <<
", this causes an error so leaving program now...\n\n";
427 throw cet::exception(
"MUSUNGen") <<
"\nMinimum angle is bigger than maximum angle....causes "
428 "an error so leaving program now....\n\n";
430 throw cet::exception(
"MUSUNGen")
431 <<
"\nPhimax has to be less than " << 2 * M_PI <<
", but was entered as " <<
fPhimax
432 <<
", this cause an error so leaving program now...\n\n";
434 throw cet::exception(
"MUSUNGen")
435 <<
"\nPhimin has to be more than 0, but was entered as " <<
fPhimin
436 <<
", this causes an error so leaving program now...\n\n";
438 throw cet::exception(
"MUSUNGen") <<
"\nMinimum angle is bigger than maximum angle....causes "
439 "an error so leaving program now....\n\n";
441 throw cet::exception(
"MUSUNGen")
442 <<
"\nMinimum energy is bigger than maximum energy....causes an error so leaving program "
446 art::ServiceHandle<geo::Geometry const> geo;
447 run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
460 std::cout <<
"Material - SURF rock" << std::endl;
462 std::cout <<
"Parameters for muon spectrum are from LVD best fit" << std::endl;
466 std::cout <<
"Azimuthal angle range = " <<
fPhimin <<
" - " << fPhimax <<
" degrees"
469 <<
" (cm^2 s)^(-1) or s^(-1) (for muons on the surface)" << std::endl;
505 std::cout <<
"Mean zenith angle (theta) = " <<
st /
NEvents <<
" degrees" << std::endl;
506 std::cout <<
"Mean azimuthal angle (phi)= " <<
sp /
NEvents <<
" degrees" << std::endl;
516 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
519 truth.SetOrigin(simb::kSingleParticle);
525 MF_LOG_DEBUG(
"MUSUN") << truth;
527 truthcol->push_back(truth);
529 evt.put(std::move(truthcol));
537 MUSUN::SampleOne(
unsigned int i, simb::MCTruth& mct, CLHEP::HepRandomEngine& engine)
539 CLHEP::RandFlat flat(engine);
540 CLHEP::RandGaussQ gauss(engine);
557 if (
phi >= 360.)
phi -= 360.;
562 int trackid = -1 * (i + 1);
563 std::string primary(
"primary");
569 double pdgfire = flat.fire();
572 static TDatabasePDG pdgt;
573 TParticlePDG* pdgp = pdgt.GetParticle(
PdgCode);
574 if (pdgp) m = pdgp->Mass();
599 double ss = sh1 + sv1 + sv2;
600 double xfl1 = flat.fire();
601 if (xfl1 <= sh1 / ss) {
606 else if (xfl1 <= (sh1 + sv1) / ss) {
636 simb::MCParticle part(trackid,
PdgCode, primary);
637 part.AddTrajectoryPoint(pos, pvec);
695 int lineNumber = 0, index = 0;
696 char inputLine[10000];
697 std::string fROOTfile;
699 for (
int a = 0;
a < 121; ++
a)
700 for (
int b = 0; b < 62; ++b)
701 for (
int c = 0; c < 50; ++c)
703 for (
int a = 0;
a < 23401; ++
a)
705 for (
int a = 0;
a < 360; ++
a)
706 for (
int b = 0; b < 91; ++b)
708 for (
int a = 0;
a < 360; ++
a)
709 for (
int b = 0; b < 91; ++b)
712 std::ostringstream File1LocStream;
714 std::string File1Loc = File1LocStream.str();
715 cet::search_path sp1(
"FW_SEARCH_PATH");
716 if (sp1.find_file(fInputFile1, fROOTfile)) File1Loc = fROOTfile;
719 throw cet::exception(
"MUSUNGen")
720 <<
"\nFile1 " << fInputFile1 <<
" not found in FW_SEARCH_PATH or at " <<
fInputDir
723 while (file1.good()) {
725 file1.getline(inputLine, 9999);
727 token = strtok(inputLine,
" ");
728 while (token != NULL) {
730 fmu[index][lineNumber] = atof(token);
731 token = strtok(NULL,
" ");
742 std::ostringstream File2LocStream;
744 std::string File2Loc = File2LocStream.str();
745 cet::search_path sp2(
"FW_SEARCH_PATH");
746 if (sp2.find_file(fInputFile2, fROOTfile)) File2Loc = fROOTfile;
747 std::ifstream file2(File2Loc.c_str(), std::ios::binary |
std::ios::in);
749 throw cet::exception(
"MUSUNGen")
750 <<
"\nFile2 " << fInputFile2 <<
" not found in FW_SEARCH_PATH or at " << fInputDir
753 int i1 = 0, i2 = 0, i3 = 0;
755 while (file2.good()) {
757 file2.read((
char*)(&readVal),
sizeof(
float));
758 spmu[i1][i2][i3] = readVal;
772 for (
int i = 0; i < 120; i++)
773 for (
int j = 0; j < 62; j++)
774 for (
int k = 0;
k < 51;
k++)
776 spmu[1][1][0] = 0.000853544;
779 std::ostringstream File3LocStream;
781 std::string File3Loc = File3LocStream.str();
782 cet::search_path sp3(
"FW_SEARCH_PATH");
783 if (sp3.find_file(fInputFile3, fROOTfile)) File3Loc = fROOTfile;
786 throw cet::exception(
"MUSUNGen")
787 <<
"\nFile3 " << fInputFile3 <<
" not found in FW_SEARCH_PATH or at " << fInputDir
790 lineNumber = index = 0;
791 while (file3.good()) {
793 file3.getline(inputLine, 9999);
795 token = strtok(inputLine,
" ");
796 while (token != NULL) {
798 depth[index][lineNumber] = atof(token);
799 token = strtok(NULL,
" ");
819 ph1 = M_PI / 180. * phi1;
820 ph2 = M_PI / 180. * phi2;
825 double theta = theta1;
829 while (theta < theta2 - dc / 2.) {
831 double theta0 = M_PI / 180. *
theta;
832 double cc = cos(theta0);
833 double ash = s_hor * cc;
834 double asv01 = s_ver1 * sqrt(1. - cc * cc);
835 double asv02 = s_ver2 * sqrt(1. - cc * cc);
836 int ic1 = (theta + 0.999);
838 if (ic2 > 91) ic2 = 91;
839 if (ic1 < 1) ic1 = 1;
843 while (phi < phi2 - dp / 2.) {
851 double asv1 = asv01 * fabs(cos(phi0));
852 double asv2 = asv02 * fabs(sin(phi0));
853 double asv0 = ash + asv1 + asv2;
855 if (figflag == 1) fl = asv0;
856 int ip1 = (phi + 0.999);
858 if (ip2 > 360) ip2 = 1;
859 if (ip1 < 1) ip1 = 360;
862 for (
int ii = 0; ii < 4; ii++) {
865 if (ip1 == 360 && (ii == 1 || ii == 3)) iip = -359;
866 if (
fmu[ip1 + iip - 1][ic1 + iic - 1] < 0) {
867 if (pow(10.,
fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 > 1
e-6) {
868 std::cout <<
"Looking at fmu [ " << ip1 <<
" + " << iip <<
" - 1 (" << ip1 + iip - 1
869 <<
") ] [ " << ic1 <<
" + " << iic <<
" - 1 (" << ic1 + iic - 1 <<
") ] ."
870 <<
"\nChanging sp1 from " << sp1 <<
" to "
871 << sp1 + pow(10.,
fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 <<
"..........."
872 << sp1 <<
" + 10 ^ (" <<
fmu[ip1 + iip - 1][ic1 + iic - 1] <<
") / 4 "
875 sp1 = sp1 + pow(10.,
fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4;
882 sc = sc + sp1 * fl * dp * M_PI / 180. * sin(theta0) * dc * M_PI / 180.;
889 theta = theta + dc / 2.;
893 for (
int ipc1 = 0; ipc1 < ipc; ipc1++)
905 CLHEP::HepRandomEngine& engine)
907 CLHEP::RandFlat flat(engine);
908 CLHEP::RandGaussQ gauss(engine);
910 #if 0 // this code is disabled for good
911 double xfl = flat.fire();
912 int loIndex = 0, hiIndex = 32400;
913 int i = (loIndex+hiIndex)/2;
914 bool foundIndex =
false;
915 if( xfl <
fnmu[loIndex] ) {
918 }
else if ( xfl >
fnmu[hiIndex] ) {
921 }
else if ( xfl >
fnmu[i-1] && xfl <=
fnmu[i] )
923 while( !foundIndex ) {
928 i = (loIndex + hiIndex)/2;
930 if( xfl >
fnmu[i-1] && xfl <=
fnmu[i] )
934 double xfl = flat.fire();
936 while (xfl >
fnmu[i])
939 int ic = (i - 2) / 360;
940 int ip = i - 2 - ic * 360;
943 theta =
the1 + ((double)ic + xfl);
945 phi =
ph1 + ((double)ip + xfl);
946 if (phi > 360) phi = phi - 360;
949 int ic1 = cos(M_PI / 180. * theta) * 50.;
950 if (ic1 < 0) ic1 = 0;
951 if (ic1 > 50) ic1 = 50;
952 int ip1 = dep / 200. - 16;
953 if (ip1 < 0) ip1 = 0;
954 if (ip1 > 61) ip1 = 61;
958 loIndex = 0, hiIndex = 120;
959 int j = (loIndex+hiIndex)/2;
961 if( xfl <
spmu[loIndex][ip1][ic1] ) {
964 }
else if ( xfl >
spmu[hiIndex][ip1][ic1] ) {
967 }
else if ( xfl >
spmu[j-1][ip1][ic1] && xfl <=
spmu[j][ip1][ic1] )
969 while( !foundIndex ) {
970 if( xfl <
spmu[j][ip1][ic1] )
974 j = (loIndex + hiIndex)/2;
976 if( xfl >
spmu[j-1][ip1][ic1] && xfl <=
spmu[j][ip1][ic1] )
981 while (xfl >
spmu[j][ip1][ic1])
985 double En1 = 0.05 * (j - 1);
986 double En2 = 0.05 * (j);
987 E = pow(10., En1 + (En2 - En1) * flat.fire());
996 DEFINE_ART_MODULE(
MUSUN)
MUSUN(fhicl::ParameterSet const &pset)
double fEmax
Maximum Kinetic Energy (GeV)
void beginRun(art::Run &run)
double fZmax
Maximum Z position.
double fYmax
Maximum Y position.
double fCavernAngle
Angle of the detector from the North to the East.
double fEmin
Minimum Kinetic Energy (GeV)
double fThetamin
Minimum theta.
std::string fInputDir
Input Directory.
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
double fT0
Central t position (ns) in world coordinates.
int fTDist
How to distribute t (gaus, or uniform)
int fPDG
PDG code of particles to generate.
void endRun(art::Run &run)
double fPhimax
Maximum phi.
std::string fInputFile3
Input File 3.
double fPhimin
Minimum phi.
double fChargeRatio
Charge ratio of particle / anti-particle.
fEngine(art::ServiceHandle< rndm::NuRandomService >() ->createEngine(*this, pset,"Seed"))
std::string fInputFile2
Input File 2.
int figflag
If want sampled from sphere or parallelepiped.
module to produce single or multiple specified particles in the detector
void initialization(double theta1, double theta2, double phi1, double phi2, int figflag, double s_hor, double s_ver1, double s_ver2, double &FI)
void SampleOne(unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
double fYmin
Minimum Y position.
double fZmin
Minimum Z position.
if &&[-z"$BASH_VERSION"] then echo Attempting to switch to bash bash shellSwitch exit fi &&["$1"= 'shellSwitch'] shift declare a IncludeDirectives for Dir in
double fSigmaT
Variation in t position (ns)
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
double fThetamax
Maximum theta.
double fXmin
Minimum X position.
art::ServiceHandle< art::TFileService > tfs
createEngine fChargeRatio
void produce(art::Event &evt)
void sampling(double &E, double &theta, double &phi, double &dep, CLHEP::HepRandomEngine &engine)
createEngine fRockDensity
std::string fInputFile1
Input File 1.
art framework interface to geometry description
BEGIN_PROLOG could also be cout
double fXmax
Maximum X position.
createEngine fCavernAngle