LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CRYHelper.cxx
Go to the documentation of this file.
1 #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 
81  RNGWrapper<CLHEP::HepRandomEngine>::set(&engine, &CLHEP::HepRandomEngine::flat);
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  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  }
211 
225  void CRYHelper::WorldBox(double* xlo, double* xhi,
226  double* ylo, double* yhi,
227  double* zlo, double* zhi) const
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
246 
264  void CRYHelper::ProjectToBoxEdge(const double xyz[],
265  const double dxyz[],
266  double &xlo, double &xhi,
267  double &ylo, double &yhi,
268  double &zlo, double &zhi,
269  double xyzout[])
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  }
315 }
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
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:258
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
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:264
STL namespace.
std::string fLatitude
Latitude of detector need space after value.
Definition: CRYHelper.h:62
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Particle class.
CRYGenerator * fGen
The CRY generator.
Definition: CRYHelper.h:57
Float_t y2[n_points_geant4]
Definition: compare.C:26
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Definition: CRYHelper.cxx:97
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
Interface to the CRY cosmic ray generator.
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::string fAltitude
Altitude of detector need space after value.
Definition: CRYHelper.h:63
Float_t d
Definition: plot.C:235
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 Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
#define MF_LOG_DEBUG(id)
Physics generators for neutrinos, cosmic rays, and others.
Definition: CRYHelper.cxx:33
Event generator information.
Definition: MCTruth.h:32
Float_t x2[n_points_geant4]
Definition: compare.C:26
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
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Cosmic rays.
Definition: MCTruth.h:24