A module to check the results from the Monte Carlo generator. More...

Inheritance diagram for evgen::CORSIKAGen:

Public Types
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

template<typename UserConfig , typename KeysToIgnore = void>
using	Table = ProducerBase::Table< UserConfig, KeysToIgnore >

Public Member Functions
	CORSIKAGen (fhicl::ParameterSet const &pset)

virtual	~CORSIKAGen ()

void	produce (art::Event &evt)

void	beginJob ()

void	beginRun (art::Run &run)

void	reconfigure (fhicl::ParameterSet const &p)

template<typename PROD , BranchType B = InEvent>
ProductID	getProductID (std::string const &instanceName={}) const

template<typename PROD , BranchType B>
ProductID	getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const

bool	modifiesEvent () const

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename T , art::BranchType BT>
art::ProductToken< T >	consumes (InputTag const &it)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , art::BranchType BT>
art::ViewToken< T >	consumesView (InputTag const &it)

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename T , art::BranchType BT>
art::ProductToken< T >	mayConsume (InputTag const &it)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , art::BranchType BT>
art::ViewToken< T >	mayConsumeView (InputTag const &it)

base_engine_t &	createEngine (seed_t seed)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)

seed_t	get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)

Static Public Member Functions
static cet::exempt_ptr< Consumer >	non_module_context ()

Protected Member Functions
CurrentProcessingContext const *	currentContext () const

void	validateConsumedProduct (BranchType const bt, ProductInfo const &pi)

void	prepareForJob (fhicl::ParameterSet const &pset)

void	showMissingConsumes () const

Private Member Functions
void	openDBs ()

void	populateNShowers ()

void	populateTOffset ()

void	GetSample (simb::MCTruth &)

double	wrapvar (const double var, const double low, const double high)

double	wrapvarBoxNo (const double var, const double low, const double high, int &boxno)

void	ProjectToBoxEdge (const double xyz[], const double dxyz[], const double xlo, const double xhi, const double ylo, const double yhi, const double zlo, const double zhi, double xyzout[])

Private Attributes
int	fShowerInputs =0
	Number of shower inputs to process from. More...

std::vector< double >	fNShowersPerEvent
	Number of showers to put in each event of duration fSampleTime; one per showerinput. More...

std::vector< int >	fMaxShowers

double	fShowerBounds [6] ={0.,0.,0.,0.,0.,0.}
	Boundaries of area over which showers are to be distributed. More...

double	fToffset_corsika =0.
	Timing offset to account for propagation time through atmosphere, populated from db. More...

ifdh_ns::ifdh *	fIFDH =0
	(optional) flux file handling More...

double	fProjectToHeight =0.
	Height to which particles will be projected [cm]. More...

std::vector< std::string >	fShowerInputFiles
	Set of CORSIKA shower data files to use. More...

std::vector< double >	fShowerFluxConstants
	Set of flux constants to be associated with each shower data file. More...

double	fSampleTime =0.
	Duration of sample [s]. More...

double	fToffset =0.
	Time offset of sample, defaults to zero (no offset) [s]. More...

std::vector< double >	fBuffBox
	Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm]. More...

double	fShowerAreaExtension =0.
	Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]. More...

sqlite3 *	fdb [5]
	Pointers to sqlite3 database object, max of 5. More...

double	fRandomXZShift =0.
	Each shower will be shifted by a random amount in xz so that showers won't repeatedly sample the same space [cm]. More...

Detailed Description

A module to check the results from the Monte Carlo generator.

Definition at line 49 of file CORSIKAGen_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

evgen::CORSIKAGen::CORSIKAGen ( fhicl::ParameterSet const & pset )

explicit

Definition at line 99 of file CORSIKAGen_module.cc.

     : fProjectToHeight(p.get< double >("ProjectToHeight",0.)),
       fShowerInputFiles(p.get< std::vector< std::string > >("ShowerInputFiles")),
       fShowerFluxConstants(p.get< std::vector< double > >("ShowerFluxConstants")),
       fSampleTime(p.get< double >("SampleTime",0.)),
       fToffset(p.get< double >("TimeOffset",0.)),
       fBuffBox(p.get< std::vector< double > >("BufferBox",{0.0, 0.0, 0.0, 0.0, 0.0, 0.0})),
       fShowerAreaExtension(p.get< double >("ShowerAreaExtension",0.)),
       fRandomXZShift(p.get< double >("RandomXZShift",0.))
   {
     
     if(fShowerInputFiles.size() != fShowerFluxConstants.size() || fShowerInputFiles.size()==0 || fShowerFluxConstants.size()==0)
       throw cet::exception("CORSIKAGen") << "ShowerInputFiles and ShowerFluxConstants have different or invalid sizes!"<<"\n";
     fShowerInputs=fShowerInputFiles.size();
     
     if(fSampleTime==0.) throw cet::exception("CORSIKAGen") << "SampleTime not set!";
     
     if(fProjectToHeight==0.) mf::LogInfo("CORSIKAGen")<<"Using 0. for fProjectToHeight!"
     ;
     // create a default random engine; obtain the random seed from NuRandomService,
     // unless overridden in configuration with key "Seed"
     art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this, "HepJamesRandom", "gen", p, { "Seed", "SeedGenerator" });
     art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this, "HepJamesRandom", "pois", p, "SeedPoisson");
 
     this->reconfigure(p);
     
     this->openDBs();
     this->populateNShowers();
     this->populateTOffset();
     
     produces< std::vector<simb::MCTruth> >();
     produces< sumdata::RunData, art::InRun >();
     
   }

evgen::CORSIKAGen::~CORSIKAGen ( )

virtual

Definition at line 134 of file CORSIKAGen_module.cc.

References fdb, fIFDH, and fShowerInputs.

                          {
     for(int i=0; i<fShowerInputs; i++){
       sqlite3_close(fdb[i]);
     }
     //cleanup temp files
     fIFDH->cleanup();
   }

Member Function Documentation

void evgen::CORSIKAGen::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 514 of file CORSIKAGen_module.cc.

                            {
     art::ServiceHandle<art::TFileService> tfs;
 
 
   }

void evgen::CORSIKAGen::beginRun ( art::Run & run )

virtual

Reimplemented from art::EDProducer.

Definition at line 520 of file CORSIKAGen_module.cc.

References geo::GeometryCore::DetectorName(), and art::Run::put().

                                       {
 
     // grab the geometry object to see what geometry we are using
     art::ServiceHandle<geo::Geometry> geo;
 
     std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
 
     run.put(std::move(runcol));
 
     return;
   }

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::consumes ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::consumes ( InputTag const & it )

inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::consumesMany ( )

inherited

Definition at line 162 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::consumesView ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::consumesView ( InputTag const & it )

inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t seed )

inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

 {
   return rng()->createEngine(placeholder_schedule_id(), seed);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make
	)

inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make,
		label_t const &	engine_label
	)

inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
 }

CurrentProcessingContext const * art::EDProducer::currentContext ( ) const

protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

   {
     return current_context_.get();
   }

EngineCreator::seed_t EngineCreator::get_seed_value	(	fhicl::ParameterSet const &	pset,
		char const	key[] = `"seed"`,
		seed_t const	implicit_seed = `-1`
	)

inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

 {
   auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
   return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
 }

template<typename PROD , BranchType B>

ProductID art::EDProducer::getProductID ( std::string const & instanceName = {} ) const

inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

   {
     return ProducerBase::getProductID<PROD, B>(moduleDescription_,
                                                instanceName);
   }

template<typename PROD , BranchType B>

ProductID art::ProducerBase::getProductID	(	ModuleDescription const &	moduleDescription,
		std::string const &	instanceName
	)		const

inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

   {
     auto const& pd =
       get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
     return pd.productID();
   }

void evgen::CORSIKAGen::GetSample ( simb::MCTruth & mctruth )

private

Definition at line 372 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), fdb, fMaxShowers, fNShowersPerEvent, fProjectToHeight, fRandomXZShift, fSampleTime, fShowerBounds, fShowerInputs, fToffset, fToffset_corsika, art::RandomNumberGenerator::getEngine(), LOG_DEBUG, ProjectToBoxEdge(), art::EngineCreator::rng(), geo::GeometryCore::WorldBox(), wrapvar(), wrapvarBoxNo(), x, x1, x2, y1, y2, and z.

Referenced by produce().

                                                 {
     //for each input, randomly pull fNShowersPerEvent[i] showers from the Particles table
     //and randomly place them in time (between -fSampleTime/2 and fSampleTime/2)
     //wrap their positions based on the size of the area under consideration
     //based on http://nusoft.fnal.gov/larsoft/doxsvn/html/CRYHelper_8cxx_source.html (Sample)
     
     //query from sqlite db with select * from particles where shower in (select id from showers ORDER BY substr(id*0.51123124141,length(id)+2) limit 100000) ORDER BY substr(shower*0.51123124141,length(shower)+2);
     //where 0.51123124141 is a random seed to allow randomly selecting rows and should be randomly generated for each query
     //the inner order by is to select randomly from the possible shower id's
     //the outer order by is to make sure the shower numbers are ordered randomly (without this, the showers always come out ordered by shower number
     //and 100000 is the number of showers to be selected at random and needs to be less than the number of showers in the showers table 
     
     //TDatabasePDG is for looking up particle masses
     static TDatabasePDG* pdgt = TDatabasePDG::Instance();
     
     //db variables
     sqlite3_stmt *statement;
     const TString kStatement("select shower,pdg,px,py,pz,x,z,t,e from particles where shower in (select id from showers ORDER BY substr(id*%f,length(id)+2) limit %d) ORDER BY substr(shower*%f,length(shower)+2)");
 
     //get rng engine
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine("gen");
     CLHEP::RandFlat flat(engine);
 
     CLHEP::HepRandomEngine &engine_pois = rng->getEngine("pois");
     CLHEP::RandPoissonQ randpois(engine_pois);
 
     // get geometry and figure where to project particles to, based on CRYHelper
     art::ServiceHandle<geo::Geometry> geom;
     double x1, x2;
     double y1, y2;
     double z1, z2;
     geom->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);  
 
     // make the world box slightly smaller so that the projection to 
     // the edge avoids possible rounding errors later on with Geant4
     double fBoxDelta=1.e-5;
     x1 += fBoxDelta;
     x2 -= fBoxDelta;
     y1 += fBoxDelta;
     y2 = fProjectToHeight;
     z1 += fBoxDelta;
     z2 -= fBoxDelta;
                
     //populate mctruth
     int ntotalCtr=0; //count number of particles added to mctruth
     int lastShower=0; //keep track of last shower id so that t can be randomized on every new shower
     int nShowerCntr=0; //keep track of how many showers are left to be added to mctruth
     int nShowerQry=0; //number of showers to query from db
     int shower,pdg;
     double px,py,pz,x,z,tParticleTime,etot,showerTime=0.,showerTimex=0.,showerTimez=0.,showerXOffset=0.,showerZOffset=0.,t;
     for(int i=0; i<fShowerInputs; i++){
       nShowerCntr=randpois.fire(fNShowersPerEvent[i]);
       mf::LogInfo("CORSIKAGEN") << " Shower input " << i << " with mean " << fNShowersPerEvent[i] << " generating " << nShowerCntr;
 
       while(nShowerCntr>0){
         //how many showers should we query?
         if(nShowerCntr>fMaxShowers[i]){
           nShowerQry=fMaxShowers[i]; //take the group size
         }else{
           nShowerQry=nShowerCntr; //take the rest that are needed
         }
         //build and do query to get nshowers
         double thisrnd=flat(); //need a new random number for each query
         TString kthisStatement=TString::Format(kStatement.Data(),thisrnd,nShowerQry,thisrnd);
         LOG_DEBUG("CORSIKAGen")<<"Executing: "<<kthisStatement;
         if ( sqlite3_prepare(fdb[i], kthisStatement.Data(), -1, &statement, 0 ) == SQLITE_OK ){
           int res=0;
           //loop over database rows, pushing particles into mctruth object
           while(1){
             res = sqlite3_step(statement);
             if ( res == SQLITE_ROW ){
               shower=sqlite3_column_int(statement,0);
               if(shower!=lastShower){
                 //each new shower gets its own random time and position offsets
                 showerTime=1e9*(flat()*fSampleTime); //converting from s to ns
                 showerTimex=1e9*(flat()*fSampleTime); //converting from s to ns
                 showerTimez=1e9*(flat()*fSampleTime); //converting from s to ns
                 //and a random offset in both z and x controlled by the fRandomXZShift parameter
                 showerXOffset=flat()*fRandomXZShift - (fRandomXZShift/2);
                 showerZOffset=flat()*fRandomXZShift - (fRandomXZShift/2);
               } 
               pdg=sqlite3_column_int(statement,1);
               //get mass for this particle
               double m = 0.; // in GeV
               TParticlePDG* pdgp = pdgt->GetParticle(pdg);
               if (pdgp) m = pdgp->Mass();
               
               //Note: position/momentum in db have north=-x and west=+z, rotate so that +z is north and +x is west
               //get momentum components
               px=sqlite3_column_double(statement,4);//uboone x=Particlez
               py=sqlite3_column_double(statement,3);
               pz=-sqlite3_column_double(statement,2);//uboone z=-Particlex
               etot=sqlite3_column_double(statement,8);
               
               //get/calculate position components
               int boxnoX=0,boxnoZ=0;
               x=wrapvarBoxNo(sqlite3_column_double(statement,6)+showerXOffset,fShowerBounds[0],fShowerBounds[1],boxnoX);
               z=wrapvarBoxNo(-sqlite3_column_double(statement,5)+showerZOffset,fShowerBounds[4],fShowerBounds[5],boxnoZ);
               tParticleTime=sqlite3_column_double(statement,7); //time offset, includes propagation time from top of atmosphere
               //actual particle time is particle surface arrival time
               //+ shower start time
               //+ global offset (fcl parameter, in s)
               //- propagation time through atmosphere
               //+ boxNo{X,Z} time offset to make grid boxes have different shower times
               t=tParticleTime+showerTime+(1e9*fToffset)-fToffset_corsika + showerTimex*boxnoX + showerTimez*boxnoZ;
               //wrap surface arrival so that it's in the desired time window
               t=wrapvar(t,(1e9*fToffset),1e9*(fToffset+fSampleTime));
               
               simb::MCParticle p(ntotalCtr,pdg,"primary",-200,m,1);
               
               //project back to wordvol/fProjectToHeight
               double xyzo[3];
               double x0[3]={x,fShowerBounds[3],z};
               double dx[3]={px,py,pz};
               this->ProjectToBoxEdge(x0, dx, x1, x2, y1, y2, z1, z2, xyzo);
                             
               TLorentzVector pos(xyzo[0],xyzo[1],xyzo[2],t);// time needs to be in ns to match GENIE, etc
               TLorentzVector mom(px,py,pz,etot);
               p.AddTrajectoryPoint(pos,mom);
               mctruth.Add(p);
               ntotalCtr++;
               lastShower=shower;
             }else if ( res == SQLITE_DONE ){
               break;
             }else{
               throw cet::exception("CORSIKAGen") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
             }
           }
         }else{
           throw cet::exception("CORSIKAGen") << "Error preparing statement: (" <<kthisStatement<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
         }
         nShowerCntr=nShowerCntr-nShowerQry;
       }
     }
   }

template<typename T , BranchType = InEvent>

ProductToken<T> art::Consumer::mayConsume ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::mayConsume ( InputTag const & it )

inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::mayConsumeMany ( )

inherited

Definition at line 205 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const & )

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const & it )

inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

bool art::ProducerBase::modifiesEvent ( ) const

inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

     {
       return true;
     }

cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )

staticinherited

Definition at line 76 of file Consumer.cc.

 {
   static Consumer invalid{InvalidTag{}};
   return &invalid;
 }

void evgen::CORSIKAGen::openDBs ( )

private

Definition at line 181 of file CORSIKAGen_module.cc.

References fdb, fIFDH, fShowerInputFiles, fShowerInputs, art::RandomNumberGenerator::getEngine(), LOG_DEBUG, and art::EngineCreator::rng().

                           {
     //choose files based on fShowerInputFiles, copy them with ifdh, open them
     // for c2: statement is unused
     //sqlite3_stmt *statement;
     //get rng engine
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine("gen");
     CLHEP::RandFlat flat(engine);
     
     //setup ifdh object
     if ( ! fIFDH ) fIFDH = new ifdh_ns::ifdh;
     const char* ifdh_debug_env = std::getenv("IFDH_DEBUG_LEVEL");
     if ( ifdh_debug_env ) {
       mf::LogInfo("CORSIKAGen") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env<<"\n";
       fIFDH->set_debug(ifdh_debug_env);
     }
     
     //get ifdh path for each file in fShowerInputFiles, put into selectedflist
     //if 1 file returned, use that file
     //if >1 file returned, randomly select one file
     //if 0 returned, make exeption for missing files
     std::vector<std::pair<std::string,long>> selectedflist;
     for(int i=0; i<fShowerInputs; i++){
       if(fShowerInputFiles[i].find("*")==std::string::npos){ 
         //if there are no wildcards, don't call findMatchingFiles
         selectedflist.push_back(std::make_pair(fShowerInputFiles[i],0));
         mf::LogInfo("CorsikaGen") << "Selected"<<selectedflist.back().first<<"\n";
           }else{
         //use findMatchingFiles
         std::vector<std::pair<std::string,long>> flist;
                 std::string path(gSystem->DirName(fShowerInputFiles[i].c_str()));
                 std::string pattern(gSystem->BaseName(fShowerInputFiles[i].c_str()));
                 flist = fIFDH->findMatchingFiles(path,pattern);
                 unsigned int selIndex=-1;
                 if(flist.size()==1){ //0th element is the search path:pattern
                         selIndex=0;
                 }else if(flist.size()>1){
                         selIndex= (unsigned int) (flat()*(flist.size()-1)+0.5); //rnd with rounding, dont allow picking the 0th element
                 }else{
                         throw cet::exception("CORSIKAGen") << "No files returned for path:pattern: "<<path<<":"<<pattern<<std::endl;
                 }
                 selectedflist.push_back(flist[selIndex]);
                 mf::LogInfo("CorsikaGen") << "For "<<fShowerInputFiles[i]<<":"<<pattern
         <<"\nFound "<< flist.size() << " candidate files"
         <<"\nChoosing file number "<< selIndex << "\n"
         <<"\nSelected "<<selectedflist.back().first<<"\n";
      }
 
     }
     
     //do the fetching, store local filepaths in locallist
     std::vector<std::string> locallist;
     for(unsigned int i=0; i<selectedflist.size(); i++){
       mf::LogInfo("CorsikaGen")
         << "Fetching: "<<selectedflist[i].first<<" "<<selectedflist[i].second<<"\n";
       std::string fetchedfile(fIFDH->fetchInput(selectedflist[i].first));
       LOG_DEBUG("CorsikaGen") << "Fetched; local path: "<<fetchedfile;
       locallist.push_back(fetchedfile);
     }
     
     //open the files in fShowerInputFilesLocalPaths with sqlite3
     for(unsigned int i=0; i<locallist.size(); i++){
       //prepare and execute statement to attach db file
       int res=sqlite3_open(locallist[i].c_str(),&fdb[i]);
       if (res!= SQLITE_OK)
         throw cet::exception("CORSIKAGen") << "Error opening db: (" <<locallist[i].c_str()<<") ("<<res<<"): " << sqlite3_errmsg(fdb[i]) << "; memory used:<<"<<sqlite3_memory_used()<<"/"<<sqlite3_memory_highwater(0)<<"\n";
       else
         mf::LogInfo("CORSIKAGen")<<"Attached db "<< locallist[i]<<"\n";
     }
   }

