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 "/cvmfs/nova.opensciencegrid.org/externals/nutools/v3_01_06/source/nutools/EventGeneratorBase/CRY/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 gen_flatrecord::config, febshutoff_auto::datapath, exit(), fAltitude, fGen, CLHEP::HepRandomEngine::flat(), fLatitude, fSetup, fSubBoxL, fhicl::ParameterSet::get(), cet::getenv(), evgb::RNGWrapper< T >::set(), and string.

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 
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
Definition: config.py:1
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
virtual double flat()=0
static void set(T *object, double(T::*func)(void))
Definition: CRYHelper.h:84
std::string getenv(std::string const &name)
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
exit(0)
double fToffset
Shift in time of particles (s)
Definition: CRYHelper.h:59
enum BeamMode string
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 268 of file CRYHelper.cxx.

References d, dx, dy, dz, fBoxDelta, and MECModelEnuComparisons::i.

Referenced by Sample().

274  {
275  // make the world box slightly smaller so that the projection to
276  // the edge avoids possible rounding errors later on with Geant4
277  xlo += fBoxDelta;
278  xhi -= fBoxDelta;
279  ylo += fBoxDelta;
280  yhi -= fBoxDelta;
281  zlo += fBoxDelta;
282  zhi -= fBoxDelta;
283 
284  // Make sure we're inside the box!
285  if(xyz[0] < xlo || xyz[0] > xhi ||
286  xyz[1] < ylo || xyz[1] > yhi ||
287  xyz[2] < zlo || xyz[2] > zhi)
288  throw cet::exception("CRYHelper") << "Projection to edge is outside"
289  << " bounds of world box:\n "
290  << "\tx: " << xyz[0] << " ("
291  << xlo << "," << xhi << ")\n"
292  << "\ty: " << xyz[1] << " ("
293  << ylo << "," << yhi << ")\n"
294  << "\tz: " << xyz[2] << " ("
295  << zlo << "," << zhi << ")";
296 
297  // Compute the distances to the x/y/z walls
298  double dx = 99.E99;
299  double dy = 99.E99;
300  double dz = 99.E99;
301  if (dxyz[0] > 0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
302  else if (dxyz[0] < 0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
303  if (dxyz[1] > 0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
304  else if (dxyz[1] < 0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
305  if (dxyz[2] > 0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
306  else if (dxyz[2] < 0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
307 
308  // Choose the shortest distance
309  double d = 0.0;
310  if (dx < dy && dx < dz) d = dx;
311  else if (dy < dz && dy < dx) d = dy;
312  else if (dz < dx && dz < dy) d = dz;
313 
314  // Make the step
315  for (int i = 0; i < 3; ++i) {
316  xyzout[i] = xyz[i] + dxyz[i]*d;
317  }
318  }
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double dy[NP][NC]
double fBoxDelta
Definition: CRYHelper.h:65
double dx[NP][NC]
double dz[NP][NC]
Float_t d
Definition: plot.C:236
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, cet::getenv(), MECModelEnuComparisons::i, simb::kCosmicRay, evgb::kCosmicRayGenerator, simb::kCRY, LOG_DEBUG, m, submit_concat_project::parts, make_root_from_grid_output::pdg, elec2geo::pos, ProjectToBoxEdge(), simb::MCTruth::SetGeneratorInfo(), simb::MCTruth::SetOrigin(), std::sqrt(), string, confusionMatrixTree::t, registry_explorer::v, WorldBox(), x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by evgen::CosmicsGen::produce(), and evgen::CosmicPionGen::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  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  std::string cryVersion;
206  if (auto v = std::getenv("CRY_VERSION"))
207  cryVersion = v;
208  mctruth.SetGeneratorInfo(simb::Generator_t::kCRY, cryVersion, {});
209 
210  /// \todo Check if this time slice passes selection criteria
211  if (w) *w = 1.0;
212  return fGen->timeSimulated()-tstart;
213 
214  }
double fSampleTime
Amount of time to sample (seconds)
Definition: CRYHelper.h:58
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
Definition: CRYHelper.cxx:229
Float_t y1[n_points_granero]
Definition: compare.C:5
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
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:268
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
std::string getenv(std::string const &name)
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 SetGeneratorInfo(simb::Generator_t generator, const std::string &genVersion, const std::unordered_map< std::string, std::string > &genConfig)
Definition: MCTruth.h:82
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
enum BeamMode string
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 229 of file CRYHelper.cxx.

References fWorldVolume, x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by Sample().

232  {
233  const TGeoShape* s = gGeoManager->GetVolume(fWorldVolume.c_str())->GetShape();
234  if(!s)
235  throw cet::exception("CRYHelper") << "No TGeoShape found for world volume";
236 
237  if (xlo || xhi) {
238  double x1, x2;
239  s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
240  }
241  if (ylo || yhi) {
242  double y1, y2;
243  s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
244  }
245  if (zlo || zhi) {
246  double z1, z2;
247  s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
248  }
249  }// end of WorldBox
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const XML_Char * s
Definition: expat.h:262
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61

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: