diff --git a/include/SimCore/Generators/GenieGenerator.h b/include/SimCore/Generators/GenieGenerator.h index 02a8132..8bd928c 100644 --- a/include/SimCore/Generators/GenieGenerator.h +++ b/include/SimCore/Generators/GenieGenerator.h @@ -80,6 +80,7 @@ class GenieGenerator : public simcore::PrimaryGenerator { std::vector targets_; std::vector abundances_; std::vector position_; + std::vector beam_size_; double target_thickness_; double time_; std::vector direction_; diff --git a/python/generators.py.in b/python/generators.py.in index 8e55197..6948e9b 100644 --- a/python/generators.py.in +++ b/python/generators.py.in @@ -219,6 +219,10 @@ class genie(simcfg.PrimaryGenerator) : Unit vector direction of incoming electron tune: str Name of GENIE tune to use + target_thickness : double + Thickness of target to use for generation + beam_size : list of double + uniform beam size width to use Examples -------- @@ -240,6 +244,7 @@ class genie(simcfg.PrimaryGenerator) : abundances = [], time = 0.0, position = [ 0.0, 0.0, 0.0 ], + beam_size = [ 0.0, 0.0 ], direction = [ 0.0, 0.0, 1.0 ], tune = 'default', spline_file = '', @@ -253,6 +258,7 @@ class genie(simcfg.PrimaryGenerator) : self.abundances = abundances self.time = time self.position = position + self.beam_size = beam_size self.direction = direction self.tune = tune self.spline_file = spline_file diff --git a/src/SimCore/Generators/GenieGenerator.cxx b/src/SimCore/Generators/GenieGenerator.cxx index 7ffcc6d..2b85cb8 100644 --- a/src/SimCore/Generators/GenieGenerator.cxx +++ b/src/SimCore/Generators/GenieGenerator.cxx @@ -72,6 +72,7 @@ void GenieGenerator::fillConfig(const framework::config::Parameters& p) time_ = p.getParameter("time");// * ns; position_ = p.getParameter< std::vector >("position"); // mm + beam_size_ = p.getParameter< std::vector >("beam_size"); // mm direction_ = p.getParameter< std::vector >("direction"); target_thickness_ = p.getParameter("target_thickness"); //mm @@ -119,6 +120,30 @@ bool GenieGenerator::validateConfig() target_thickness_ = std::abs(target_thickness_); } + if(beam_size_.size()!=2) + { + if (beam_size_.size()==0){ + std::cout << "beam size not set. Using zero." + << std::endl; + beam_size_.resize(2); beam_size_[0] = 0.0; beam_size_[1] = 0.0; + } + else{ + std::cout << "beam size is set, but does not have size 2." + << " " << beam_size_.size() + << std::endl; + ret=false; + } + } + else if(beam_size_[0] < 0 || beam_size_[1] < 0 ) + { + std::cout << "Beam size set as negative value? " + << "(" << beam_size_[0] << "," << beam_size_[1] << ")" + << std::endl; + std::cout << "Changing to positive." << std::endl; + beam_size_[0] = std::abs(beam_size_[0]); + beam_size_[1] = std::abs(beam_size_[1]); + } + //normalize abundances double abundance_sum=0; for( auto a : abundances_){ @@ -319,10 +344,12 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) } + auto x_pos = position_[0] + (G4Random::getTheGenerator()->flat()-0.5)*beam_size_[0]; + auto y_pos = position_[1] + (G4Random::getTheGenerator()->flat()-0.5)*beam_size_[1]; auto z_pos = position_[2] + (G4Random::getTheGenerator()->flat()-0.5)*target_thickness_; if(verbosity_>=1) std::cout << "Generating interaction at (x,y,z)=" - << "(" << position_[0] << "," << position_[1] << "," << z_pos << ")" + << "(" << x_pos << "," << y_pos << "," << z_pos << ")" << std::endl; genie::InitialState initial_state(targets_.at(nucl_target_i),11); @@ -332,13 +359,17 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) //setup the initial election TParticle initial_e; initial_e.SetPdgCode(11); - auto elec_i_p = std::sqrt(energy_*energy_ - initial_e.GetMass()*initial_e.GetMass()); + double elec_i_p = std::sqrt(energy_*energy_ - (double)initial_e.GetMass()*(double)initial_e.GetMass()); initial_e.SetMomentum(elec_i_p*direction_[0], elec_i_p*direction_[1], elec_i_p*direction_[2], energy_); TLorentzVector e_p4; initial_e.Momentum(e_p4); + if(verbosity_>=1) + std::cout << "Generating interation with (px,py,pz,e)=" + << "(" << e_p4.Px() << "," << e_p4.Py() << "," << e_p4.Pz() << "," << e_p4.E() << ")" + << std::endl; //calculate total xsec if( n_events_by_target_[nucl_target_i]==0) @@ -359,8 +390,8 @@ void GenieGenerator::GeneratePrimaryVertex(G4Event* event) //setup the primary vertex now G4PrimaryVertex* vertex = new G4PrimaryVertex(); - vertex->SetPosition(position_[0], - position_[1], + vertex->SetPosition(x_pos, + y_pos, z_pos); vertex->SetWeight(genie_event->Weight()); diff --git a/src/SimCore/Simulator.cxx b/src/SimCore/Simulator.cxx index 22332f4..bab2183 100644 --- a/src/SimCore/Simulator.cxx +++ b/src/SimCore/Simulator.cxx @@ -155,10 +155,10 @@ void Simulator::FillGHepParticles(const genie::EventRecord* record, my_particle.fFirstDaugher = particle->FirstDaughter(); my_particle.fLastDaughter = particle->LastDaughter(); - my_particle.fP_x = particle->P4()->X(); - my_particle.fP_y = particle->P4()->Y(); - my_particle.fP_z = particle->P4()->Z(); - my_particle.fP_t = particle->P4()->T(); + my_particle.fP_x = particle->P4()->Px(); + my_particle.fP_y = particle->P4()->Py(); + my_particle.fP_z = particle->P4()->Pz(); + my_particle.fP_t = particle->P4()->E(); my_particle.fX_x = particle->X4()->X(); my_particle.fX_y = particle->X4()->Y();