evgendp::Gen311 Class Reference

Inheritance diagram for evgendp::Gen311:

Public Member Functions
	Gen311 (fhicl::ParameterSet const &p)
virtual	~Gen311 ()
	Gen311 (Gen311 const &)=delete
	Gen311 (Gen311 &&)=delete
Gen311 &	operator= (Gen311 const &)=delete
Gen311 &	operator= (Gen311 &&)=delete
void	produce (art::Event &e) override
void	reconfigure (fhicl::ParameterSet const &p)
void	beginJob () override
void	beginRun (art::Run &run) override
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	openDBs ()
void	populateNShowers ()
void	populateTOffset ()
bool	InTPC (const simb::MCParticle particle)
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	GetSample (simb::MCTruth &)

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	fRandomYZShift =0.
	Each shower will be shifted by a random amount in xz so that showers won't repeatedly sample the same space [cm]. More...
std::vector< double >	fActiveVolumeCut
	Active volume cut. More...
std::vector< double >	fLeadingParticlesList
bool	fUseIFDH
	< List of pdg codes for particles that can be potentially used as triggers More...
int	fRun
int	fEvent
double	fTheta
double	fPhi
double	fMomX
double	fMomY
double	fMomZ
double	fMomT
double	fStartX
double	fStartY
double	fStartZ
TTree *	fTree
art::ServiceHandle< geo::Geometry >	geom
CLHEP::HepRandomEngine &	fGenEngine
CLHEP::HepRandomEngine &	fPoisEngine

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 55 of file Gen311_module.cc.

Constructor & Destructor Documentation

evgendp::Gen311::Gen311 ( fhicl::ParameterSet const &  p )

explicit

Definition at line 247 of file Gen311_module.cc.

    : EDProducer{p},
      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.)),
      fRandomYZShift(p.get< double >("RandomYZShift",0.)),
      fActiveVolumeCut(p.get< std::vector< double > >("ActiveVolumeCut")),
      fLeadingParticlesList(p.get< std::vector< double > >("LeadingParticlesList")),
      fUseIFDH(p.get< bool >("UseIFDH")),
      // create a default random engine; obtain the random seed from NuRandomService,
      // unless overridden in configuration with key "Seed"
      fGenEngine(art::ServiceHandle<rndm::NuRandomService>()
                 ->createEngine(*this, "HepJamesRandom", "gen", p, { "Seed", "SeedGenerator" })),
      fPoisEngine(art::ServiceHandle<rndm::NuRandomService>()
                  ->createEngine(*this, "HepJamesRandom", "pois", p, "SeedPoisson"))
   {

    if(fShowerInputFiles.size() != fShowerFluxConstants.size() || fShowerInputFiles.size()==0 || fShowerFluxConstants.size()==0)
      throw cet::exception("Gen311") << "ShowerInputFiles and ShowerFluxConstants have different or invalid sizes!"<<"\n";
    fShowerInputs=fShowerInputFiles.size();

    if(fSampleTime==0.) throw cet::exception("Gen311") << "SampleTime not set!";

    if(fProjectToHeight==0.) mf::LogInfo("Gen311")<<"Using 0. for fProjectToHeight!"
    ;

    this->reconfigure(p);

    this->openDBs();
    this->populateNShowers();
    this->populateTOffset();

    produces< std::vector<simb::MCTruth> >();
    produces< sumdata::RunData, art::InRun >();
   }

evgendp::Gen311::~Gen311 ( )

virtual

Definition at line 287 of file Gen311_module.cc.

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

evgendp::Gen311::Gen311 ( Gen311 const &  )

delete

evgendp::Gen311::Gen311 ( Gen311 &&  )

delete

Member Function Documentation

void evgendp::Gen311::beginJob ( )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 300 of file Gen311_module.cc.

                              {


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

    fTree = tfs->make<TTree>("entries", "entries tree");
    fTree->Branch("Run", &fRun, "Run/I");
    fTree->Branch("Event", &fEvent, "Event/I");
    fTree->Branch("Theta", &fTheta, "Theta/D");
    fTree->Branch("Phi", &fPhi, "Phi/D");
    fTree->Branch("MomX", &fMomX, "MomX/D");
    fTree->Branch("MomY", &fMomY, "MomY/D");
    fTree->Branch("MomZ", &fMomZ, "MomZ/D");
    fTree->Branch("MomT", &fMomT, "MomT/D");
    fTree->Branch("StartX", &fStartX, "StartX/D");
    fTree->Branch("StartY", &fStartY, "StartY/D");
    fTree->Branch("StartZ", &fStartZ, "StartZ/D");

  }

void evgendp::Gen311::beginRun ( art::Run &  run )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 320 of file Gen311_module.cc.

                                         {
   // 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;
  }

void evgendp::Gen311::GetSample ( simb::MCTruth &  mctruth )

private

Definition at line 617 of file Gen311_module.cc.

                                                   {
    //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

    //dummy holder for particles object
    std::map< int, std::vector<simb::MCParticle> >  ParticleMap;
    std::map< int, int >  ShowerTrkIDMap;

    //define the trigger object
    Trigger trg;
    trg.SetCryoBuffer( fBuffBox );
    trg.SetTPCBuffer( fActiveVolumeCut );

    //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)");

    CLHEP::RandFlat flat(fGenEngine);
    CLHEP::RandPoissonQ randpois(fPoisEngine);

    // 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;

    geom->WorldBox(&y1, &y2, &x1, &x2, &z1, &z2);
    y1 += fBoxDelta;
    y2 -= fBoxDelta;
    x1 += fBoxDelta;
    x2 = 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,etot, y=0,z,t; //tParticleTime,,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("Gen311") << " 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);
        MF_LOG_DEBUG("Gen311")<<"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);

              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();

              //get momentum components
              //NR rotation from CORSIKA ref system to LArSoft ref system:
              //LArsoft: px, py, pz;
              //uBoone px', py', pz'; << from the sqlite database instances
              // px = py'
              // py = pz'
              // pz = px'
              px = sqlite3_column_double(statement,3);
              py = sqlite3_column_double(statement,4);
              pz = sqlite3_column_double(statement,2);
              etot=sqlite3_column_double(statement,8);
              y = sqlite3_column_double(statement,6);
              z = sqlite3_column_double(statement,5);
              t = sqlite3_column_double(statement,7); //time offset, includes propagation time from top of atmosphere

              TLorentzVector pos(fShowerBounds[1], y, z,t);// time needs to be in ns to match GENIE, etc
              TLorentzVector mom(px,py,pz,etot);

              simb::MCParticle p(ntotalCtr,pdg,"primary",-200,m,1);
              p.AddTrajectoryPoint(pos,mom);

              //if it is muon, fill up the list of possible triggers
              std::vector<double>::iterator pdgIt;
              pdgIt = find( fLeadingParticlesList.begin(), fLeadingParticlesList.end(), p.PdgCode() );
              if( pdgIt != fLeadingParticlesList.end() ){ trg.AddMuon( p ); }

              ParticleMap[ shower ].push_back(p);
              ShowerTrkIDMap[ ntotalCtr ] = shower; //<< use the unique trackID to trace the source shower back

              ntotalCtr++;
            }else if ( res == SQLITE_DONE ){
              break;
            }else{
              throw cet::exception("Gen311") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
            }
          }
        }else{
          throw cet::exception("Gen311") << "Error preparing statement: (" <<kthisStatement<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
        }
        nShowerCntr=nShowerCntr-nShowerQry;
      }
    } //end loop over showers

    trg.MakeTrigger(fGenEngine); //<<--Select a muon to use as trigger

    double showerTime=0, showerTimey=0, showerTimez=0, showerYOffset=0, showerZOffset=0;
    int boxnoY=0, boxnoZ=0;

    for( auto ParticleMapIt : ParticleMap ){

      if( ParticleMapIt.first == ShowerTrkIDMap[ trg.GetTriggerId() ] ){

        showerTime =1e9*trg.GetTriggerOffsetT(); //converting from s to ns
        showerTimey=1e9*trg.GetTriggerOffsetT();
        showerTimez=1e9*trg.GetTriggerOffsetT();

        //and a random offset in both z and x controlled by the fRandomYZShift parameter
        showerYOffset=1e9*trg.GetTriggerOffsetY();
        showerZOffset=1e9*trg.GetTriggerOffsetZ();

      }else{

        showerTime =1e9*(flat()*fSampleTime); //converting from s to ns
        showerTimey=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 fRandomYZShift parameter
        showerYOffset=flat()*fRandomYZShift - (fRandomYZShift/2);
        showerZOffset=flat()*fRandomYZShift - (fRandomYZShift/2);
      }

      for( auto particleIt : ParticleMapIt.second ){

        simb::MCParticle particle(particleIt.TrackId(),particleIt.PdgCode(),"primary",-200,particleIt.Mass(),1);

        double x0[3]={0.};
        double xyzo[3]={0.};
        double dx[3]={0.};

        if( particleIt.TrackId() == trg.GetTriggerId() ){

          x0[0] = trg.GetTriggerPos().X();
          x0[1] = trg.GetTriggerPos().Y();
          x0[2] = trg.GetTriggerPos().Z();
          dx[0] = trg.GetTriggerMom().X();
          dx[1] = trg.GetTriggerMom().Y();
          dx[2] = trg.GetTriggerMom().Z();

          t = trg.GetTriggerPos().T();

        }else{

          y = wrapvarBoxNo(particleIt.Position().Y()+showerYOffset,fShowerBounds[2],fShowerBounds[3],boxnoY);
          z = wrapvarBoxNo(particleIt.Position().Z()+showerZOffset,fShowerBounds[4],fShowerBounds[5],boxnoZ);

          t=particleIt.Position().T()+showerTime+(1e9*fToffset)-fToffset_corsika + showerTimey*boxnoY + showerTimez*boxnoZ;
          t=wrapvar(t,(1e9*fToffset),1e9*(fToffset+fSampleTime));

          x0[0]= fShowerBounds[1];
          x0[1]= y;
          x0[2]= z;
          dx[0]= particleIt.Momentum().X();
          dx[1]= particleIt.Momentum().Y();
          dx[2]= particleIt.Momentum().Z();

        }

        trg.ProjectToBoxEdge(x0, dx, x1, x2, y1, y2, z1, z2, xyzo);
        TLorentzVector pos(xyzo[0],xyzo[1],xyzo[2], t);
        TLorentzVector mom(dx[0],dx[1],dx[2], particleIt.Momentum().T() );
        particle.AddTrajectoryPoint( pos, mom );

        mctruth.Add( particle );
      } //end loop over particles
    }//loop over showers
  } //end GetSample

