CRYHelper.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file CRYHelper.cxx
3 /// \brief Implementation of an interface to the CRY cosmic-ray generator.
4 ///
5 /// \version $Id: CRYHelper.cxx,v 1.27 2012-10-15 20:46:42 brebel Exp $
6 /// \author messier@indiana.edu
7 ////////////////////////////////////////////////////////////////////////
8 #include <cmath>
9 #include <iostream>
10 
11 // CRY include files
12 #include "CRYSetup.h"
13 #include "CRYParticle.h"
14 #include "CRYGenerator.h"
15 
16 // ROOT include files
17 #include "TDatabasePDG.h"
18 #include "TLorentzVector.h"
19 #include "TGeoManager.h"
20 
21 // Framework includes
23 #include "fhiclcpp/ParameterSet.h"
25 #include "cetlib_except/exception.h"
26 
27 // NuTools include files
32 
33 namespace evgb{
34 
35  //......................................................................
37  {
38  }
39 
40  //......................................................................
42  CLHEP::HepRandomEngine& engine,
43  std::string const& worldVol)
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  }
88 
89  //......................................................................
91  {
92  delete fGen;
93  delete fSetup;
94  }
95 
96  //......................................................................
97  double CRYHelper::Sample(simb::MCTruth& mctruth,
98  double const& surfaceY,
99  double const& detectorLength,
100  double* w,
101  double rantime)
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  }
215 
216  ///----------------------------------------------------------------
217  ///
218  /// Return the ranges of x,y and z for the "world volume" that the
219  /// entire geometry lives in. If any pointers are 0, then those
220  /// coordinates are ignored.
221  ///
222  /// \param xlo : On return, lower bound on x positions
223  /// \param xhi : On return, upper bound on x positions
224  /// \param ylo : On return, lower bound on y positions
225  /// \param yhi : On return, upper bound on y positions
226  /// \param zlo : On return, lower bound on z positions
227  /// \param zhi : On return, upper bound on z positions
228  ///
229  void CRYHelper::WorldBox(double* xlo, double* xhi,
230  double* ylo, double* yhi,
231  double* zlo, double* zhi) const
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
250 
251  ///----------------------------------------------------------------
252  /// Project along a direction from a particular starting point to the
253  /// edge of a box
254  ///
255  /// \param xyz - The starting x,y,z location. Must be inside box.
256  /// \param dxyz - Direction vector
257  /// \param xlo - Low edge of box in x
258  /// \param xhi - Low edge of box in x
259  /// \param ylo - Low edge of box in y
260  /// \param yhi - Low edge of box in y
261  /// \param zlo - Low edge of box in z
262  /// \param zhi - Low edge of box in z
263  /// \param xyzout - On output, the position at the box edge
264  ///
265  /// Note: It should be safe to use the same array for input and
266  /// output.
267  ///
268  void CRYHelper::ProjectToBoxEdge(const double xyz[],
269  const double dxyz[],
270  double &xlo, double &xhi,
271  double &ylo, double &yhi,
272  double &zlo, double &zhi,
273  double xyzout[])
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  }
319 }
320 ////////////////////////////////////////////////////////////////////////
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
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
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
CRYSetup * fSetup
CRY configuration.
Definition: CRYHelper.h:56
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
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::string fLatitude
Latitude of detector need space after value.
Definition: CRYHelper.h:62
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
const XML_Char * s
Definition: expat.h:262
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Definition: CRYHelper.cxx:97
double dy[NP][NC]
std::string fWorldVolume
Name of the world volume.
Definition: CRYHelper.h:61
double fBoxDelta
Definition: CRYHelper.h:65
double dx[NP][NC]
virtual double flat()=0
static void set(T *object, double(T::*func)(void))
Definition: CRYHelper.h:84
double dz[NP][NC]
Interface to the CRY cosmic ray generator.
std::string getenv(std::string const &name)
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fAltitude
Altitude of detector need space after value.
Definition: CRYHelper.h:63
Float_t d
Definition: plot.C:236
bool fSingleEventMode
flag to turn on producing only a single cosmic ray
Definition: CRYHelper.h:67
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
void SetGeneratorInfo(simb::Generator_t generator, const std::string &genVersion, const std::unordered_map< std::string, std::string > &genConfig)
Definition: MCTruth.h:82
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
exit(0)
Physics generators for neutrinos, cosmic rays, and others.
Definition: CRYHelper.cxx:33
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
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