void evgen::CORSIKAGen::populateNShowers ( )

private

Definition at line 290 of file CORSIKAGen_module.cc.

References bounds(), geo::GeometryCore::CryostatBoundaries(), fBuffBox, fdb, fMaxShowers, fNShowersPerEvent, fRandomXZShift, fSampleTime, fShowerAreaExtension, fShowerBounds, fShowerFluxConstants, fShowerInputs, fToffset, and geo::GeometryCore::Ncryostats().

                                    {
     //populate vector of the number of showers per event based on:
       //AREA the showers are being distributed over
       //TIME of the event (fSampleTime)
       //flux constants that determine the overall normalizations (fShowerFluxConstants)
       //Energy range over which the sample was generated (ERANGE_*)
       //power spectrum over which the sample was generated (ESLOPE)
   
     
     //compute shower area based on the maximal x,z dimensions of cryostat boundaries + fShowerAreaExtension
     art::ServiceHandle<geo::Geometry> geom;
     for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
       double bounds[6] = {0.};
       geom->CryostatBoundaries(bounds, c); 
       for (unsigned int bnd = 0; bnd<6; bnd++){
         mf::LogVerbatim("CORSIKAGen")<<"Cryo Boundary: "<<bnd<<"="<<bounds[bnd]<<" ( + Buffer="<<fBuffBox[bnd]<<")\n";
         if(fabs(bounds[bnd])>fabs(fShowerBounds[bnd])){
           fShowerBounds[bnd]=bounds[bnd];
         }
       } 
     }
     //add on fShowerAreaExtension without being clever
     fShowerBounds[0] = fShowerBounds[0] - fShowerAreaExtension;
     fShowerBounds[1] = fShowerBounds[1] + fShowerAreaExtension;
     fShowerBounds[4] = fShowerBounds[4] - fShowerAreaExtension;
     fShowerBounds[5] = fShowerBounds[5] + fShowerAreaExtension;
     
     double showersArea=(fShowerBounds[1]/100-fShowerBounds[0]/100)*(fShowerBounds[5]/100-fShowerBounds[4]/100);
     
     mf::LogInfo("CORSIKAGen")
       <<  "Area extended by : "<<fShowerAreaExtension
       <<"\nShowers to be distributed betweeen: x="<<fShowerBounds[0]<<","<<fShowerBounds[1]
                              <<" & z="<<fShowerBounds[4]<<","<<fShowerBounds[5]
       <<"\nShowers to be distributed with random XZ shift: "<<fRandomXZShift<<" cm"
       <<"\nShowers to be distributed over area: "<<showersArea<<" m^2"
       <<"\nShowers to be distributed over time: "<<fSampleTime<<" s"
       <<"\nShowers to be distributed with time offset: "<<fToffset<<" s"
       <<"\nShowers to be distributed at y: "<<fShowerBounds[3]<<" cm"
       ;
     
     //db variables
     sqlite3_stmt *statement;
     const std::string kStatement("select erange_high,erange_low,eslope,nshow from input");
     double upperLimitOfEnergyRange=0.,lowerLimitOfEnergyRange=0.,energySlope=0.,oneMinusGamma=0.,EiToOneMinusGamma=0.,EfToOneMinusGamma=0.;
     
     /*
     // get random number generator engine... currently unused here
     art::ServiceHandle<art::RandomNumberGenerator> rng;
     CLHEP::HepRandomEngine &engine = rng->getEngine("gen");
     CLHEP::RandFlat flat(engine);
     */
     
     for(int i=0; i<fShowerInputs; i++){
         //build and do query to get run info from databases
       //  double thisrnd=flat();//need a new random number for each query
         if ( sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0 ) == SQLITE_OK ){
           int res=0;
           res = sqlite3_step(statement);
           if ( res == SQLITE_ROW ){
             upperLimitOfEnergyRange=sqlite3_column_double(statement,0);
             lowerLimitOfEnergyRange=sqlite3_column_double(statement,1);
             energySlope = sqlite3_column_double(statement,2);
             fMaxShowers.push_back(sqlite3_column_int(statement,3));
             oneMinusGamma = 1 + energySlope;
             EiToOneMinusGamma = pow(lowerLimitOfEnergyRange, oneMinusGamma);
             EfToOneMinusGamma = pow(upperLimitOfEnergyRange, oneMinusGamma);
             mf::LogVerbatim("CORSIKAGen")<<"For showers input "<< i<<" found e_hi="<<upperLimitOfEnergyRange<<", e_lo="<<lowerLimitOfEnergyRange<<", slope="<<energySlope<<", k="<<fShowerFluxConstants[i]<<"\n";
           }else{
             throw cet::exception("CORSIKAGen") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
           }         
         }else{
           throw cet::exception("CORSIKAGen") << "Error preparing statement: (" <<kStatement<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
         }
       
       //this is computed, how?
       double NShowers=( M_PI * showersArea * fShowerFluxConstants[i] * (EfToOneMinusGamma - EiToOneMinusGamma) / oneMinusGamma )*fSampleTime; 
       fNShowersPerEvent.push_back(NShowers);
       mf::LogVerbatim("CORSIKAGen")<<"For showers input "<< i
                                <<" the number of showers per event is "<<(int)NShowers<<"\n";
     }
   }

void evgen::CORSIKAGen::populateTOffset ( )

private

Definition at line 263 of file CORSIKAGen_module.cc.

References fdb, fShowerInputs, and fToffset_corsika.

                                   {
     //populate TOffset_corsika by finding minimum ParticleTime from db file
     
     sqlite3_stmt *statement;
     const std::string kStatement("select min(t) from particles");
     double t=0.;
     
     for(int i=0; i<fShowerInputs; i++){
         //build and do query to get run min(t) from each db
         if ( sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0 ) == SQLITE_OK ){
           int res=0;
           res = sqlite3_step(statement);
           if ( res == SQLITE_ROW ){
             t=sqlite3_column_double(statement,0);
             mf::LogInfo("CORSIKAGen")<<"For showers input "<< i<<" found particles.min(t)="<<t<<"\n";
             if (i==0 || t<fToffset_corsika) fToffset_corsika=t;
           }else{
             throw cet::exception("CORSIKAGen") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
           }         
         }else{
           throw cet::exception("CORSIKAGen") << "Error preparing statement: (" <<kStatement<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
         }
     }
     
     mf::LogInfo("CORSIKAGen")<<"Found corsika timeoffset [ns]: "<< fToffset_corsika<<"\n";
   }

void art::Consumer::prepareForJob ( fhicl::ParameterSet const & pset )

protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

 {
   if (!moduleContext_)
     return;
 
   pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
   for (auto& consumablesPerBranch : consumables_) {
     cet::sort_all(consumablesPerBranch);
   }
 }

void evgen::CORSIKAGen::produce ( art::Event & evt )

virtual

Implements art::EDProducer.

Definition at line 532 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), bounds(), geo::GeometryCore::CryostatBoundaries(), DEFINE_ART_MODULE, fBuffBox, simb::MCTruth::GetParticle(), GetSample(), simb::kCosmicRay, simb::MCParticle::Momentum(), geo::GeometryCore::Ncryostats(), simb::MCTruth::NParticles(), simb::MCParticle::Position(), art::Event::put(), and simb::MCTruth::SetOrigin().

                                        {
     std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
 
     art::ServiceHandle<geo::Geometry> geom;
     
     simb::MCTruth truth;
     truth.SetOrigin(simb::kCosmicRay);
     
     simb::MCTruth pretruth;
     GetSample(pretruth);
     mf::LogInfo("CORSIKAGen")<<"GetSample number of particles returned: "<<pretruth.NParticles()<<"\n";
     // loop over particles in the truth object
     for(int i = 0; i < pretruth.NParticles(); ++i){
       simb::MCParticle particle = pretruth.GetParticle(i);
       const TLorentzVector& v4 = particle.Position();
       const TLorentzVector& p4 = particle.Momentum();
       double x0[3] = {v4.X(),  v4.Y(),  v4.Z() };
       double dx[3] = {p4.Px(), p4.Py(), p4.Pz()};
 
       // now check if the particle goes through any cryostat in the detector
       // if so, add it to the truth object.
       for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
         double bounds[6] = {0.};
         geom->CryostatBoundaries(bounds, c);
     
         //add a buffer box around the cryostat bounds to increase the acceptance and account for scattering
         //By default, the buffer box has zero size
         for (unsigned int cb=0; cb<6; cb++)
            bounds[cb] = bounds[cb]+fBuffBox[cb];
         
         //calculate the intersection point with each cryostat surface
         bool intersects_cryo = false;
         for (int bnd=0; bnd!=6; ++bnd) {
           if (bnd<2) {
             double p2[3] = {bounds[bnd],  x0[1] + (dx[1]/dx[0])*(bounds[bnd] - x0[0]), x0[2] + (dx[2]/dx[0])*(bounds[bnd] - x0[0])};
             if ( p2[1] >= bounds[2] && p2[1] <= bounds[3] && 
                  p2[2] >= bounds[4] && p2[2] <= bounds[5] ) {
               intersects_cryo = true;
               break;
             }
           }
           else if (bnd>=2 && bnd<4) {
             double p2[3] = {x0[0] + (dx[0]/dx[1])*(bounds[bnd] - x0[1]), bounds[bnd], x0[2] + (dx[2]/dx[1])*(bounds[bnd] - x0[1])};
             if ( p2[0] >= bounds[0] && p2[0] <= bounds[1] && 
                  p2[2] >= bounds[4] && p2[2] <= bounds[5] ) {
               intersects_cryo = true;
         break;
             }
           }
           else if (bnd>=4) {
             double p2[3] = {x0[0] + (dx[0]/dx[2])*(bounds[bnd] - x0[2]), x0[1] + (dx[1]/dx[2])*(bounds[bnd] - x0[2]), bounds[bnd]};
             if ( p2[0] >= bounds[0] && p2[0] <= bounds[1] && 
                  p2[1] >= bounds[2] && p2[1] <= bounds[3] ) {
               intersects_cryo = true;
         break;
             }
           }
         }
 
         if (intersects_cryo){
           truth.Add(particle);
           break; //leave loop over cryostats to avoid adding particle multiple times  
         }// end if particle goes into a cryostat
       }// end loop over cryostats in the detector
 
     }// loop on particles
     
     mf::LogInfo("CORSIKAGen")<<"Number of particles from getsample crossing cryostat + bounding box: "<<truth.NParticles()<<"\n";
 
     truthcol->push_back(truth);
     evt.put(std::move(truthcol));
   
     return;
   }// end produce

void evgen::CORSIKAGen::ProjectToBoxEdge	(	const double	xyz[],
		const double	dxyz[],
		const double	xlo,
		const double	xhi,
		const double	ylo,
		const double	yhi,
		const double	zlo,
		const double	zhi,
		double	xyzout[]
	)

private

Definition at line 142 of file CORSIKAGen_module.cc.

References d.

