32 #include "NovaDAQConventions/DAQConventions.h" 93 double zlo,
double zhi);
97 double xlo,
double xhi,
98 double ylo,
double yhi,
99 double zlo,
double zhi,
100 const double surface_y,
105 double xlo,
double xhi,
106 double ylo,
double yhi,
107 double zlo,
double zhi,
108 double xyzout[],
double& dlout);
154 produces< std::vector<simb::MCTruth> >();
155 produces< sumdata::SubRunData, art::InSubRun >();
156 produces< sumdata::RunData, art::InRun >();
163 produces< sumdata::CosmicExposure, art::InSubRun >();
184 run.
put(std::move(runcol));
201 subrun.
put(std::move(ce));
207 subrun.
put(std::move(sd));
214 std::unique_ptr< std::vector<simb::MCTruth> > truthcol(
new std::vector<simb::MCTruth>);
224 double total_sampletime = 0.;
228 bool skippedSet =
true;
229 int total_skipped = -1;
230 while ( skippedSet ){
245 total_sampletime += sampletime;
250 this->
Enhance( truth, &skippedSet, rand_num);
262 truthcol->push_back(truth);
265 evt.
put(std::move(truthcol));
281 double nNeutron = -5.18;
282 double E0Neutron = 4;
284 double nPhoton = -2.63;
287 double nProton = -3.91;
306 double event_score = 0.001;
310 for (
int part_Idx = 0; part_Idx < nparticles; ++part_Idx){
316 const TLorentzVector& fourMomentum = particle.
Momentum();
317 double energy = fourMomentum.E();
319 double score =
pow((E0/energy),n);
320 if ( score > event_score ) event_score = score;
322 if ( event_score > 1 )
break;
326 if ( event_score>1 ) event_score = 1;
327 event_w = 1.0/event_score;
331 if ( random_number<event_score ) {
332 *skipThisSet =
false;
345 part.SetWeight(event_w);
373 if ( x1==0 && x2==0 &&
377 <<
"No Geometry Box Coordinates Set\n" 378 << __FILE__ <<
":" << __LINE__ <<
"\n";
387 double px = particle.
Px() / particle.
P();
388 double py = particle.
Py() / particle.
P();
389 double pz = particle.
Pz() / particle.
P();
391 double vx = particle.
EndX();
392 double vy = particle.
EndY();
393 double vz = particle.
EndZ();
407 double xyz[3] = { vx, vy, vz};
408 double dxyz[3] = { px, py, pz};
420 part.AddTrajectoryPoint(TLorentzVector(vx, vy, vz, particle.
T()),
455 geo->
WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
457 if ( x1==0 && x2==0 &&
461 <<
"No World Box Coordinates Set\n" 462 << __FILE__ <<
":" << __LINE__ <<
"\n";
480 double px = particle.
Px() / particle.
P();
481 double py = particle.
Py() / particle.
P();
482 double pz = particle.
Pz() / particle.
P();
484 double vx = particle.
Vx();
485 double vy = particle.
Vy();
486 double vz = particle.
Vz();
488 double xyz[3] = { vx, vy, vz};
489 double dxyz[3] = { px, py, pz};
496 TLorentzVector
pos(xyzo[0], xyzo[1], xyzo[2], particle.
T());
497 TLorentzVector mom(particle.
Px(), particle.
Py(), particle.
Pz(), particle.
E());
499 patsurface.AddTrajectoryPoint(
pos, mom);
501 truth_after_surface_projection.
Add(patsurface);
505 truth = truth_after_surface_projection;
534 if ( x1==0 && x2==0 &&
538 <<
"No Geometry Box Coordinates Set\n" 539 << __FILE__ <<
":" << __LINE__ <<
"\n";
563 double px = particle.
Px() / particle.
P();
564 double py = particle.
Py() / particle.
P();
565 double pz = particle.
Pz() / particle.
P();
567 double vx = particle.
Vx();
568 double vy = particle.
Vy();
569 double vz = particle.
Vz();
571 double xyz[3] = { vx, vy, vz};
572 double dxyz[3] = { px, py, pz};
576 this->
ProjectToBox(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo, dlout);
579 double m = particle.
Mass();
580 double etot = particle.
E();
581 double ke = etot -
m;
583 double X = 2.5 * dlout;
584 double ksi = 250000.;
588 if(KE_after_rock < 0)
continue;
589 double E_after_rock = KE_after_rock +
m;
591 double P_after_rock =
sqrt(E_after_rock*E_after_rock - m*m);
592 double Px = particle.
Px()/particle.
P() * P_after_rock;
593 double Py = particle.
Py()/particle.
P() * P_after_rock;
594 double Pz = particle.
Pz()/particle.
P() * P_after_rock;
596 TLorentzVector
pos(xyzo[0], xyzo[1], xyzo[2], particle.
T());
597 TLorentzVector mom(Px, Py, Pz, E_after_rock);
599 pafterrock.AddTrajectoryPoint(
pos, mom);
611 truth_after_bigbox_projection.
Add(pafterrock);
616 truth = truth_after_bigbox_projection;
623 std::vector<simb::MCTruth>* truthcol){
645 truthcol->push_back(ptruth);
657 double zlo,
double zhi) {
660 if(xyz[0]>=xlo && xyz[0]<=xhi &&
661 xyz[1]>=ylo && xyz[1]<=yhi &&
662 xyz[2]>=zlo && xyz[2]<=zhi)
return true;
666 double dx = 0.,
dy = 0.,
dz = 0.;
674 if( (
dy > 0 && dxyz[1] < 0) ||
675 (
dy < 0 && dxyz[1] > 0) ){
676 double dl =
fabs(
dy / dxyz[1] );
677 double x = xyz[0] + dxyz[0] * dl;
678 double z = xyz[2] + dxyz[2] * dl;
681 if( x>=xlo && x<=xhi &&
682 z>=zlo && z<=zhi)
return true;
689 if( (
dy > 0 && dxyz[1] < 0) ||
690 (
dy < 0 && dxyz[1] > 0) ){
691 double dl =
fabs(
dy / dxyz[1] );
692 double x = xyz[0] + dxyz[0] * dl;
693 double z = xyz[2] + dxyz[2] * dl;
696 if( x>=xlo && x<=xhi &&
697 z>=zlo && z<=zhi)
return true;
705 if( (
dz > 0 && dxyz[2] < 0) ||
706 (
dz < 0 && dxyz[2] > 0) ){
707 double dl =
fabs(
dz / dxyz[2] );
708 double x = xyz[0] + dxyz[0] * dl;
709 double y = xyz[1] + dxyz[1] * dl;
712 if( x>=xlo && x<=xhi &&
713 y>=ylo && y<=yhi)
return true;
720 if( (
dz > 0 && dxyz[2] < 0) ||
721 (
dz < 0 && dxyz[2] > 0) ){
722 double dl =
fabs(
dz / dxyz[2] );
723 double x = xyz[0] + dxyz[0] * dl;
724 double y = xyz[1] + dxyz[1] * dl;
727 if( x>=xlo && x<=xhi &&
728 y>=ylo && y<=yhi)
return true;
735 if( (dx > 0 && dxyz[0] < 0) ||
736 (dx < 0 && dxyz[0] > 0) ){
737 double dl =
fabs( dx / dxyz[0] );
738 double y = xyz[1] + dxyz[1] * dl;
739 double z = xyz[2] + dxyz[2] * dl;
742 if( z>=zlo && z<=zhi &&
743 y>=ylo && y<=yhi)
return true;
751 if( (dx > 0 && dxyz[0] < 0) ||
752 (dx < 0 && dxyz[0] > 0) ){
753 double dl =
fabs( dx / dxyz[0] );
754 double y = xyz[1] + dxyz[1] * dl;
755 double z = xyz[2] + dxyz[2] * dl;
758 if( z>=zlo && z<=zhi &&
759 y>=ylo && y<=yhi)
return true;
786 double zlo,
double zhi,
787 const double surface_y,
794 xyz[1] >= surface_y);
797 assert(surface_y >= ylo && surface_y <= yhi);
800 double dl = - (xyz[1] - surface_y) / dxyz[1];
803 for (
int i=0;
i<3; ++
i)
804 xyzout[
i] = xyz[
i] + dxyz[
i]*dl;
809 assert(xyzout[0]>=xlo && xyzout[0]<=xhi);
810 assert(xyzout[1]>=ylo && xyzout[1]<=yhi);
811 assert(xyzout[2]>=zlo && xyzout[2]<=zhi);
821 double zlo,
double zhi,
822 double xyzout[],
double& dlout) {
837 if( (dy > 0 && dxyz[1] < 0) ||
838 (dy < 0 && dxyz[1] > 0) ){
839 double dl =
fabs( dy / dxyz[1] );
840 double x = xyz[0] + dxyz[0] * dl;
841 double z = xyz[2] + dxyz[2] * dl;
844 if( x>=xlo && x<=xhi &&
860 if( (dy > 0 && dxyz[1] < 0) ||
861 (dy < 0 && dxyz[1] > 0) ){
862 double dl =
fabs( dy / dxyz[1] );
863 double x = xyz[0] + dxyz[0] * dl;
864 double z = xyz[2] + dxyz[2] * dl;
867 if( x>=xlo && x<=xhi &&
882 if( (dz > 0 && dxyz[2] < 0) ||
883 (dz < 0 && dxyz[2] > 0) ){
884 double dl =
fabs( dz / dxyz[2] );
885 double x = xyz[0] + dxyz[0] * dl;
886 double y = xyz[1] + dxyz[1] * dl;
889 if( x>=xlo && x<=xhi &&
903 if( (dz > 0 && dxyz[2] < 0) ||
904 (dz < 0 && dxyz[2] > 0) ){
905 double dl =
fabs( dz / dxyz[2] );
906 double x = xyz[0] + dxyz[0] * dl;
907 double y = xyz[1] + dxyz[1] * dl;
910 if( x>=xlo && x<=xhi &&
924 if( (dx > 0 && dxyz[0] < 0) ||
925 (dx < 0 && dxyz[0] > 0) ){
926 double dl =
fabs( dx / dxyz[0] );
927 double y = xyz[1] + dxyz[1] * dl;
928 double z = xyz[2] + dxyz[2] * dl;
931 if( z>=zlo && z<=zhi &&
946 if( (dx > 0 && dxyz[0] < 0) ||
947 (dx < 0 && dxyz[0] > 0) ){
948 double dl =
fabs( dx / dxyz[0] );
949 double y = xyz[1] + dxyz[1] * dl;
950 double z = xyz[2] + dxyz[2] * dl;
953 if( z>=zlo && z<=zhi &&
double E(const int i=0) const
double fExposure
Livetime this subrun, in seconds.
Interface to the CRY cosmic-ray generator.
const TLorentzVector & Position(const int i=0) const
double Py(const int i=0) const
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
double fProjectedZOffset
How far above the BigBox ND muons should be started.
void SetOrigin(simb::Origin_t origin)
Float_t x1[n_points_granero]
void ProjectToBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[], double &dlout)
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
unsigned int GetRandomNumberSeed()
double Px(const int i=0) const
base_engine_t & createEngine(seed_t seed)
bool fEnhance
Flag to enhace high energy particles.
::xsd::cxx::tree::exception< char > exception
art::ProductID put(std::unique_ptr< PROD > &&)
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
void Add(simb::MCParticle &part)
void ProjectMuonsToDetectorBigBox(simb::MCTruth &truth)
int fCycle
MC cycle generation number.
CosmicsGen(fhicl::ParameterSet const &pset)
ProductID put(std::unique_ptr< PROD > &&product)
double sd(Eigen::VectorXd x)
evgb::CRYHelper * fCRYHelp
CRY generator object.
base_engine_t & getEngine() const
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Interface to the CRY cosmic ray generator.
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
void DetectorBigBoxCut(simb::MCTruth &truth)
double P(const int i=0) const
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
fvar< T > exp(const fvar< T > &x)
double T(const int i=0) const
Near Detector in the NuMI cavern.
int fEnhancePDG
PDG of enhanced particles (neutron, proton or photon)
void ProjectToSurface(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, const double surface_y, double xyzout[])
const simb::MCParticle & GetParticle(int i) const
double Vx(const int i=0) const
void produce(art::Event &evt)
bool fGenSingleEvents
generate one cosmic per record
ProductID put(std::unique_ptr< PROD > &&)
const TLorentzVector & Momentum(const int i=0) const
void ProjectCosmicsToSurface(simb::MCTruth &truth)
double Pz(const int i=0) const
virtual void endSubRun(art::SubRun &run)
double Vz(const int i=0) const
assert(nhit_max >=nhit_nbins)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
virtual void beginSubRun(art::SubRun &run)
A module to produce cosmic ray events.
bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
virtual void beginRun(art::Run &run)
void SeparateMCTruth(simb::MCTruth &truth, std::vector< simb::MCTruth > *truthcol)
Method to take the MCTruth object and beak it into pieces for each individual Cosmic Ray...
Event generator information.
Module to generate only pions from cosmic rays.
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
void Enhance(simb::MCTruth &truth, bool *skipThisSet, double random_number)
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
std::string FileBaseName() const