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

Inheritance diagram for evgen::EventGeneratorTest:

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

void	analyze (art::Event const &evt)

Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDAnalyzer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept

	Analyzer (fhicl::ParameterSet const &pset)

template<typename Config >
	Analyzer (Table< Config > const &config)

void	doBeginJob ()

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::Observer
	~Observer () noexcept

	Observer (Observer const &)=delete

	Observer (Observer &&)=delete

Observer &	operator= (Observer const &)=delete

Observer &	operator= (Observer &&)=delete

void	registerProducts (ProductDescriptions &, ModuleDescription const &)

void	fillDescriptions (ModuleDescription const &)

std::string const &	processName () const

bool	wantAllEvents () const

bool	wantEvent (Event const &e)

fhicl::ParameterSetID	selectorConfig () const

Handle< TriggerResults >	getTriggerResults (Event const &e) const

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
fhicl::ParameterSet	GENIEParameterSet (std::string fluxType, bool usePOTPerSpill)

void	GENIETest (fhicl::ParameterSet const &pset)

void	GENIEHistogramFluxTest ()

void	GENIESimpleFluxTest ()

void	GENIEMonoFluxTest ()

void	GENIEAtmoFluxTest ()

void	GENIENtupleFluxTest ()

fhicl::ParameterSet	CRYParameterSet ()

void	CRYTest ()

bool	IntersectsDetector (simb::MCParticle const &part)

void	ProjectToSurface (TLorentzVector pos, TLorentzVector mom, int axis, double surfaceLoc, double *xyz)

Private Attributes
double	fTotalGENIEPOT

double	fTotalGENIEInteractions

double	fTotalCRYSpills
	number of spills to use when testing CRY More...

std::string	fTopVolume
	Top Volume used by GENIE. More...

std::string	fGeometryFile
	location of Geometry GDML file to test More...

double	fCryDetLength
	length of detector to test CRY, units of cm More...

double	fCryDetWidth
	width of detector to test CRY, units of cm More...

double	fCryDetHeight
	height of detector to test CRY, units of cm More...

CLHEP::HepRandomEngine &	fEngine
	art-owned engine used in generation random numbers More...

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >

using	ModuleType = EDAnalyzer

Protected Member Functions inherited from art::Observer
	Observer (fhicl::ParameterSet const &config)

	Observer (std::vector< std::string > const &paths, fhicl::ParameterSet const &config)

detail::ProcessAndEventSelectors &	processAndEventSelectors ()

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

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

Definition at line 46 of file EventGeneratorTest_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 93 of file EventGeneratorTest_module.cc.

     : EDAnalyzer             (pset)
     , fTotalGENIEPOT         ( pset.get< double      >("TotalGENIEPOT",          5e18))
     , fTotalGENIEInteractions( pset.get< double      >("TotalGENIEInteractions", 100) )
     , fTotalCRYSpills        ( pset.get< double      >("TotalCRYSpills",         1000))
     , fTopVolume             ( pset.get< std::string >("TopVolume"                   ))
     , fGeometryFile          ( pset.get< std::string >("GeometryFile"                ))
     , fCryDetLength          (1000.)
     , fCryDetWidth           (500.)
     , fCryDetHeight          (500.)
     , fEngine{createEngine(pset.get< int >("Seed", evgb::GetRandomNumberSeed()))}
   {}

Member Function Documentation

void evgen::EventGeneratorTest::analyze ( art::Event const & evt )

virtual

Implements art::EDAnalyzer.

Definition at line 107 of file EventGeneratorTest_module.cc.

   {
     mf::LogWarning("EventGeneratorTest") << "testing GENIE...";
     mf::LogWarning("EventGeneratorTest") << "\t histogram flux...";
     this->GENIEHistogramFluxTest();
     mf::LogWarning("EventGeneratorTest") << "\t \t done."
                                          << "\t simple flux...";
     this->GENIESimpleFluxTest();
     mf::LogWarning("EventGeneratorTest") << "\t \t done."
                                          << "\t atmo flux...";
     this->GENIEAtmoFluxTest();
     mf::LogWarning("EventGeneratorTest") << "\t \t done."
                                          << "\t mono flux...";
     this->GENIEMonoFluxTest();
     mf::LogWarning("EventGeneratorTest") << "\t \t done.\n"
                                          << "GENIE tests done";
 
     mf::LogWarning("EventGeneratorTest") << "testing CRY...";
     this->CRYTest();
     mf::LogWarning("EventGeneratorTest") << "\t CRY test done.";
   }

fhicl::ParameterSet evgen::EventGeneratorTest::CRYParameterSet ( )

private