Referenced by GetSample().

                                                           {
                                         
     
     //we want to project backwards, so take mirror of momentum
     const double dxyz[3]={-indxyz[0],-indxyz[1],-indxyz[2]};
       
     // Compute the distances to the x/y/z walls
     double dx = 99.E99;
     double dy = 99.E99;
     double dz = 99.E99;
     if      (dxyz[0] > 0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
     else if (dxyz[0] < 0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
     if      (dxyz[1] > 0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
     else if (dxyz[1] < 0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
     if      (dxyz[2] > 0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
     else if (dxyz[2] < 0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
     
     
     // Choose the shortest distance
     double d = 0.0;
     if      (dx < dy && dx < dz) d = dx;
     else if (dy < dz && dy < dx) d = dy;
     else if (dz < dx && dz < dy) d = dz;
     
     // Make the step
     for (int i = 0; i < 3; ++i) {
       xyzout[i] = xyz[i] + dxyz[i]*d;
     }
     
   }

void evgen::CORSIKAGen::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 509 of file CORSIKAGen_module.cc.

void art::Consumer::showMissingConsumes ( ) const

protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

 {
   if (!moduleContext_)
     return;
 
   // If none of the branches have missing consumes statements, exit early.
   if (std::all_of(cbegin(missingConsumes_),
                   cend(missingConsumes_),
                   [](auto const& perBranch) { return perBranch.empty(); }))
     return;
 
   constexpr cet::HorizontalRule rule{60};
   mf::LogPrint log{"MTdiagnostics"};
   log << '\n'
       << rule('=') << '\n'
       << "The following consumes (or mayConsume) statements are missing from\n"
       << module_context(moduleDescription_) << '\n'
       << rule('-') << '\n';
 
   cet::for_all_with_index(
     missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
       for (auto const& pi : perBranch) {
         log << "  "
             << assemble_consumes_statement(static_cast<BranchType>(i), pi)
             << '\n';
       }
     });
   log << rule('=');
 }

void art::Consumer::validateConsumedProduct	(	BranchType const	bt,
		ProductInfo const &	pi
	)

protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

 {
   // Early exits if consumes tracking has been disabled or if the
   // consumed product is an allowed consumable.
   if (!moduleContext_)
     return;
 
   if (cet::binary_search_all(consumables_[bt], pi))
     return;
 
   if (requireConsumes_) {
     throw Exception(errors::ProductRegistrationFailure,
                     "Consumer: an error occurred during validation of a "
                     "retrieved product\n\n")
       << "The following consumes (or mayConsume) statement is missing from\n"
       << module_context(moduleDescription_) << ":\n\n"
       << "  " << assemble_consumes_statement(bt, pi) << "\n\n";
   }
 
   missingConsumes_[bt].insert(pi);
 }

double evgen::CORSIKAGen::wrapvar	(	const double	var,
		const double	low,
		const double	high
	)

private

Definition at line 252 of file CORSIKAGen_module.cc.

Referenced by GetSample().

                                                                                   {
     //wrap variable so that it's always between low and high
     return (var - (high - low) * floor(var/(high-low))) + low;
   }

double evgen::CORSIKAGen::wrapvarBoxNo	(	const double	var,
		const double	low,
		const double	high,
		int &	boxno
	)

private

Definition at line 257 of file CORSIKAGen_module.cc.

Referenced by GetSample().

                                                                                                    {
     //wrap variable so that it's always between low and high
     boxno=int(floor(var/(high-low)));
     return (var - (high - low) * floor(var/(high-low))) + low;
   }

Member Data Documentation

std::vector<double> evgen::CORSIKAGen::fBuffBox

private

Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm].

Definition at line 90 of file CORSIKAGen_module.cc.

Referenced by populateNShowers(), and produce().

sqlite3* evgen::CORSIKAGen::fdb[5]

private

Pointers to sqlite3 database object, max of 5.

Definition at line 92 of file CORSIKAGen_module.cc.

Referenced by GetSample(), openDBs(), populateNShowers(), populateTOffset(), and ~CORSIKAGen().

ifdh_ns::ifdh* evgen::CORSIKAGen::fIFDH =0

private

(optional) flux file handling

Definition at line 82 of file CORSIKAGen_module.cc.

Referenced by openDBs(), and ~CORSIKAGen().

std::vector<int> evgen::CORSIKAGen::fMaxShowers

private

Definition at line 79 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

std::vector<double> evgen::CORSIKAGen::fNShowersPerEvent

private

Number of showers to put in each event of duration fSampleTime; one per showerinput.

Definition at line 78 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fProjectToHeight =0.

private

Height to which particles will be projected [cm].

Definition at line 85 of file CORSIKAGen_module.cc.

Referenced by GetSample().

double evgen::CORSIKAGen::fRandomXZShift =0.

private

Each shower will be shifted by a random amount in xz so that showers won't repeatedly sample the same space [cm].

Definition at line 93 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fSampleTime =0.

private

Duration of sample [s].

Definition at line 88 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fShowerAreaExtension =0.

private

Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm].

Definition at line 91 of file CORSIKAGen_module.cc.

Referenced by populateNShowers().

double evgen::CORSIKAGen::fShowerBounds[6] ={0.,0.,0.,0.,0.,0.}

private

Boundaries of area over which showers are to be distributed.

Definition at line 80 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

std::vector< double > evgen::CORSIKAGen::fShowerFluxConstants

private

Set of flux constants to be associated with each shower data file.

Definition at line 87 of file CORSIKAGen_module.cc.

Referenced by populateNShowers().

std::vector< std::string > evgen::CORSIKAGen::fShowerInputFiles

private

Set of CORSIKA shower data files to use.

Definition at line 86 of file CORSIKAGen_module.cc.

Referenced by openDBs().

int evgen::CORSIKAGen::fShowerInputs =0

private

Number of shower inputs to process from.

Definition at line 77 of file CORSIKAGen_module.cc.

Referenced by GetSample(), openDBs(), populateNShowers(), populateTOffset(), and ~CORSIKAGen().

double evgen::CORSIKAGen::fToffset =0.

private

Time offset of sample, defaults to zero (no offset) [s].

Definition at line 89 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fToffset_corsika =0.

private

Timing offset to account for propagation time through atmosphere, populated from db.

Definition at line 81 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateTOffset().

The documentation for this class was generated from the following file:

/cvmfs/larsoft.opensciencegrid.org/products/larsim/v07_09_00/source/larsim/EventGenerator/CORSIKA/CORSIKAGen_module.cc

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation