diff --git a/sbncode/EventGenerator/MultiPart/CMakeLists.txt b/sbncode/EventGenerator/MultiPart/CMakeLists.txt index 752ac5ed7..47833dbee 100644 --- a/sbncode/EventGenerator/MultiPart/CMakeLists.txt +++ b/sbncode/EventGenerator/MultiPart/CMakeLists.txt @@ -1,69 +1,3 @@ -cet_build_plugin( MultiPartRain art::module - LIBRARIES - larcorealg::Geometry - larcore::Geometry_Geometry_service - lardata::Utilities - larevt::Filters - lardataobj::RawData - larevt::CalibrationDBI_IOVData - larevt::CalibrationDBI_Providers - lardataobj::RecoBase - lardata::ArtDataHelper - larsim::Simulation nug4::ParticleNavigation - lardataobj::Simulation - nusimdata::SimulationBase - lardata::Utilities - lardataobj::MCBase - larcoreobj::SummaryData - nusimdata::SimulationBase - nurandom::RandomUtils_NuRandomService_service - lardataobj::MCBase - CLHEP::CLHEP - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art_root_io::tfile_support ROOT::Core - art_root_io::TFileService_service - art::Persistency_Common - art::Persistency_Provenance - art::Utilities - canvas::canvas - ROOT::EG - ) - -cet_build_plugin( MultiPartVertex art::module - LIBRARIES - larcorealg::Geometry - larcore::Geometry_Geometry_service - lardata::Utilities - larevt::Filters - lardataobj::RawData - larevt::CalibrationDBI_IOVData - larevt::CalibrationDBI_Providers - lardataobj::RecoBase - lardata::ArtDataHelper - larsim::Simulation nug4::ParticleNavigation - lardataobj::Simulation - nusimdata::SimulationBase - lardata::Utilities - lardataobj::MCBase - larcoreobj::SummaryData - nusimdata::SimulationBase - nurandom::RandomUtils_NuRandomService_service - lardataobj::MCBase - CLHEP::CLHEP - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art_root_io::tfile_support ROOT::Core - art_root_io::TFileService_service - art::Persistency_Common - art::Persistency_Provenance - art::Utilities - canvas::canvas - ROOT::EG - ) - install_headers() install_fhicl() install_source() diff --git a/sbncode/EventGenerator/MultiPart/MultiPartRain_module.cc b/sbncode/EventGenerator/MultiPart/MultiPartRain_module.cc deleted file mode 100644 index 303746118..000000000 --- a/sbncode/EventGenerator/MultiPart/MultiPartRain_module.cc +++ /dev/null @@ -1,521 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// Class: MultiPartRain -// Module Type: producer -// File: MultiPartRain_module.cc -// -// Generated at Tue Dec 13 15:48:59 2016 by Kazuhiro Terao using artmod -// from cetpkgsupport v1_11_00. -//////////////////////////////////////////////////////////////////////// - -#include "art/Framework/Core/EDProducer.h" -#include "art/Framework/Core/ModuleMacros.h" -#include "art/Framework/Principal/Event.h" -#include "art/Framework/Principal/Handle.h" -#include "art/Framework/Principal/Run.h" -#include "art/Framework/Principal/SubRun.h" -//#include "art/Utilities/InputTag.h" -#include "fhiclcpp/ParameterSet.h" -#include "messagefacility/MessageLogger/MessageLogger.h" - -#include -#include -#include -#include "cetlib/pow.h" - -#include "CLHEP/Random/RandFlat.h" -#include "CLHEP/Random/RandGauss.h" -#include "CLHEP/Random/RandGeneral.h" - -#include "TRandom.h" -#include "nurandom/RandomUtils/NuRandomService.h" -#include "larcore/Geometry/Geometry.h" -#include "larcore/CoreUtils/ServiceUtil.h" -#include "larcoreobj/SummaryData/RunData.h" -#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" - -#include "nusimdata/SimulationBase/MCTruth.h" -#include "nusimdata/SimulationBase/MCParticle.h" - -#include "TLorentzVector.h" -#include "TDatabasePDG.h" - -struct PartGenParam { - std::vector pdg; - std::vector mass; - std::array multi; - std::array kerange; - bool use_mom; - bool direct_inward; - double weight; -}; - -class MultiPartRain; - -class MultiPartRain : public art::EDProducer { -public: - explicit MultiPartRain(fhicl::ParameterSet const & p); - // The destructor generated by the compiler is fine for classes - // without bare pointers or other resource use. - ~MultiPartRain(); - - // Plugins should not be copied or assigned. - MultiPartRain(MultiPartRain const &) = delete; - MultiPartRain(MultiPartRain &&) = delete; - MultiPartRain & operator = (MultiPartRain const &) = delete; - MultiPartRain & operator = (MultiPartRain &&) = delete; - - // Required functions. - void produce(art::Event & e) override; - - void beginRun(art::Run& run) override; - - void GenPosition(double& x, double& y, double& z); - - std::array extractDirection() const; - void GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz, - const double x, const double y, const double z); - - std::vector GenParticles() const; - -private: - - CLHEP::HepRandomEngine& fFlatEngine; - std::unique_ptr fFlatRandom; - std::unique_ptr fNormalRandom; - std::unique_ptr fCosmicAngleRandom; - - // exception thrower - void abort(const std::string msg) const; - - // array of particle info for generation - std::vector _param_v; - - // g4 time of generation - double _t0; - double _t0_sigma; - - // g4 position - std::array _xrange; - std::array _yrange; - std::array _zrange; - - // Whether to use uniform or cosmic-like angle distribution - bool _cosmic_distribution; - - // TPC array - std::vector > _tpc_v; - - // multiplicity constraint - size_t _multi_min; - size_t _multi_max; - - // verbosity flag - unsigned short _debug; -}; - -void MultiPartRain::abort(const std::string msg) const -{ - std::cerr << "\033[93m" << msg.c_str() << "\033[00m" << std::endl; - throw std::exception(); -} - -MultiPartRain::~MultiPartRain() -{ } - -MultiPartRain::MultiPartRain(fhicl::ParameterSet const & p) -: EDProducer(p) -, fFlatEngine(art::ServiceHandle()->registerAndSeedEngine( - createEngine(0, "HepJamesRandom", "GenRain"), "HepJamesRandom", "GenRain")) - //, fFlatEngine(art::ServiceHandle()->createEngine(*this, "HepJamesRandom", "Gen", p, "Seed")) -// : -// Initialize member data here. -{ - - // - // Random engine initialization - // - fFlatRandom = std::make_unique(fFlatEngine,0,1); - fNormalRandom = std::make_unique(fFlatEngine); - - // x^2 distribution for fCosmicAngleRandom - // will be used to draw cos(theta) - const int nbins(100); - double parent[nbins]; - for (size_t idx = 0; idx < nbins; ++idx) { - //parent[idx] = cet::square(cos(((float) idx)/nbins * util::pi())); - parent[idx] = cet::square(((float) idx)/nbins); - } - fCosmicAngleRandom = std::make_unique(fFlatEngine, parent, nbins); - - produces< std::vector >(); - produces< sumdata::RunData, art::InRun >(); - - _debug = p.get("DebugMode",0); - - _t0 = p.get("G4Time"); - _t0_sigma = p.get("G4TimeJitter"); - if(_t0_sigma < 0) this->abort("Cannot have a negative value for G4 time jitter"); - - _multi_min = p.get("MultiMin"); - _multi_max = p.get("MultiMax"); - _cosmic_distribution = p.get("CosmicDistribution", false); - - _tpc_v = p.get > >("TPCRange"); - auto const xrange = p.get > ("XRange"); - auto const yrange = p.get > ("YRange"); - auto const zrange = p.get > ("ZRange"); - auto const part_cfg = p.get("ParticleParameter"); - - _param_v.clear(); - auto const pdg_v = part_cfg.get > > ("PDGCode"); - auto const minmult_v = part_cfg.get > ("MinMulti"); - auto const maxmult_v = part_cfg.get > ("MaxMulti"); - auto const weight_v = part_cfg.get > ("ProbWeight"); - - auto kerange_v = part_cfg.get > > ("KERange"); - auto momrange_v = part_cfg.get > > ("MomRange"); - - if( (kerange_v.empty() && momrange_v.empty()) || - (!kerange_v.empty() && !momrange_v.empty()) ) { - this->abort("Only one of KERange or MomRange must be empty!"); - } - - bool use_mom = false; - if(kerange_v.empty()){ - kerange_v = momrange_v; - use_mom = true; - } - // sanity check - if( pdg_v.size() != kerange_v.size() || - pdg_v.size() != minmult_v.size() || - pdg_v.size() != maxmult_v.size() || - pdg_v.size() != weight_v.size() ) - this->abort("configuration parameters have incompatible lengths!"); - - // further sanity check (1 more depth for some double-array) - for(auto const& r : pdg_v ) { if( r.empty() ) this->abort("PDG code not given!"); } - for(auto const& r : kerange_v) { if( r.size()!=2) this->abort("Incompatible legnth @ KE vector!"); } - - for(size_t idx=0; idx maxmult_v[idx]) this->abort("Particle MinMulti > Particle MaxMulti!"); - if(minmult_v[idx] > _multi_max) this->abort("Particle MinMulti > overall MultiMax!"); - if(minmult_v[idx] > _multi_min) - _multi_min = minmult_v[idx]; - } - if(_multi_max < _multi_min) this->abort("Overall MultiMax <= overall MultiMin!"); - - if(!xrange.empty() && xrange.size() >2) this->abort("Incompatible legnth @ X vector!" ); - if(!yrange.empty() && yrange.size() >2) this->abort("Incompatible legnth @ Y vector!" ); - if(!zrange.empty() && zrange.size() >2) this->abort("Incompatible legnth @ Z vector!" ); - - // slight modification from mpv: define the overall volume across specified TPC IDs + range options - double xmin,xmax,ymin,ymax,zmin,zmax; - xmin = ymin = zmin = 1.e20; - xmax = ymax = zmax = -1.e20; - // Implementation of required member function here. - auto geop = lar::providerFrom(); - for(auto const& tpc_id : _tpc_v) { - assert(tpc_id.size() == 2); - size_t cid = tpc_id[0]; - size_t tid = tpc_id[1]; - auto const& cryostat = geop->Cryostat(geo::CryostatID(cid)); - assert(cryostat.HasTPC(tid)); - - auto const& tpc = cryostat.TPC(tid); - auto const& tpcabox = tpc.ActiveBoundingBox(); - xmin = std::min(tpcabox.MinX(), xmin); - ymin = std::min(tpcabox.MinY(), ymin); - zmin = std::min(tpcabox.MinZ(), zmin); - xmax = std::max(tpcabox.MaxX(), xmax); - ymax = std::max(tpcabox.MaxY(), ymax); - zmax = std::max(tpcabox.MaxZ(), zmax); - - if(_debug) { - std::cout << "Using Cryostat " << tpc_id[0] << " TPC " << tpc_id[1] - << " ... X " << xmin << " => " << xmax - << " ... Y " << ymin << " => " << ymax - << " ... Z " << zmin << " => " << zmax - << std::endl; - } - } - - // range register - if(xrange.size()==1) { _xrange[0] = _xrange[1] = xrange[0]; } - if(yrange.size()==1) { _yrange[0] = _yrange[1] = yrange[0]; } - if(zrange.size()==1) { _zrange[0] = _zrange[1] = zrange[0]; } - if(xrange.size()==2) { _xrange[0] = xrange[0]; _xrange[1] = xrange[1]; } - if(yrange.size()==2) { _yrange[0] = yrange[0]; _yrange[1] = yrange[1]; } - if(zrange.size()==2) { _zrange[0] = zrange[0]; _zrange[1] = zrange[1]; } - - _xrange[0] = xmin + _xrange[0]; - _xrange[1] = xmax - _xrange[1]; - _yrange[0] = ymin + _yrange[0]; - _yrange[1] = ymax - _yrange[1]; - _zrange[0] = zmin + _zrange[0]; - _zrange[1] = zmax - _zrange[1]; - - // check - assert(_xrange[0] <= _xrange[1]); - assert(_yrange[0] <= _yrange[1]); - assert(_zrange[0] <= _zrange[1]); - - if(_debug>0) { - std::cout<<"Particle generation world boundaries..."< " << _xrange[1] << std::endl - <<"Y " << _yrange[0] << " => " << _yrange[1] << std::endl - <<"Z " << _zrange[0] << " => " << _zrange[1] << std::endl; - } - - - // register - //auto db = new TDatabasePDG; - auto db = TDatabasePDG::Instance(); - bool direct_inward = p.get("DirectInward", true); - if (_cosmic_distribution && direct_inward) { - std::cout << "Cannot satisfy both CosmicDistribution and Inward Directioning at the same time" << std::endl; - throw std::exception(); - } - for(size_t idx=0; idxGetParticle(param.pdg[i])->Mass(); - - // sanity check - if(kerange[0]<0 || kerange[1]<0) - this->abort("You provided negative energy? Fuck off Mr. Trump."); - - // overall range check - if(param.kerange[0] > param.kerange[1]) this->abort("KE range has no phase space..."); - - if(_debug>0) { - std::cout << "Generating particle (PDG"; - for(auto const& pdg : param.pdg) std::cout << " " << pdg; - std::cout << ")" << std::endl - << (param.use_mom ? " KE range ....... " : " Mom range ....... ") - << param.kerange[0] << " => " << param.kerange[1] << " MeV" << std::endl - << std::endl; - } - - _param_v.push_back(param); - } -} - -void MultiPartRain::beginRun(art::Run& run) -{ - // grab the geometry object to see what geometry we are using - art::ServiceHandle geo; - - std::unique_ptr runData(new sumdata::RunData(geo->DetectorName())); - - run.put(std::move(runData), art::fullRun()); - - return; -} - -std::vector MultiPartRain::GenParticles() const { - - std::vector result; - std::vector gen_count_v(_param_v.size(),0); - std::vector weight_v(_param_v.size(),0); - for(size_t idx=0; idx<_param_v.size(); ++idx) - weight_v[idx] = _param_v[idx].weight; - - int num_part = (int)(fFlatRandom->fire(_multi_min,_multi_max+1-1.e-10)); - - while(num_part) { - - double total_weight = 0; - for(auto const& v : weight_v) total_weight += v; - - double rval = 0; - rval = fFlatRandom->fire(0,total_weight); - - size_t idx = 0; - for(idx = 0; idx < weight_v.size(); ++idx) { - rval -= weight_v[idx]; - if(rval <=0.) break; - } - - // register to the output - result.push_back(idx); - - // if generation count exceeds max, set probability weight to be 0 - gen_count_v[idx] += 1; - if(gen_count_v[idx] >= _param_v[idx].multi[1]) - weight_v[idx] = 0.; - - --num_part; - } - return result; -} - -void MultiPartRain::GenPosition(double& x, double& y, double& z) { - - x = fFlatRandom->fire(_xrange[0],_xrange[1]); - y = fFlatRandom->fire(_yrange[0],_yrange[1]); - z = fFlatRandom->fire(_zrange[0],_zrange[1]); - - if(_debug>0) { - std::cout << "Generating a rain particle at (" - << x << "," << y << "," << z << ")" << std::endl; - } -} - -std::array MultiPartRain::extractDirection() const { - double px, py, pz; - if (_cosmic_distribution) { - double phi = fFlatRandom->fire(0, 2 * util::pi()); - // Zenith Angle Theta is in [pi/2, pi] - double costheta = - 1. * fCosmicAngleRandom->fire(); - double sintheta = sqrt(1 - costheta * costheta); - pz = cos(phi) * sintheta; - px = sin(phi) * sintheta; - py = costheta; - } else { - px = fNormalRandom->fire(0, 1); - py = fNormalRandom->fire(0, 1); - pz = fNormalRandom->fire(0, 1); - double p = std::hypot(px, py, pz); - px = px / p; - py = py / p; - pz = pz / p; - } - - std::array result = {px, py, pz}; - return result; -} - -void MultiPartRain::GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz, - const double x, const double y, const double z) { - - double tot_energy = 0; - if(param.use_mom) - tot_energy = sqrt(cet::square(fFlatRandom->fire(param.kerange[0],param.kerange[1])) + cet::square(mass)); - else - tot_energy = fFlatRandom->fire(param.kerange[0],param.kerange[1]) + mass; - - double mom_mag = sqrt(cet::square(tot_energy) - cet::square(mass)); - - /* Generating unit vector with uniform distribution - * in direction = over the sphere. - * - * It is sufficient to draw a normal variable in - * each direction and normalize. - * - * https://mathworld.wolfram.com/SpherePointPicking.html - */ - - if(!param.direct_inward) { - std::cout<<"No inward directioning..."< p = extractDirection(); - px = p[0]; py = p[1]; pz = p[2]; - }else{ - double sign_x = ( (x - _xrange[0]) < (_xrange[1] - x) ? 1.0 : -1.0 ); - double sign_y = ( (y - _yrange[0]) < (_yrange[1] - y) ? 1.0 : -1.0 ); - double sign_z = ( (z - _zrange[0]) < (_zrange[1] - z) ? 1.0 : -1.0 ); - - if(_debug) { - std::cout<<"Generating XYZ direction sign for a particle at (" << x << "," << y << "," << z << ")" << std::endl - <<"X: " << _xrange[0] << " => " << _xrange[1] << " ... " << sign_x << std::endl - <<"Y: " << _yrange[0] << " => " << _yrange[1] << " ... " << sign_y << std::endl - <<"Z: " << _zrange[0] << " => " << _zrange[1] << " ... " << sign_z << std::endl - << std::endl; - } - - while(1) { - std::array p = extractDirection(); - px = p[0]; py = p[1]; pz = p[2]; - if( (px * sign_x) >= 0. && - (py * sign_y) >= 0. && - (pz * sign_z) >= 0. ) - break; - } - - } - //std::cout<<"LOGME,"<1) - std::cout << " Direction : (" << px << "," << py << "," << pz << ")" << std::endl - << " Momentum : " << mom_mag << " [MeV/c]" << std::endl - << " Energy : " << tot_energy << " [MeV/c^2]" << std::endl; - px *= mom_mag; - py *= mom_mag; - pz *= mom_mag; - -} - -void MultiPartRain::produce(art::Event & e) -{ - if(_debug>0) std::cout << "Processing a new event..." << std::endl; - - std::unique_ptr< std::vector > mctArray(new std::vector); - - simb::MCTruth mct; - - mct.SetOrigin(simb::kCosmicRay); - - std::vector part_v; - - auto const param_idx_v = GenParticles(); - if(_debug) - std::cout << "Generating" << param_idx_v.size() << " particles..." << std::endl; - - for(size_t idx=0; idxfire(0,param.pdg.size()-1.e-10)); - auto const& pdg = param.pdg[pdg_index]; - auto const& mass = param.mass[pdg_index]; - if(_debug) std::cout << " " << idx << "th instance PDG " << pdg << std::endl; - double x, y, z; - GenPosition(x,y,z); - double g4_time = fFlatRandom->fire(_t0 - _t0_sigma/2., _t0 + _t0_sigma/2.); - TLorentzVector pos(x,y,z,g4_time); - std::cout<<"Generating momentum..."<push_back(mct); - - e.put(std::move(mctArray)); -} - -DEFINE_ART_MODULE(MultiPartRain) diff --git a/sbncode/EventGenerator/MultiPart/MultiPartVertex_module.cc b/sbncode/EventGenerator/MultiPart/MultiPartVertex_module.cc deleted file mode 100644 index 5a4e3cf08..000000000 --- a/sbncode/EventGenerator/MultiPart/MultiPartVertex_module.cc +++ /dev/null @@ -1,709 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// Class: MultiPartVertex -// Module Type: producer -// File: MultiPartVertex_module.cc -// -// Generated at Tue Dec 13 15:48:59 2016 by Kazuhiro Terao using artmod -// from cetpkgsupport v1_11_00. -//////////////////////////////////////////////////////////////////////// - -#include "art/Framework/Core/EDProducer.h" -#include "art/Framework/Core/ModuleMacros.h" -#include "art/Framework/Principal/Event.h" -#include "art/Framework/Principal/Handle.h" -#include "art/Framework/Principal/Run.h" -#include "art/Framework/Principal/SubRun.h" -//#include "art/Utilities/InputTag.h" -#include "fhiclcpp/ParameterSet.h" -#include "messagefacility/MessageLogger/MessageLogger.h" - -#include -#include -#include "cetlib/pow.h" - -#include "CLHEP/Random/RandFlat.h" -#include "CLHEP/Random/RandGauss.h" -#include "TRandom.h" -#include "nurandom/RandomUtils/NuRandomService.h" -#include "larcore/Geometry/Geometry.h" -#include "larcore/CoreUtils/ServiceUtil.h" -#include "larcoreobj/SummaryData/RunData.h" -#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" - -#include "nusimdata/SimulationBase/MCTruth.h" -#include "nusimdata/SimulationBase/MCParticle.h" - -#include "TLorentzVector.h" -#include "TDatabasePDG.h" - -struct PartGenParam { - std::vector pdg; - std::vector mass; - std::array multi; - std::array kerange; - bool use_mom; - double weight; -}; - -class MultiPartVertex; - -class MultiPartVertex : public art::EDProducer { -public: - explicit MultiPartVertex(fhicl::ParameterSet const & p); - // The destructor generated by the compiler is fine for classes - // without bare pointers or other resource use. - ~MultiPartVertex(); - - // Plugins should not be copied or assigned. - MultiPartVertex(MultiPartVertex const &) = delete; - MultiPartVertex(MultiPartVertex &&) = delete; - MultiPartVertex & operator = (MultiPartVertex const &) = delete; - MultiPartVertex & operator = (MultiPartVertex &&) = delete; - - // Required functions. - void produce(art::Event & e) override; - - void beginRun(art::Run& run) override; - - void GenPosition(double& x, double& y, double& z); - void GenBoost(double& bx, double& by, double& bz); - - std::array extractDirection() const; - void GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz); - void GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz, bool& same_range); - - void GenMomentumSF(const double& sf, const double& mass, const double& p, double &p_sf); - std::vector GenParticles() const; - -private: - - CLHEP::HepRandomEngine& fFlatEngine; - std::unique_ptr fFlatRandom; - std::unique_ptr fNormalRandom; - - // exception thrower - void abort(const std::string msg) const; - - // array of particle info for generation - std::vector _param_v; - - // g4 time of generation - double _t0; - double _t0_sigma; - - // g4 position - std::array _xrange; - std::array _yrange; - std::array _zrange; - - // TPC range - std::vector > _tpc_v; - - // multiplicity constraint - size_t _multi_min; - size_t _multi_max; - - // Range of boosts - std::array _gamma_beta_range; - - // verbosity flag - unsigned short _debug; - - unsigned short _use_boost; - - unsigned short _revert; -}; - -void MultiPartVertex::abort(const std::string msg) const -{ - std::cerr << "\033[93m" << msg.c_str() << "\033[00m" << std::endl; - throw std::exception(); -} - -MultiPartVertex::~MultiPartVertex() -{ } - -MultiPartVertex::MultiPartVertex(fhicl::ParameterSet const & p) - : EDProducer(p) - , fFlatEngine(art::ServiceHandle()->registerAndSeedEngine( - createEngine(0, "HepJamesRandom", "GenVertex"), "HepJamesRandom", "GenVertex")) - //, fFlatEngine(art::ServiceHandle()->createEngine(*this, "HepJamesRandom", "Gen", p, "Seed")) - // Initialize member data here. -{ - - // - // Random engine initialization - // - fFlatRandom = std::make_unique(fFlatEngine,0,1); - fNormalRandom = std::make_unique(fFlatEngine); - - produces< std::vector >(); - produces< sumdata::RunData, art::InRun >(); - - _debug = p.get("DebugMode",0); - _revert = p.get("Revert",0); - _use_boost = p.get("UseBoost",0); - _t0 = p.get("G4Time"); - _t0_sigma = p.get("G4TimeJitter"); - if(_t0_sigma < 0) this->abort("Cannot have a negative value for G4 time jitter"); - - _multi_min = p.get("MultiMin"); - _multi_max = p.get("MultiMax"); - - _tpc_v = p.get > >("TPCRange"); - auto const xrange = p.get > ("XRange"); - auto const yrange = p.get > ("YRange"); - auto const zrange = p.get > ("ZRange"); - - auto const part_cfg = p.get("ParticleParameter"); - _gamma_beta_range = p.get >("GammaBetaRange", {0.0,0.0}); // _gamma_beta denotes gamma * beta - - _param_v.clear(); - auto const pdg_v = part_cfg.get > > ("PDGCode"); - auto const minmult_v = part_cfg.get > ("MinMulti"); - auto const maxmult_v = part_cfg.get > ("MaxMulti"); - auto const weight_v = part_cfg.get > ("ProbWeight"); - - auto kerange_v = part_cfg.get > > ("KERange"); - auto momrange_v = part_cfg.get > > ("MomRange"); - - if( (kerange_v.empty() && momrange_v.empty()) || - (!kerange_v.empty() && !momrange_v.empty()) ) { - this->abort("Only one of KERange or MomRange must be empty!"); - } - - bool use_mom = false; - if(kerange_v.empty()){ - kerange_v = momrange_v; - use_mom = true; - } - // sanity check - if( pdg_v.size() != kerange_v.size() || - pdg_v.size() != minmult_v.size() || - pdg_v.size() != maxmult_v.size() || - pdg_v.size() != weight_v.size() ) - this->abort("configuration parameters have incompatible lengths!"); - - // further sanity check (1 more depth for some double-array) - for(auto const& r : pdg_v ) { if( r.empty() ) this->abort("PDG code not given!"); } - for(auto const& r : kerange_v) { if( r.size()!=2) this->abort("Incompatible legnth @ KE vector!"); } - if(_gamma_beta_range[0] > _gamma_beta_range[1]) this->abort("Incompatible boost range!"); - - size_t multi_min = 0; - for(size_t idx=0; idx maxmult_v[idx]) this->abort("Particle MinMulti > Particle MaxMulti!"); - if(minmult_v[idx] > _multi_max) this->abort("Particle MinMulti > overall MultiMax!"); - multi_min += minmult_v[idx]; - } - _multi_min = std::max(_multi_min, multi_min); - if(_multi_max < _multi_min) this->abort("Overall MultiMax <= overall MultiMin!"); - - /* - if(!xrange.empty() && xrange.size() >2) this->abort("Incompatible legnth @ X vector!" ); - if(!yrange.empty() && yrange.size() >2) this->abort("Incompatible legnth @ Y vector!" ); - if(!zrange.empty() && zrange.size() >2) this->abort("Incompatible legnth @ Z vector!" ); - - // range register - if(xrange.size()==1) { _xrange[0] = _xrange[1] = xrange[0]; } - if(xrange.size()==2) { _xrange[0] = xrange[0]; _xrange[1] = xrange[1]; } - if(yrange.size()==1) { _yrange[0] = _yrange[1] = yrange[0]; } - if(yrange.size()==2) { _yrange[0] = yrange[0]; _yrange[1] = yrange[1]; } - if(zrange.size()==1) { _zrange[0] = _zrange[1] = zrange[0]; } - if(zrange.size()==2) { _zrange[0] = zrange[0]; _zrange[1] = zrange[1]; } - */ - - if(!xrange.empty() && xrange.size() >2) this->abort("Incompatible legnth @ X vector!" ); - if(!yrange.empty() && yrange.size() >2) this->abort("Incompatible legnth @ Y vector!" ); - if(!zrange.empty() && zrange.size() >2) this->abort("Incompatible legnth @ Z vector!" ); - - // slight modification from mpv: define the overall volume across specified TPC IDs + range options - double xmin,xmax,ymin,ymax,zmin,zmax; - xmin = ymin = zmin = 1.e20; - xmax = ymax = zmax = -1.e20; - // Implementation of required member function here. - auto geop = lar::providerFrom(); - for(auto const& tpc_id : _tpc_v) { - assert(tpc_id.size() == 2); - size_t cid = tpc_id[0]; - size_t tid = tpc_id[1]; - auto const& cryostat = geop->Cryostat(geo::CryostatID(cid)); - assert(cryostat.HasTPC(tid)); - - auto const& tpc = cryostat.TPC(tid); - auto const& tpcabox = tpc.ActiveBoundingBox(); - xmin = std::min(tpcabox.MinX(), xmin); - ymin = std::min(tpcabox.MinY(), ymin); - zmin = std::min(tpcabox.MinZ(), zmin); - xmax = std::max(tpcabox.MaxX(), xmax); - ymax = std::max(tpcabox.MaxY(), ymax); - zmax = std::max(tpcabox.MaxZ(), zmax); - - if(_debug) { - std::cout << "Using Cryostat " << tpc_id[0] << " TPC " << tpc_id[1] - << " ... X " << xmin << " => " << xmax - << " ... Y " << ymin << " => " << ymax - << " ... Z " << zmin << " => " << zmax - << std::endl; - } - } - - // range register - if(xrange.size()==1) { _xrange[0] = _xrange[1] = xrange[0]; } - if(yrange.size()==1) { _yrange[0] = _yrange[1] = yrange[0]; } - if(zrange.size()==1) { _zrange[0] = _zrange[1] = zrange[0]; } - if(xrange.size()==2) { _xrange[0] = xrange[0]; _xrange[1] = xrange[1]; } - if(yrange.size()==2) { _yrange[0] = yrange[0]; _yrange[1] = yrange[1]; } - if(zrange.size()==2) { _zrange[0] = zrange[0]; _zrange[1] = zrange[1]; } - - _xrange[0] = xmin + _xrange[0]; - _xrange[1] = xmax - _xrange[1]; - _yrange[0] = ymin + _yrange[0]; - _yrange[1] = ymax - _yrange[1]; - _zrange[0] = zmin + _zrange[0]; - _zrange[1] = zmax - _zrange[1]; - - // check - assert(_xrange[0] <= _xrange[1]); - assert(_yrange[0] <= _yrange[1]); - assert(_zrange[0] <= _zrange[1]); - - if(_debug>0) { - std::cout<<"Particle generation world boundaries..."< " << _xrange[1] << std::endl - <<"Y " << _yrange[0] << " => " << _yrange[1] << std::endl - <<"Z " << _zrange[0] << " => " << _zrange[1] << std::endl; - } - - - // register - //auto db = new TDatabasePDG; - auto db = TDatabasePDG::Instance(); - for(size_t idx=0; idxGetParticle(param.pdg[i])->Mass(); - - // sanity check - if(kerange[0]<0 || kerange[1]<0) - this->abort("You provided negative energy? Fuck off Mr. Trump."); - - // overall range check - if(param.kerange[0] > param.kerange[1]) this->abort("KE range has no phase space..."); - - if(_debug>0) { - std::cout << "Generating particle (PDG"; - for(auto const& pdg : param.pdg) std::cout << " " << pdg; - std::cout << ")" << std::endl - << (param.use_mom ? " KE range ....... " : " Mom range ....... ") - << param.kerange[0] << " => " << param.kerange[1] << " MeV" << std::endl - << std::endl; - } - - _param_v.push_back(param); - } -} - -void MultiPartVertex::beginRun(art::Run& run) -{ - // grab the geometry object to see what geometry we are using - art::ServiceHandle geo; - - std::unique_ptr runData(new sumdata::RunData(geo->DetectorName())); - - run.put(std::move(runData), art::fullRun()); - - return; -} - -std::vector MultiPartVertex::GenParticles() const { - - std::vector result; - std::vector gen_count_v(_param_v.size(),0); - - int num_part = (int)(fFlatRandom->fire(_multi_min,_multi_max+1-1.e-10)); - - // generate min multiplicity first - std::vector weight_v(_param_v.size(),0); - for(size_t idx=0; idx<_param_v.size(); ++idx) { - weight_v[idx] = _param_v[idx].weight; - for(size_t ctr=0; ctr<_param_v[idx].multi[0]; ++ctr) { - result.push_back(idx); - gen_count_v[idx] += 1; - num_part -= 1; - } - if(gen_count_v[idx] >= _param_v[idx].multi[1]) - weight_v[idx] = 0.; - } - - assert(num_part >= 0); - - while(num_part) { - - double total_weight = 0; - for(auto const& v : weight_v) total_weight += v; - - double rval = 0; - rval = fFlatRandom->fire(0,total_weight); - - size_t idx = 0; - for(idx = 0; idx < weight_v.size(); ++idx) { - rval -= weight_v[idx]; - if(rval <=0.) break; - } - - // register to the output - result.push_back(idx); - - // if generation count exceeds max, set probability weight to be 0 - gen_count_v[idx] += 1; - if(gen_count_v[idx] >= _param_v[idx].multi[1]) - weight_v[idx] = 0.; - - --num_part; - } - return result; -} - - -void MultiPartVertex::GenPosition(double& x, double& y, double& z) { - - x = fFlatRandom->fire(_xrange[0],_xrange[1]); - y = fFlatRandom->fire(_yrange[0],_yrange[1]); - z = fFlatRandom->fire(_zrange[0],_zrange[1]); - - if(_debug>0) { - std::cout << "Generating a rain particle at (" - << x << "," << y << "," << z << ")" << std::endl; - } -} - -void MultiPartVertex::GenBoost(double& bx, double& by, double& bz) { - - double gbmag = 0.; // saves gamma * beta - if ( _gamma_beta_range[1] > _gamma_beta_range[0]) { - gbmag = fFlatRandom->fire(_gamma_beta_range[0],_gamma_beta_range[1]); - } - else { - gbmag = _gamma_beta_range[0]; - } - double ct = fFlatRandom->fire(-1.0,1.0); - double st = TMath::Sqrt(1.0 - ct*ct); - double phi = fFlatRandom->fire(0.0,2.0 * M_PI); - double bmag = gbmag / TMath::Sqrt(1.0 + gbmag*gbmag); - - bx = bmag * st * TMath::Cos(phi); - by = bmag * st * TMath::Sin(phi); - bz = bmag * ct; - - if(_debug>0) { - std::cout << "Generated boost parameter gamma beta = " << gbmag << ", cos(theta) = " << ct - << ", phi = " << phi << std::endl; - std::cout << "Generating a boost (" - << bx << "," << by << "," << bz << ")" << std::endl; - } -} - -/* -void MultiPartVertex::GenPosition(double& x, double& y, double& z) { - - auto const& tpc_id = _tpc_v.at((size_t)(fFlatRandom->fire(0,_tpc_v.size()))); - // Implementation of required member function here. - auto geop = lar::providerFrom(); - size_t cid = tpc_id[0]; - size_t tid = tpc_id[1]; - auto const& cryostat = geop->Cryostat(cid); - if(!cryostat.HasTPC(tid)) { - std::cerr<< "\033[93mTPC " << tid << " not found... in cryostat " << cid << "\033[00m" << std::endl; - throw std::exception(); - } - - auto const& tpc = cryostat.TPC(tid); - auto const& tpcabox = tpc.ActiveBoundingBox(); - double xmin = tpcabox.MinX() + _xrange[0]; - double xmax = tpcabox.MaxX() - _xrange[1]; - double ymin = tpcabox.MinY() + _yrange[0]; - double ymax = tpcabox.MaxY() - _yrange[1]; - double zmin = tpcabox.MinZ() + _zrange[0]; -double zmax = tpcabox.MaxZ() - _zrange[1]; -x = fFlatRandom->fire(xmin,xmax); -y = fFlatRandom->fire(ymin,ymax); -z = fFlatRandom->fire(zmin,zmax); - -} -*/ - -std::array MultiPartVertex::extractDirection() const { - double ct = fFlatRandom->fire(-1.0,1.0); - double st = TMath::Sqrt(1.0 - ct*ct); - double phi = fFlatRandom->fire(0.0,2.0*M_PI); - double px = st * TMath::Cos(phi); - double py = st * TMath::Sin(phi); - double pz = ct; - std::array result = { px, py, pz }; - return result; -} -/* -std::array MultiPartVertex::extractDirection() const { - double px, py, pz; - px = fNormalRandom->fire(0, 1); - py = fNormalRandom->fire(0, 1); - pz = fNormalRandom->fire(0, 1); - double p = std::hypot(px, py, pz); - px = px / p; - py = py / p; - pz = pz / p; - std::array result = { px, py, pz }; - return result; -} -*/ - -void MultiPartVertex::GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz) { - - double tot_energy = 0; - if(param.use_mom) - tot_energy = sqrt(cet::square(fFlatRandom->fire(param.kerange[0],param.kerange[1])) + cet::square(mass)); - else - tot_energy = fFlatRandom->fire(param.kerange[0],param.kerange[1]) + mass; - - double mom_mag = sqrt(cet::square(tot_energy) - cet::square(mass)); - - std::array p = extractDirection(); - px = p[0]; py = p[1]; pz = p[2]; - - if(_debug>1) - std::cout << " Direction : (" << px << "," << py << "," << pz << ")" << std::endl - << " Momentum : " << mom_mag << " [MeV/c]" << std::endl - << " Energy : " << tot_energy << " [MeV/c^2]" << std::endl; - px *= mom_mag; - py *= mom_mag; - pz *= mom_mag; - -} - -void MultiPartVertex::GenMomentum(const PartGenParam& param, const double& mass, double& px, double& py, double& pz, bool& same_range) { - - if (!same_range){ - GenMomentum(param, mass, px, py, pz); - } - else{ - double tot_energy = 0; - if(param.use_mom) - tot_energy = sqrt(cet::square(fFlatRandom->fire(param.kerange[0],param.kerange[1])) + cet::square(mass)); - else - tot_energy = fFlatRandom->fire(param.kerange[0],1) + mass; - - double mom_mag = sqrt(cet::square(tot_energy) - cet::square(mass)); - - std::array p = extractDirection(); - px = p[0]; py = p[1]; pz = p[2]; - - if(_debug>1) - std::cout << " Direction : (" << px << "," << py << "," << pz << ")" << std::endl - << " Momentum : " << mom_mag << " [MeV/c]" << std::endl - << " Energy : " << tot_energy << " [MeV/c^2]" << std::endl; - px *= mom_mag; - py *= mom_mag; - pz *= mom_mag; - } -} - - -void MultiPartVertex::GenMomentumSF(const double& sf, const double& m, const double& p, double& p_sf) -{ - p_sf = sqrt( sf*( 2*(sf-1)*cet::square(m)-2*(sf-1)*m*sqrt(cet::square(m)+cet::square(p))+sf*cet::square(p) ) )/p; -} - -void MultiPartVertex::produce(art::Event & e) -{ - if(_debug>0) std::cout << "Processing a new event..." << std::endl; - - std::unique_ptr< std::vector > mctArray(new std::vector); - - double g4_time = fFlatRandom->fire(_t0 - _t0_sigma/2., _t0 + _t0_sigma/2.); - - double x, y, z; - GenPosition(x,y,z); - TLorentzVector pos(x,y,z,g4_time); - - double bx,by,bz; - GenBoost(bx,by,bz); - - simb::MCTruth mct; - - mct.SetOrigin(simb::kBeamNeutrino); - - std::vector part_v; - - std::vector E_vec; - std::vector mass_vec; - std::vector pdg_vec; - std::vector px_vec; - std::vector py_vec; - std::vector pz_vec; - - auto const param_idx_v = GenParticles(); - if(_debug) - std::cout << "Event Vertex @ (" << x << "," << y << "," << z << ") ... " << param_idx_v.size() << " particles..." << std::endl; - - for(size_t idx=0; idxfire(0,param.pdg.size()-1.e-10)); - auto const& pdg = param.pdg[pdg_index]; - auto const& mass = param.mass[pdg_index]; - if(_debug) std::cout << " " << idx << "th instance PDG " << pdg << std::endl; - - if (_revert) { - GenMomentum(param,mass,px,py,pz); - } - else { - GenMomentum(param,mass,px,py,pz,same_range); - } - // moving boost - // - //TLorentzVector mom(px,py,pz,sqrt(cet::square(px)+cet::square(py)+cet::square(pz)+cet::square(mass))); - // - /* - if(_use_boost){ - mom.Boost(bx,by,bz); - px=mom.X(); - py=mom.Y(); - pz=mom.Z(); - }*/ - - E_vec.push_back(sqrt(cet::square(px)+cet::square(py)+cet::square(pz)+cet::square(mass))); - mass_vec.push_back(mass); - pdg_vec.push_back(pdg); - px_vec.push_back(px); - py_vec.push_back(py); - pz_vec.push_back(pz); - - } - - if(_revert){ - for(size_t idx=0; idx1.){ - - for(size_t idx=0; idxfire(0,1); - double total_KE_sf = gen_total_KE/total_KE; - if(_debug) std::cout << "gen_total_KE: "<< gen_total_KE << " , KE scale factor : " << total_KE_sf << std::endl; - - for(size_t idx=0; idxpush_back(mct); - - e.put(std::move(mctArray)); -} - -DEFINE_ART_MODULE(MultiPartVertex) diff --git a/sbncode/EventGenerator/MultiPart/gen_mpvmpr.fcl b/sbncode/EventGenerator/MultiPart/gen_mpvmpr.fcl index 2b046e5e0..6193d8937 100644 --- a/sbncode/EventGenerator/MultiPart/gen_mpvmpr.fcl +++ b/sbncode/EventGenerator/MultiPart/gen_mpvmpr.fcl @@ -63,8 +63,8 @@ physics.producers.generator.ZRange : [30,30] physics.producers.generator.TPCRange : [[0,0],[0,1],[0,2],[0,3]] #[1,0],[1,1],[1,2],[1,3]] physics.producers.generator.MultiMax : 7 #6 physics.producers.generator.GammaBetaRange : [0.0, 3.0] # [0.0,20.0] # Boost, GammaBeta = Gamma * Beta -physics.producers.generator.Revert: 0 #1 -physics.producers.generator.UseBoost : 0# 1 +physics.producers.generator.Revert: false #1 +physics.producers.generator.UseBoost : false# 1 physics.producers.generator.MultiMin : 2 physics.producers.generator.ParticleParameter.PDGCode : [[-11,11,-13,13], [111], [211,-211], [2212], [22]] physics.producers.generator.ParticleParameter.MinMulti : [ 0, 0, 0, 0, 0]