bool evgendp::Gen311::InTPC ( const simb::MCParticle  particle )

private

void evgendp::Gen311::openDBs ( )

private

Definition at line 329 of file Gen311_module.cc.

                             {
    //choose files based on fShowerInputFiles, copy them with ifdh, open them
    //sqlite3_stmt *statement;
    CLHEP::RandFlat flat(fGenEngine);

    if(fUseIFDH){
      //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("CORSIKAGendp") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env<<"\n";
        fIFDH->set_debug(ifdh_debug_env);
      }
    }

    //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("CORSIKAGendp") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env<<"\n";
      fIFDH->set_debug(ifdh_debug_env);
    }

    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("Gen311") << "Selected"<<selectedflist.back().first<<"\n";
      }else{
        //read all the file mathching the pattern
        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()));

        if(fUseIFDH){
          //acess the files using ifdh
          flist = fIFDH->findMatchingFiles(path,pattern);
        }else{
          //access the files using drent
          auto wildcardPosition = pattern.find("*");
          pattern = pattern.substr( 0, wildcardPosition );

           DIR *dir;
           struct dirent *ent;
           int index=0;
           if ((dir = opendir( path.c_str() )) != NULL) {
               while ((ent = readdir (dir)) != NULL) {
                 index++;
                 std::pair<std::string,long> name;
                 std::string parsename(ent->d_name);

                 if( parsename.substr(0, wildcardPosition) == pattern ){
                   name.first= path+"/"+parsename;
                   name.second = index;
                   flist.push_back( name );
                }
              }
              closedir (dir);
           }else {
               throw cet::exception("Gen311") << "Can't open directory with pattern: "<<path<<":"<<pattern<<std::endl;
           }
        }

        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("Gen311") << "No files returned for path:pattern: "<<path<<":"<<pattern<<std::endl;
        }
        selectedflist.push_back(flist[selIndex]);
        mf::LogInfo("Gen311") << "For "<<fShowerInputFiles[i]<<":"<<pattern
        <<"\nFound "<< flist.size() << " candidate files"
        <<"\nChoosing file number "<< selIndex << "\n"
        <<"\nSelected "<<selectedflist.back().first<<"\n";
     }
    }

    //open the files in fShowerInputFilesLocalPaths with sqlite3
    std::vector<std::string> locallist;
    for(unsigned int i=0; i<selectedflist.size(); i++){
      mf::LogInfo("Gen311")
        << "Fetching: "<<selectedflist[i].first<<" "<<selectedflist[i].second<<"\n";
      std::string fetchedfile(selectedflist[i].first);
      MF_LOG_DEBUG("Gen311") << "Fetched; local path: "<<fetchedfile;
      locallist.push_back(fetchedfile);
    }

    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("Gen311") << "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("Gen311")<<"Attached db "<< locallist[i]<<"\n";
    }
  }//end openDBs

Gen311& evgendp::Gen311::operator= ( Gen311 const &  )

delete

Gen311& evgendp::Gen311::operator= ( Gen311 &&  )

delete

void evgendp::Gen311::populateNShowers ( )

private

Definition at line 531 of file Gen311_module.cc.

                                      {
    //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("Gen311")<<"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[2] = fShowerBounds[2] - fShowerAreaExtension;
    fShowerBounds[3] = fShowerBounds[3] + fShowerAreaExtension;
    fShowerBounds[4] = fShowerBounds[4] - fShowerAreaExtension;
    fShowerBounds[5] = fShowerBounds[5] + fShowerAreaExtension;

    double showersArea=(fShowerBounds[3]/100-fShowerBounds[2]/100)*(fShowerBounds[5]/100-fShowerBounds[4]/100);

    mf::LogInfo("Gen311")
      <<  "Area extended by : "<<fShowerAreaExtension
      <<"\nShowers to be distributed betweeen: y="<<fShowerBounds[2]<<","<<fShowerBounds[3]
                             <<" & z="<<fShowerBounds[4]<<","<<fShowerBounds[5]
      <<"\nShowers to be distributed with random YZ shift: "<<fRandomYZShift<<" 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 x: "<<fShowerBounds[1]<<" 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.;

    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("Gen311")<<"For showers input "<< i<<" found e_hi="<<upperLimitOfEnergyRange<<", e_lo="<<lowerLimitOfEnergyRange<<", slope="<<energySlope<<", k="<<fShowerFluxConstants[i]<<"\n";
          }else{
            throw cet::exception("Gen311") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
          }
        }else{
          throw cet::exception("Gen311") << "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("Gen311")<<"For showers input "<< i
                               <<" the number of showers per event is "<<(int)NShowers<<"\n";
    }
  }

void evgendp::Gen311::populateTOffset ( )

private

Definition at line 504 of file Gen311_module.cc.

                                     {
    //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("Gen311")<<"For showers input "<< i<<" found particles.min(t)="<<t<<"\n";
            if (i==0 || t<fToffset_corsika) fToffset_corsika=t;
          }else{
            throw cet::exception("Gen311") << "Unexpected sqlite3_step return value: (" <<res<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
          }
        }else{
          throw cet::exception("Gen311") << "Error preparing statement: (" <<kStatement<<"); "<<"ERROR:"<<sqlite3_errmsg(fdb[i])<<"\n";
        }
    }

    mf::LogInfo("Gen311")<<"Found corsika timeoffset [ns]: "<< fToffset_corsika<<"\n";
  }

void evgendp::Gen311::produce ( art::Event &  e )

overridevirtual

Implements art::EDProducer.

Definition at line 1144 of file Gen311_module.cc.

                                         {
    //Fetch the particle generated by GetParticle

    fRun = e.id().run();
    fEvent = e.id().event();

    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("CORSIKAGendp")<<"GetSample number of particles returned: "<<pretruth.NParticles()<<"\n";
    // loop over particles in the truth object

      //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(int i = 0; i < pretruth.NParticles(); ++i){
      simb::MCParticle particle = pretruth.GetParticle(i);

      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("Gen311")<<"Number of particles from getsample crossing cryostat + bounding box: "<<truth.NParticles()<<"\n";

    truthcol->push_back(truth);
    e.put(std::move(truthcol));

    return;
  }

void evgendp::Gen311::reconfigure

(

fhicl::ParameterSet const &

p

)

Definition at line 295 of file Gen311_module.cc.

                                                           {

    return;
  }

double evgendp::Gen311::wrapvar ( const double  var, const double  low, const double  high  )

private

Definition at line 606 of file Gen311_module.cc.

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

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

private

Definition at line 611 of file Gen311_module.cc.

                                                                                                      {
    //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> evgendp::Gen311::fActiveVolumeCut

private

Active volume cut.

Definition at line 91 of file Gen311_module.cc.

std::vector<double> evgendp::Gen311::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 87 of file Gen311_module.cc.

sqlite3* evgendp::Gen311::fdb[5]

private

Pointers to sqlite3 database object, max of 5.

Definition at line 89 of file Gen311_module.cc.

int evgendp::Gen311::fEvent

private

Definition at line 96 of file Gen311_module.cc.

CLHEP::HepRandomEngine& evgendp::Gen311::fGenEngine

private

Definition at line 110 of file Gen311_module.cc.

ifdh_ns::ifdh* evgendp::Gen311::fIFDH =0

private

(optional) flux file handling

Definition at line 79 of file Gen311_module.cc.

std::vector< double > evgendp::Gen311::fLeadingParticlesList

private

Definition at line 92 of file Gen311_module.cc.

std::vector<int> evgendp::Gen311::fMaxShowers

private

Definition at line 76 of file Gen311_module.cc.

double evgendp::Gen311::fMomT

private

Definition at line 103 of file Gen311_module.cc.

double evgendp::Gen311::fMomX

private

Definition at line 100 of file Gen311_module.cc.

double evgendp::Gen311::fMomY

private

Definition at line 101 of file Gen311_module.cc.

double evgendp::Gen311::fMomZ

private

Definition at line 102 of file Gen311_module.cc.

std::vector<double> evgendp::Gen311::fNShowersPerEvent

private

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

Definition at line 75 of file Gen311_module.cc.

double evgendp::Gen311::fPhi

private

Definition at line 99 of file Gen311_module.cc.

CLHEP::HepRandomEngine& evgendp::Gen311::fPoisEngine

private

Definition at line 111 of file Gen311_module.cc.

double evgendp::Gen311::fProjectToHeight =0.

private

Height to which particles will be projected [cm].

Definition at line 82 of file Gen311_module.cc.

double evgendp::Gen311::fRandomYZShift =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 90 of file Gen311_module.cc.

int evgendp::Gen311::fRun

private

Definition at line 95 of file Gen311_module.cc.

double evgendp::Gen311::fSampleTime =0.

private

Duration of sample [s].

Definition at line 85 of file Gen311_module.cc.

double evgendp::Gen311::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 88 of file Gen311_module.cc.

double evgendp::Gen311::fShowerBounds[6] ={0.,0.,0.,0.,0.,0.}

private

Boundaries of area over which showers are to be distributed.

Definition at line 77 of file Gen311_module.cc.

std::vector< double > evgendp::Gen311::fShowerFluxConstants

private

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

Definition at line 84 of file Gen311_module.cc.

std::vector< std::string > evgendp::Gen311::fShowerInputFiles

private

Set of CORSIKA shower data files to use.

Definition at line 83 of file Gen311_module.cc.

int evgendp::Gen311::fShowerInputs =0

private

Number of shower inputs to process from.

Definition at line 74 of file Gen311_module.cc.

double evgendp::Gen311::fStartX

private

Definition at line 104 of file Gen311_module.cc.

double evgendp::Gen311::fStartY

private

Definition at line 105 of file Gen311_module.cc.

double evgendp::Gen311::fStartZ

private

Definition at line 106 of file Gen311_module.cc.

double evgendp::Gen311::fTheta

private

Definition at line 98 of file Gen311_module.cc.

double evgendp::Gen311::fToffset =0.

private

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

Definition at line 86 of file Gen311_module.cc.

double evgendp::Gen311::fToffset_corsika =0.

private

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

Definition at line 78 of file Gen311_module.cc.

TTree* evgendp::Gen311::fTree

private

Definition at line 107 of file Gen311_module.cc.

bool evgendp::Gen311::fUseIFDH

private

< List of pdg codes for particles that can be potentially used as triggers

< use ifdh protocol

Definition at line 93 of file Gen311_module.cc.

art::ServiceHandle<geo::Geometry> evgendp::Gen311::geom

private

Definition at line 109 of file Gen311_module.cc.

---

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

• dunesim/dunesim/EventGenerator/CORSIKAprotodunedp/Gen311_module.cc