Interface to the CRY cosmic-ray generator. More...

#include "CRYHelper.h"

Public Member Functions
	CRYHelper ()

	CRYHelper (fhicl::ParameterSet const &pset, CLHEP::HepRandomEngine &engine, std::string const &worldVol="vWorld")

	~CRYHelper ()

double	Sample (simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)

Private Member Functions
void	WorldBox (double xlo_cm, double xhi_cm, double ylo_cm, double yhi_cm, double zlo_cm, double zhi_cm) const

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

Private Attributes
CRYSetup *	fSetup
	CRY configuration. More...

CRYGenerator *	fGen
	The CRY generator. More...

double	fSampleTime
	Amount of time to sample (seconds) More...

double	fToffset
	Shift in time of particles (s) More...

double	fEthresh
	Cut on kinetic energy (GeV) More...

std::string	fWorldVolume
	Name of the world volume. More...

std::string	fLatitude
	Latitude of detector need space after value. More...

std::string	fAltitude
	Altitude of detector need space after value. More...

std::string	fSubBoxL
	Length of subbox (m) need space after value. More...

double	fBoxDelta

bool	fSingleEventMode
	flag to turn on producing only a single cosmic ray More...

Detailed Description

Interface to the CRY cosmic-ray generator.

Definition at line 26 of file CRYHelper.h.

Constructor & Destructor Documentation

evgb::CRYHelper::CRYHelper ( )

Definition at line 36 of file CRYHelper.cxx.

37 {

38 }

evgb::CRYHelper::CRYHelper	(	fhicl::ParameterSet const &	pset,
		CLHEP::HepRandomEngine &	engine,
		std::string const &	worldVol = `"vWorld"`
	)

explicit

Definition at line 41 of file CRYHelper.cxx.

References fAltitude, fGen, fLatitude, fSetup, fSubBoxL, fhicl::ParameterSet::get(), and evgb::RNGWrapper< T >::set().

     : fSampleTime     (pset.get< double      >("SampleTime")            )
     , fToffset        (pset.get< double      >("TimeOffset")            )
     , fEthresh        (pset.get< double      >("EnergyThreshold")       )
     , fWorldVolume    (worldVol)
     , fLatitude       (pset.get< std::string >("Latitude")              )
     , fAltitude       (pset.get< std::string >("Altitude")              )
     , fSubBoxL        (pset.get< std::string >("SubBoxLength")          )
     , fBoxDelta       (pset.get< double      >("WorldBoxDelta", 1.e-5)  )
     , fSingleEventMode(pset.get< bool        >("GenSingleEvents", false))
   {    
     // Construct the CRY generator
     std::string config("date 1-1-2014 ");
 
     // all particles are turned on by default.  have to have trailing space if 
     // configured in .fcl file
     config += pset.get< std::string >("GammaSetting",    "returnGammas    1 ");
     config += pset.get< std::string >("ElectronSetting", "returnElectrons 1 ");
     config += pset.get< std::string >("MuonSetting",     "returnMuons     1 ");
     config += pset.get< std::string >("PionSetting",     "returnPions     1 ");
     config += pset.get< std::string >("NeutronSetting",  "returnNeutrons  1 ");
     config += pset.get< std::string >("ProtonSetting",   "returnProtons   1 ");
     config += fLatitude;
     config += fAltitude;
     config += fSubBoxL;
 
     // Find the pointer to the CRY data tables
     std::string crydatadir;
     const char* datapath = getenv("CRYDATAPATH");
     if( datapath != 0) crydatadir = datapath;
     else{
       mf::LogError("CRYHelper") << "no variable CRYDATAPATH set for cry data location, bail";
       exit(0);
     }
       
     // Construct the event generator object
     fSetup = new CRYSetup(config, crydatadir);
 
     RNGWrapper<CLHEP::HepRandomEngine>::set(&engine, &CLHEP::HepRandomEngine::flat);
 
     fSetup->setRandomFunction(RNGWrapper<CLHEP::HepRandomEngine>::rng);
 
     fGen = new CRYGenerator(fSetup);
 
   }  

evgb::CRYHelper::~CRYHelper ( )

Definition at line 90 of file CRYHelper.cxx.

References fGen, and fSetup.

   {
     delete fGen;
     delete fSetup;
   }

Member Function Documentation

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

private

Project along a direction from a particular starting point to the edge of a box

Parameters

xyz	- The starting x,y,z location. Must be inside box.
dxyz	- Direction vector
xlo	- Low edge of box in x
xhi	- Low edge of box in x
ylo	- Low edge of box in y
yhi	- Low edge of box in y
zlo	- Low edge of box in z
zhi	- Low edge of box in z
xyzout	- On output, the position at the box edge

Note: It should be safe to use the same array for input and output.

Definition at line 264 of file CRYHelper.cxx.

References d, and fBoxDelta.

Referenced by Sample().

   {
     // make the world box slightly smaller so that the projection to 
     // the edge avoids possible rounding errors later on with Geant4
     xlo += fBoxDelta;
     xhi -= fBoxDelta;
     ylo += fBoxDelta;
     yhi -= fBoxDelta;
     zlo += fBoxDelta;
     zhi -= fBoxDelta;
 
     // Make sure we're inside the box!
     if(xyz[0] < xlo || xyz[0] > xhi ||
        xyz[1] < ylo || xyz[1] > yhi || 
        xyz[2] < zlo || xyz[2] > zhi)
       throw cet::exception("CRYHelper") << "Projection to edge is outside"
                                         << " bounds of world box:\n "
                                         << "\tx: " << xyz[0] << " ("
                                         << xlo << "," << xhi << ")\n"
                                         << "\ty: " << xyz[1] << " ("
                                         << ylo << "," << yhi << ")\n"
                                         << "\tz: " << xyz[2] << " ("
                                         << zlo << "," << zhi << ")";
 
     // 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;
     }
   }

double evgb::CRYHelper::Sample	(	simb::MCTruth &	mctruth,
		double const &	surfaceY,
		double const &	detectorLength,
		double *	w,
		double	rantime = `0`
	)

Todo:: Check if this time slice passes selection criteria

Definition at line 97 of file CRYHelper.cxx.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), fEthresh, fGen, fSampleTime, fSingleEventMode, fToffset, simb::kCosmicRay, evgb::kCosmicRayGenerator, LOG_DEBUG, ProjectToBoxEdge(), simb::MCTruth::SetOrigin(), WorldBox(), x1, x2, y1, and y2.

Referenced by evgen::CosmicsGen::produce().

   {
     // Generator time at start of sample
     double tstart = fGen->timeSimulated();
     int    idctr = 1;
     bool particlespushed = false;
     while (1) {
       std::vector<CRYParticle*> parts;
       fGen->genEvent(&parts);
       for (unsigned int i=0; i<parts.size(); ++i) {
         
         // Take ownership of the particle from the vector
         std::unique_ptr<CRYParticle> cryp(parts[i]);
         
         // Pull out the PDG code
         int pdg = cryp->PDGid();
         
         // Get the energies of the particles
         double ke = cryp->ke()*1.0E-3; // MeV to GeV conversion
         if (ke<fEthresh) continue;
         
         double m    = 0.; // in GeV
         
         static TDatabasePDG*  pdgt = TDatabasePDG::Instance();
         TParticlePDG* pdgp = pdgt->GetParticle(pdg);
         if (pdgp) m = pdgp->Mass();
         
         double etot = ke + m;
         double ptot = etot*etot-m*m;
         if (ptot>0.0) ptot = sqrt(ptot);
         else          ptot = 0.0;
         
         // Sort out the momentum components. Remember that the NOvA
         // frame has y up and z along the beam. So uvw -> zxy
         double px = ptot * cryp->v();
         double py = ptot * cryp->w();
         double pz = ptot * cryp->u();
         
         // Particle start position. CRY distributes uniformly in x-y
         // plane at fixed z, where z is the vertical direction. This
         // requires some offsets and rotations to put the particles at
         // the surface in the geometry as well as some rotations
         // since the coordinate frame has y up and z along the
         // beam.
         double vx = cryp->y()*100.0;
         double vy = cryp->z()*100.0 + surfaceY;
         double vz = cryp->x()*100.0 + 0.5*detectorLength;
         double t  = cryp->t()-tstart + fToffset; // seconds
         if(fSingleEventMode) t  = fSampleTime*rantime; // seconds
 
         // Project backward to edge of world volume
         double xyz[3]  = { vx,  vy,  vz};
         double xyzo[3];
         double dxyz[3] = {-px, -py, -pz};
         double x1 = 0.;
         double x2 = 0.;
         double y1 = 0.;
         double y2 = 0.;
         double z1 = 0.;
         double z2 = 0.;
         this->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
         
         LOG_DEBUG("CRYHelper") << xyz[0] << " " << xyz[1] << " " << xyz[2] << " " 
                                << x1 << " " << x2 << " " 
                                << y1 << " " << y2 << " " 
                                << z1 << " " << z2;
         
         this->ProjectToBoxEdge(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo);
         
         vx = xyzo[0];
         vy = xyzo[1];
         vz = xyzo[2];
         
         // Boiler plate...
         int istatus    =  1;
         int imother1   = kCosmicRayGenerator;
         
         // Push the particle onto the stack
         std::string primary("primary");
         
         particlespushed=true;
         simb::MCParticle p(idctr,
                            pdg,
                            primary,
                            imother1,
                            m,
                            istatus);
         TLorentzVector pos(vx,vy,vz,t*1e9);// time needs to be in ns to match GENIE, etc
         TLorentzVector mom(px,py,pz,etot);
         p.AddTrajectoryPoint(pos,mom);
         
         mctruth.Add(p);
         ++idctr;
       } // Loop on particles in event
 
       // Check if we're done with this time sample
       // note that now requiring npart==1 in singlevent mode.
       
       if (fGen->timeSimulated()-tstart > fSampleTime || 
           (fSingleEventMode && particlespushed ) 
           ) break;    
     } // Loop on events simulated
     
     mctruth.SetOrigin(simb::kCosmicRay);
 
     if (w) *w = 1.0;
     return fGen->timeSimulated()-tstart;
 
   }

void evgb::CRYHelper::WorldBox	(	double *	xlo,
		double *	xhi,
		double *	ylo,
		double *	yhi,
		double *	zlo,
		double *	zhi
	)		const

private

Return the ranges of x,y and z for the "world volume" that the entire geometry lives in. If any pointers are 0, then those coordinates are ignored.

Parameters

xlo	: On return, lower bound on x positions
xhi	: On return, upper bound on x positions
ylo	: On return, lower bound on y positions
yhi	: On return, upper bound on y positions
zlo	: On return, lower bound on z positions
zhi	: On return, upper bound on z positions

Definition at line 225 of file CRYHelper.cxx.

References fWorldVolume, s, x1, x2, y1, and y2.

Referenced by Sample().

   {
     const TGeoShape* s = gGeoManager->GetVolume(fWorldVolume.c_str())->GetShape();
     if(!s)
       throw cet::exception("CRYHelper") << "No TGeoShape found for world volume";
     
     if (xlo || xhi) {
       double x1, x2;
       s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
     }
     if (ylo || yhi) {
       double y1, y2;
       s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
     }
     if (zlo || zhi) {
       double z1, z2;
       s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
     }
   }// end of WorldBox

Member Data Documentation

std::string evgb::CRYHelper::fAltitude

private

Altitude of detector need space after value.

Definition at line 63 of file CRYHelper.h.

Referenced by CRYHelper().

double evgb::CRYHelper::fBoxDelta

private

Adjustment to the size of the world box in each dimension to avoid G4 rounding errors

Definition at line 65 of file CRYHelper.h.

Referenced by ProjectToBoxEdge().

double evgb::CRYHelper::fEthresh

private

Cut on kinetic energy (GeV)

Definition at line 60 of file CRYHelper.h.

Referenced by Sample().

CRYGenerator* evgb::CRYHelper::fGen

private

The CRY generator.

Definition at line 57 of file CRYHelper.h.

Referenced by CRYHelper(), Sample(), and ~CRYHelper().

std::string evgb::CRYHelper::fLatitude

private

Latitude of detector need space after value.

Definition at line 62 of file CRYHelper.h.

Referenced by CRYHelper().

double evgb::CRYHelper::fSampleTime

private

Amount of time to sample (seconds)

Definition at line 58 of file CRYHelper.h.

Referenced by Sample().

CRYSetup* evgb::CRYHelper::fSetup

private

CRY configuration.

Definition at line 56 of file CRYHelper.h.

Referenced by CRYHelper(), and ~CRYHelper().

bool evgb::CRYHelper::fSingleEventMode

private

flag to turn on producing only a single cosmic ray

Definition at line 67 of file CRYHelper.h.

Referenced by Sample().

std::string evgb::CRYHelper::fSubBoxL

private

Length of subbox (m) need space after value.

Definition at line 64 of file CRYHelper.h.

Referenced by CRYHelper().

double evgb::CRYHelper::fToffset

private

Shift in time of particles (s)

Definition at line 59 of file CRYHelper.h.

Referenced by Sample().

std::string evgb::CRYHelper::fWorldVolume

private

Name of the world volume.

Definition at line 61 of file CRYHelper.h.

Referenced by WorldBox().

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

nutools/EventGeneratorBase/CRY/CRYHelper.h
nutools/EventGeneratorBase/CRY/CRYHelper.cxx

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation