LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
evgb::CRYHelper Class Reference

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

44  : fSampleTime (pset.get< double >("SampleTime") )
45  , fToffset (pset.get< double >("TimeOffset") )
46  , fEthresh (pset.get< double >("EnergyThreshold") )
47  , fWorldVolume (worldVol)
48  , fLatitude (pset.get< std::string >("Latitude") )
49  , fAltitude (pset.get< std::string >("Altitude") )
50  , fSubBoxL (pset.get< std::string >("SubBoxLength") )
51  , fBoxDelta (pset.get< double >("WorldBoxDelta", 1.e-5) )
52  , fSingleEventMode(pset.get< bool >("GenSingleEvents", false))
53  {
54  // Construct the CRY generator
55  std::string config("date 1-1-2014 ");
56 
57  // all particles are turned on by default. have to have trailing space if
58  // configured in .fcl file
59  config += pset.get< std::string >("GammaSetting", "returnGammas 1 ");
60  config += pset.get< std::string >("ElectronSetting", "returnElectrons 1 ");
61  config += pset.get< std::string >("MuonSetting", "returnMuons 1 ");
62  config += pset.get< std::string >("PionSetting", "returnPions 1 ");
63  config += pset.get< std::string >("NeutronSetting", "returnNeutrons 1 ");
64  config += pset.get< std::string >("ProtonSetting", "returnProtons 1 ");
65  config += fLatitude;
66  config += fAltitude;
67  config += fSubBoxL;
68 
69  // Find the pointer to the CRY data tables
70  std::string crydatadir;
71  const char* datapath = getenv("CRYDATAPATH");
72  if( datapath != 0) crydatadir = datapath;
73  else{
74  mf::LogError("CRYHelper") << "no variable CRYDATAPATH set for cry data location, bail";
75  exit(0);
76  }
77 
78  // Construct the event generator object
79  fSetup = new CRYSetup(config, crydatadir);
80 
81  RNGWrapper<CLHEP::HepRandomEngine>::set(&engine, &CLHEP::HepRandomEngine::flat);
82 
84 
85  fGen = new CRYGenerator(fSetup);
86 
87  }
double fSampleTime
Amount of time to sample (seconds)
Definition: CRYHelper.h:58
CRYSetup * fSetup
CRY configuration.
Definition: CRYHelper.h:56
std::string fLatitude
Latitude of detector need space after value.
Definition: CRYHelper.h:62
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61
double fBoxDelta
Definition: CRYHelper.h:65
static void set(T *object, double(T::*func)(void))
Definition: CRYHelper.h:84
std::string fAltitude
Altitude of detector need space after value.
Definition: CRYHelper.h:63
bool fSingleEventMode
flag to turn on producing only a single cosmic ray
Definition: CRYHelper.h:67
static double rng(void)
Definition: CRYHelper.h:88
std::string fSubBoxL
Length of subbox (m) need space after value.
Definition: CRYHelper.h:64
double fEthresh
Cut on kinetic energy (GeV)
Definition: CRYHelper.h:60
double fToffset
Shift in time of particles (s)
Definition: CRYHelper.h:59
evgb::CRYHelper::~CRYHelper ( )

Definition at line 90 of file CRYHelper.cxx.

References fGen, and fSetup.

91  {
92  delete fGen;
93  delete fSetup;
94  }
CRYSetup * fSetup
CRY configuration.
Definition: CRYHelper.h:56
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57

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