Definition at line 340 of file EventGeneratorTest_module.cc.

   {
     fhicl::ParameterSet pset;
     pset.put("SampleTime",       600e-6            );
     pset.put("TimeOffset",      -30e-6             );
     pset.put("EnergyThreshold",  50e-3             );
     pset.put("Latitude",         "latitude 41.8 "  );
     pset.put("Altitude",         "altitude 0 "     );
     pset.put("SubBoxLength",     "subboxLength 75 ");
 
     mf::LogWarning("EventGeneratorTest") << pset.to_string();
 
     return pset;
   }

void evgen::EventGeneratorTest::CRYTest ( )

private

Definition at line 356 of file EventGeneratorTest_module.cc.

   {
     // make the parameter set
     fhicl::ParameterSet pset = this->CRYParameterSet();
 
     // make the CRYHelper
     evgb::CRYHelper help(pset, fEngine);
 
     int    nspill         = 0;
     double avPartPerSpill = 0.;
     double avPartIntersectPerSpill = 0.;
     double avMuonIntersectPerSpill = 0.;
     double avEIntersectPerSpill    = 0.;
     while(nspill < TMath::Nint(fTotalCRYSpills) ){
 
       simb::MCTruth mct;
       
       help.Sample(mct, 
                   1.,
                   100.,
                   0);
 
       avPartPerSpill += mct.NParticles();
 
       // now check to see if the particles go through the 
       // detector enclosure
       for(int p = 0; p < mct.NParticles(); ++p){
         if(this->IntersectsDetector(mct.GetParticle(p)) ){
            avPartIntersectPerSpill += 1.;
            if(TMath::Abs(mct.GetParticle(p).PdgCode()) == 13)
              avMuonIntersectPerSpill += 1.;
            else if(TMath::Abs(mct.GetParticle(p).PdgCode()) == 11)
              avEIntersectPerSpill += 1.;
         }
       }
     
 
 
       ++nspill;
     }
 
     mf::LogWarning("EventGeneratorTest") << "there are " << avPartPerSpill/(1.*nspill)
                                          << " cosmic rays made per spill \n"
                                          << avPartIntersectPerSpill/(1.*nspill)
                                          << " intersect the detector per spill"
                                          << "\n\t " 
                                          << avMuonIntersectPerSpill/(1.*nspill)
                                          << " muons \n\t"
                                          << avEIntersectPerSpill/(1.*nspill)
                                          << " electrons";
 
 
     return;
   }

void evgen::EventGeneratorTest::GENIEAtmoFluxTest ( )

private

Definition at line 322 of file EventGeneratorTest_module.cc.

   {
     
     // make the parameter set
     fhicl::ParameterSet pset1 = this->GENIEParameterSet("atmo_FLUKA", false);
     
     mf::LogWarning("EventGeneratorTest") << "\t\t 1 event per spill...\n";
     
     this->GENIETest(pset1);
     
     return;
     
   }

void evgen::EventGeneratorTest::GENIEHistogramFluxTest ( )

private

Definition at line 269 of file EventGeneratorTest_module.cc.

   {
 
     mf::LogWarning("EventGeneratorTest") << "\t\t\t 1 event per spill...\n";
 
     // make the parameter set
     fhicl::ParameterSet pset1(this->GENIEParameterSet("histogram", false));
 
     this->GENIETest(pset1);
 
     mf::LogWarning("EventGeneratorTest") <<"\t\t\t events based on POT per spill...\n";
 
     fhicl::ParameterSet pset2(this->GENIEParameterSet("histogram", true));
     this->GENIETest(pset2);
 
     return;
   } 

void evgen::EventGeneratorTest::GENIEMonoFluxTest ( )

private

Definition at line 307 of file EventGeneratorTest_module.cc.

   {
 
     // make the parameter set
     fhicl::ParameterSet pset1 = this->GENIEParameterSet("mono", false);
 
     mf::LogWarning("EventGeneratorTest") << "\t\t 1 event per spill...\n";
 
     this->GENIETest(pset1);
 
     return;
 
   }

void evgen::EventGeneratorTest::GENIENtupleFluxTest ( )

private

fhicl::ParameterSet evgen::EventGeneratorTest::GENIEParameterSet	(	std::string	fluxType,
		bool	usePOTPerSpill
	)

private

Definition at line 130 of file EventGeneratorTest_module.cc.

   {
     // make a parameter set first so that we can pass it to the GENIEHelper
     // object we are going to make
     std::vector<double> beamCenter; 
     beamCenter.push_back(0.0); beamCenter.push_back(0.); beamCenter.push_back(0.0);
 
     std::vector<double> beamDir; 
     beamDir.push_back(0.); beamDir.push_back(0.); beamDir.push_back(1.);
 
     std::vector<int> flavors;
     if(fluxType.compare("atmo_FLUKA") == 0){
       flavors.push_back(14);
     }
     else{
       flavors.push_back(12); flavors.push_back(14); flavors.push_back(-12); flavors.push_back(-14);
     }
     
 
     std::vector<std::string> env;
     env.push_back("GPRODMODE"); env.push_back("YES");
     env.push_back("GEVGL");     env.push_back("Default");
 
     double potPerSpill = 5.e13;
     double eventsPerSpill = 0;
     if(!usePOTPerSpill) eventsPerSpill = 1;
 
     std::vector<std::string> fluxFiles;
     fluxFiles.push_back("samples_for_geniehelper/L010z185i_lowthr_ipndshed.root");
     if(fluxType.compare("simple_flux") == 0){
       fluxFiles.clear();
       fluxFiles.push_back("samples_for_geniehelper/gsimple_NOvA-NDOS_le010z185i_20100521_RHC_lowth_s_00001.root");
     }
     else if(fluxType.compare("atmo_FLUKA") == 0){
       fluxFiles.clear();
       // at FNAL this is installed relative to in /nusoft/data/flux
       fluxFiles.push_back("atmospheric/battistoni/sdave_numu07.dat");
     }
 
     else if(fluxType.compare("ntuple") == 0){
       throw cet::exception("EventGeneratorTest") <<"No ntuple flux file "
                                                  << "exists, bail ungracefully";
     }
 
     fhicl::ParameterSet pset;
     pset.put("FluxType",         fluxType);
     pset.put("FluxFiles",        fluxFiles);
     pset.put("BeamName",         "numi");
     pset.put("TopVolume",        fTopVolume);
     pset.put("EventsPerSpill",   eventsPerSpill);
     pset.put("POTPerSpill",      potPerSpill);
     pset.put("BeamCenter",       beamCenter);
     pset.put("BeamDirection",    beamDir);
     pset.put("GenFlavors",       flavors);    
     pset.put("Environment",      env);
     pset.put("DetectorLocation", "NOvA-ND");
 
     mf::LogWarning("EventGeneratorTest") << pset.to_string();
 
     return pset;
   }

void evgen::EventGeneratorTest::GENIESimpleFluxTest ( )

private

Definition at line 288 of file EventGeneratorTest_module.cc.

   {
 
     // make the parameter set
     mf::LogWarning("EventGeneratorTest") << "testing GENIEHelper in simple_flux mode with \n"
               << "\t 1 event per spill...\n";
 
     fhicl::ParameterSet pset1 = this->GENIEParameterSet("simple_flux", false);
     this->GENIETest(pset1);
 
     mf::LogWarning("EventGeneratorTest") <<"\t events based on POT per spill...\n";
 
     fhicl::ParameterSet pset2 = this->GENIEParameterSet("simple_flux", true);
     this->GENIETest(pset2);
 
     return;
   } 

void evgen::EventGeneratorTest::GENIETest ( fhicl::ParameterSet const & pset )

private

Definition at line 194 of file EventGeneratorTest_module.cc.

   {
     // use cet::search_path to get the Geometry file path
     cet::search_path sp("FW_SEARCH_PATH");
     std::string geometryFile = fGeometryFile;
     if( !sp.find_file(geometryFile, fGeometryFile) )
       throw cet::exception("EventGeneratorTest") << "cannot find geometry file:\n " 
                                                  << geometryFile
                                                  << "\n to test GENIE";
 
     TGeoManager::Import(geometryFile.c_str());
 
     // make the GENIEHelper object
     evgb::GENIEHelper help(pset,
                            gGeoManager,
                            geometryFile,
                            gGeoManager->FindVolumeFast(pset.get< std::string>("TopVolume").c_str())->Weight());
     help.Initialize();
 
     int interactionCount = 0;
 
     int nspill = 0;
     int spillLimit = 0;
 
     // decide if we are in POT/Spill or Events/Spill mode
     double eps = pset.get<double>("EventsPerSpill");
     if(eps > 0.) spillLimit = TMath::Nint(fTotalGENIEInteractions/eps);
     else         spillLimit = 1000;
 
     while(nspill < spillLimit){
       ++nspill;
       while( !help.Stop() ){
 
         simb::MCTruth truth;
         simb::MCFlux  flux;
         simb::GTruth  gTruth;
 
         if( help.Sample(truth, flux, gTruth) )
           ++interactionCount;
 
       } // end creation loop for this spill
 
     } // end loop over spills
 
     // count the POT used and the number of events made
     mf::LogWarning("EventGeneratorTest") << "made " << interactionCount << " interactions with " 
                                       << help.TotalExposure() << " POTs";
 
     // compare to a simple expectation
     double totalExp = 0.;
     if(help.FluxType().compare("histogram") == 0 && pset.get<double>("EventsPerSpill") == 0){
       std::vector<TH1D*> fluxhist = help.FluxHistograms();
 
       if(fluxhist.size() < 1){
         throw cet::exception("EventGeneratorTest") << "using histogram fluxes but no histograms provided!";
       }
       
       // see comments in GENIEHelper::Initialize() for how this calculation was done.
       totalExp = 1.e-38*1.e-20*help.TotalHistFlux();
       totalExp *= help.TotalExposure()*help.TotalMass()/(1.67262158e-27);
 
       mf::LogWarning("EventGeneratorTest") << "expected " << totalExp << " interactions";
       if(std::abs(interactionCount - totalExp) > 3.*std::sqrt(totalExp) ){
         throw cet::exception("EventGeneratorTest") << "generated count is more than "
                                                    << "3 sigma off expectation";
       }
 
     }// end if histogram fluxes
 
 
     return;
 
   }

bool evgen::EventGeneratorTest::IntersectsDetector ( simb::MCParticle const & part )

private

Definition at line 412 of file EventGeneratorTest_module.cc.

   {
 
     // the particle's initial position and momentum
     TLorentzVector pos = part.Position();
     TLorentzVector mom = part.Momentum();
 
     if(TMath::Abs(mom.P()) == 0){
       mf::LogWarning("EventGeneratorTest") << "particle has no momentum!!! bail";
       return false;
     }
 
     double xyz[3] = {0.};
 
     // Checking intersection with 6 planes
 
     // 1. Check intersection with the y = +halfheight plane
     this->ProjectToSurface(pos, mom, 1, 0.5*fCryDetHeight, xyz);
 
     if( TMath::Abs(xyz[0]) <= 0.5*fCryDetWidth 
         && xyz[2] > 0.
         && TMath::Abs(xyz[2]) <= fCryDetLength ) return true;
       
 
     // 2. Check intersection with the +x plane
     this->ProjectToSurface(pos, mom, 0, 0.5*fCryDetWidth, xyz);
 
     if( TMath::Abs(xyz[1]) <= 0.5*fCryDetHeight 
         && xyz[2] > 0.
         && TMath::Abs(xyz[2]) <= fCryDetLength ) return true;
 
     // 3. Check intersection with the -x plane
     this->ProjectToSurface(pos, mom, 0, -0.5*fCryDetWidth, xyz);
 
     if( TMath::Abs(xyz[1]) <= 0.5*fCryDetHeight 
         && xyz[2] > 0.
         && TMath::Abs(xyz[2]) <= fCryDetLength ) return true;
 
     // 4. Check intersection with the z=0 plane
     this->ProjectToSurface(pos, mom, 2, 0., xyz);
 
     if( TMath::Abs(xyz[0]) <= 0.5*fCryDetWidth 
         && TMath::Abs(xyz[1]) <= 0.5*fCryDetHeight ) return true;
 
     // 5. Check intersection with the z=detlength plane
     this->ProjectToSurface(pos, mom, 2, fCryDetLength, xyz);
 
     if( TMath::Abs(xyz[0]) <= 0.5*fCryDetWidth 
         && TMath::Abs(xyz[1]) <= 0.5*fCryDetHeight ) return true;
 
     return false;
   }

void evgen::EventGeneratorTest::ProjectToSurface	(	TLorentzVector	pos,
		TLorentzVector	mom,
		int	axis,
		double	surfaceLoc,
		double *	xyz
	)

private

Definition at line 466 of file EventGeneratorTest_module.cc.

   {
     double momDir = 0.;
     double posDir = 0.;
     if(axis == 0){
       momDir = mom.Px();
       posDir = pos.X();
     }
     else if(axis == 1){
       momDir = mom.Py();
       posDir = pos.X();
     }
     else if(axis == 2){
       momDir = mom.Pz();
       posDir = pos.X();
     }
     
     double ddS        = (momDir/mom.P());
     double length1Dim = (posDir - surfaceLoc);
 
     if(TMath::Abs(ddS) > 0.){
       length1Dim /= ddS;
       xyz[0] = pos.X() + length1Dim*mom.Px()/mom.P();
       xyz[1] = pos.Y() + length1Dim*mom.Py()/mom.P();
       xyz[2] = pos.Z() + length1Dim*mom.Pz()/mom.P();
     }
 
     return;
   }