Public Member Functions | Private Member Functions | Private Attributes | List of all members
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.

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 
83  fSetup->setRandomFunction(RNGWrapper<CLHEP::HepRandomEngine>::rng);
84 
85  fGen = new CRYGenerator(fSetup);
86 
87  }
evgb::CRYHelper::fSampleTime
double fSampleTime
Amount of time to sample (seconds)
Definition: CRYHelper.h:58
string
std::string string
Definition: nybbler.cc:12
evgb::CRYHelper::fSetup
CRYSetup * fSetup
CRY configuration.
Definition: CRYHelper.h:56
evgb::CRYHelper::fLatitude
std::string fLatitude
Latitude of detector need space after value.
Definition: CRYHelper.h:62
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
evgb::CRYHelper::fGen
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
evgb::CRYHelper::fWorldVolume
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61
evgb::CRYHelper::fBoxDelta
double fBoxDelta
Definition: CRYHelper.h:65
evgb::RNGWrapper::set
static void set(T *object, double(T::*func)(void))
Definition: CRYHelper.h:84
config
static Config * config
Definition: config.cpp:1054
cet::getenv
std::string getenv(std::string const &name)
Definition: getenv.cc:15
evgb::CRYHelper::fAltitude
std::string fAltitude
Altitude of detector need space after value.
Definition: CRYHelper.h:63
evgb::CRYHelper::fSingleEventMode
bool fSingleEventMode
flag to turn on producing only a single cosmic ray
Definition: CRYHelper.h:67
evgb::RNGWrapper::rng
static double rng(void)
Definition: CRYHelper.h:88
evgb::CRYHelper::fSubBoxL
std::string fSubBoxL
Length of subbox (m) need space after value.
Definition: CRYHelper.h:64
evgb::CRYHelper::fEthresh
double fEthresh
Cut on kinetic energy (GeV)
Definition: CRYHelper.h:60
make_pion_angle.flat
def flat(h_in)
Definition: make_pion_angle.py:38
evgb::CRYHelper::fToffset
double fToffset
Shift in time of particles (s)
Definition: CRYHelper.h:59
evgb::CRYHelper::~CRYHelper ( )

Definition at line 90 of file CRYHelper.cxx.

91  {
92  delete fGen;
93  delete fSetup;
94  }
evgb::CRYHelper::fSetup
CRYSetup * fSetup
CRY configuration.
Definition: CRYHelper.h:56
evgb::CRYHelper::fGen
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.

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  }
bump_copyright.d
string d
Definition: bump_copyright.py:70
evgb::CRYHelper::fBoxDelta
double fBoxDelta
Definition: CRYHelper.h:65
fhicl::exception
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.

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 
206  /// \todo Check if this time slice passes selection criteria
207  if (w) *w = 1.0;
208  return fGen->timeSimulated()-tstart;
209 
210  }
evgb::CRYHelper::fSampleTime
double fSampleTime
Amount of time to sample (seconds)
Definition: CRYHelper.h:58
evgb::CRYHelper::WorldBox
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
reco_momentum_tuples.py
py
Definition: reco_momentum_tuples.py:219
simb::MCTruth::SetOrigin
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
string
std::string string
Definition: nybbler.cc:12
evgb::CRYHelper::ProjectToBoxEdge
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
reco_momentum_tuples.px
px
Definition: reco_momentum_tuples.py:218
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
evgb::CRYHelper::fGen
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
test.p
p
Definition: test.py:223
evgb::CRYHelper::fSingleEventMode
bool fSingleEventMode
flag to turn on producing only a single cosmic ray
Definition: CRYHelper.h:67
cache_state.m
m
Definition: cache_state.py:415
evgb::CRYHelper::fEthresh
double fEthresh
Cut on kinetic energy (GeV)
Definition: CRYHelper.h:60
reco_momentum_tuples.pz
pz
Definition: reco_momentum_tuples.py:220
simb::MCTruth::Add
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
evgb::CRYHelper::fToffset
double fToffset
Shift in time of particles (s)
Definition: CRYHelper.h:59
simb::kCosmicRay
Cosmic rays.
Definition: MCTruth.h:24
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
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.

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
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
evgb::CRYHelper::fWorldVolume
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61
s
static QCString * s
Definition: config.cpp:1042
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82

Member Data Documentation

std::string evgb::CRYHelper::fAltitude
private

Altitude of detector need space after value.

Definition at line 63 of file CRYHelper.h.

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.

double evgb::CRYHelper::fEthresh
private

Cut on kinetic energy (GeV)

Definition at line 60 of file CRYHelper.h.

CRYGenerator* evgb::CRYHelper::fGen
private

The CRY generator.

Definition at line 57 of file CRYHelper.h.

std::string evgb::CRYHelper::fLatitude
private

Latitude of detector need space after value.

Definition at line 62 of file CRYHelper.h.

double evgb::CRYHelper::fSampleTime
private

Amount of time to sample (seconds)

Definition at line 58 of file CRYHelper.h.

CRYSetup* evgb::CRYHelper::fSetup
private

CRY configuration.

Definition at line 56 of file CRYHelper.h.

bool evgb::CRYHelper::fSingleEventMode
private

flag to turn on producing only a single cosmic ray

Definition at line 67 of file CRYHelper.h.

std::string evgb::CRYHelper::fSubBoxL
private

Length of subbox (m) need space after value.

Definition at line 64 of file CRYHelper.h.

double evgb::CRYHelper::fToffset
private

Shift in time of particles (s)

Definition at line 59 of file CRYHelper.h.

std::string evgb::CRYHelper::fWorldVolume
private

Name of the world volume.

Definition at line 61 of file CRYHelper.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11