270  {
271  // make the world box slightly smaller so that the projection to
272  // the edge avoids possible rounding errors later on with Geant4
273  xlo += fBoxDelta;
274  xhi -= fBoxDelta;
275  ylo += fBoxDelta;
276  yhi -= fBoxDelta;
277  zlo += fBoxDelta;
278  zhi -= fBoxDelta;
279 
280  // Make sure we're inside the box!
281  if(xyz[0] < xlo || xyz[0] > xhi ||
282  xyz[1] < ylo || xyz[1] > yhi ||
283  xyz[2] < zlo || xyz[2] > zhi)
284  throw cet::exception("CRYHelper") << "Projection to edge is outside"
285  << " bounds of world box:\n "
286  << "\tx: " << xyz[0] << " ("
287  << xlo << "," << xhi << ")\n"
288  << "\ty: " << xyz[1] << " ("
289  << ylo << "," << yhi << ")\n"
290  << "\tz: " << xyz[2] << " ("
291  << zlo << "," << zhi << ")";
292 
293  // Compute the distances to the x/y/z walls
294  double dx = 99.E99;
295  double dy = 99.E99;
296  double dz = 99.E99;
297  if (dxyz[0] > 0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
298  else if (dxyz[0] < 0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
299  if (dxyz[1] > 0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
300  else if (dxyz[1] < 0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
301  if (dxyz[2] > 0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
302  else if (dxyz[2] < 0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
303 
304  // Choose the shortest distance
305  double d = 0.0;
306  if (dx < dy && dx < dz) d = dx;
307  else if (dy < dz && dy < dx) d = dy;
308  else if (dz < dx && dz < dy) d = dz;
309 
310  // Make the step
311  for (int i = 0; i < 3; ++i) {
312  xyzout[i] = xyz[i] + dxyz[i]*d;
313  }
314  }
double fBoxDelta
Definition: CRYHelper.h:65
Float_t d
Definition: plot.C:235
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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, MF_LOG_DEBUG, ProjectToBoxEdge(), simb::MCTruth::SetOrigin(), WorldBox(), x1, x2, y1, and y2.

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

102  {
103  // Generator time at start of sample
104  double tstart = fGen->timeSimulated();
105  int idctr = 1;
106  bool particlespushed = false;
107  while (1) {
108  std::vector<CRYParticle*> parts;
109  fGen->genEvent(&parts);
110  for (unsigned int i=0; i<parts.size(); ++i) {
111 
112  // Take ownership of the particle from the vector
113  std::unique_ptr<CRYParticle> cryp(parts[i]);
114 
115  // Pull out the PDG code
116  int pdg = cryp->PDGid();
117 
118  // Get the energies of the particles
119  double ke = cryp->ke()*1.0E-3; // MeV to GeV conversion
120  if (ke<fEthresh) continue;
121 
122  double m = 0.; // in GeV
123 
124  static TDatabasePDG* pdgt = TDatabasePDG::Instance();
125  TParticlePDG* pdgp = pdgt->GetParticle(pdg);
126  if (pdgp) m = pdgp->Mass();
127 
128  double etot = ke + m;
129  double ptot = etot*etot-m*m;
130  if (ptot>0.0) ptot = sqrt(ptot);
131  else ptot = 0.0;
132 
133  // Sort out the momentum components. Remember that the NOvA
134  // frame has y up and z along the beam. So uvw -> zxy
135  double px = ptot * cryp->v();
136  double py = ptot * cryp->w();
137  double pz = ptot * cryp->u();
138 
139  // Particle start position. CRY distributes uniformly in x-y
140  // plane at fixed z, where z is the vertical direction. This
141  // requires some offsets and rotations to put the particles at
142  // the surface in the geometry as well as some rotations
143  // since the coordinate frame has y up and z along the
144  // beam.
145  double vx = cryp->y()*100.0;
146  double vy = cryp->z()*100.0 + surfaceY;
147  double vz = cryp->x()*100.0 + 0.5*detectorLength;
148  double t = cryp->t()-tstart + fToffset; // seconds
149  if(fSingleEventMode) t = fSampleTime*rantime; // seconds
150 
151  // Project backward to edge of world volume
152  double xyz[3] = { vx, vy, vz};
153  double xyzo[3];
154  double dxyz[3] = {-px, -py, -pz};
155  double x1 = 0.;
156  double x2 = 0.;
157  double y1 = 0.;
158  double y2 = 0.;
159  double z1 = 0.;
160  double z2 = 0.;
161  this->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
162 
163  MF_LOG_DEBUG("CRYHelper") << xyz[0] << " " << xyz[1] << " " << xyz[2] << " "
164  << x1 << " " << x2 << " "
165  << y1 << " " << y2 << " "
166  << z1 << " " << z2;
167 
168  this->ProjectToBoxEdge(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo);
169 
170  vx = xyzo[0];
171  vy = xyzo[1];
172  vz = xyzo[2];
173 
174  // Boiler plate...
175  int istatus = 1;
176  int imother1 = kCosmicRayGenerator;
177 
178  // Push the particle onto the stack
179  std::string primary("primary");
180 
181  particlespushed=true;
182  simb::MCParticle p(idctr,
183  pdg,
184  primary,
185  imother1,
186  m,
187  istatus);
188  TLorentzVector pos(vx,vy,vz,t*1e9);// time needs to be in ns to match GENIE, etc
189  TLorentzVector mom(px,py,pz,etot);
190  p.AddTrajectoryPoint(pos,mom);
191 
192  mctruth.Add(p);
193  ++idctr;
194  } // Loop on particles in event
195 
196  // Check if we're done with this time sample
197  // note that now requiring npart==1 in singlevent mode.
198 
199  if (fGen->timeSimulated()-tstart > fSampleTime ||
200  (fSingleEventMode && particlespushed )
201  ) break;
202  } // Loop on events simulated
203 
204  mctruth.SetOrigin(simb::kCosmicRay);
205 
207  if (w) *w = 1.0;
208  return fGen->timeSimulated()-tstart;
209 
210  }
double fSampleTime
Amount of time to sample (seconds)
Definition: CRYHelper.h:58
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
Definition: CRYHelper.cxx:225
Float_t y1[n_points_granero]
Definition: compare.C:5
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
Float_t x1[n_points_granero]
Definition: compare.C:5
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double &xlo, double &xhi, double &ylo, double &yhi, double &zlo, double &zhi, double xyzout[])
Definition: CRYHelper.cxx:264
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
Float_t y2[n_points_geant4]
Definition: compare.C:26
bool fSingleEventMode
flag to turn on producing only a single cosmic ray
Definition: CRYHelper.h:67
double fEthresh
Cut on kinetic energy (GeV)
Definition: CRYHelper.h:60
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
#define MF_LOG_DEBUG(id)
Float_t x2[n_points_geant4]
Definition: compare.C:26
Float_t w
Definition: plot.C:20
double fToffset
Shift in time of particles (s)
Definition: CRYHelper.h:59
Cosmic rays.
Definition: MCTruth.h:24
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, x1, x2, y1, and y2.

Referenced by Sample().

228  {
229  const TGeoShape* s = gGeoManager->GetVolume(fWorldVolume.c_str())->GetShape();
230  if(!s)
231  throw cet::exception("CRYHelper") << "No TGeoShape found for world volume";
232 
233  if (xlo || xhi) {
234  double x1, x2;
235  s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
236  }
237  if (ylo || yhi) {
238  double y1, y2;
239  s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
240  }
241  if (zlo || zhi) {
242  double z1, z2;
243  s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
244  }
245  }// end of WorldBox
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
Float_t y2[n_points_geant4]
Definition: compare.C:26
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61
Float_t x2[n_points_geant4]
Definition: compare.C:26
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

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: