CosmicsGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Generator for cosmic-rays
3 /// \author brebel@fnal.gov
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
6 
7 extern "C" {
8 #include <getopt.h>
9 }
10 #include <cmath> // for fabs() etc
11 #include <memory>
12 
13 // ROOT includes
14 #include "TH1.h"
15 #include "TH2.h"
16 
17 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
26 
27 // NOvA includes
28 
31 #include "Geometry/Geometry.h"
32 #include "NovaDAQConventions/DAQConventions.h"
36 #include "SummaryData/RunData.h"
37 #include "SummaryData/SubRunData.h"
38 
39 
40 namespace evgen {
41 
42  /// A module to produce cosmic ray events
43  class CosmicsGen : public art::EDProducer {
44  public:
45  explicit CosmicsGen(fhicl::ParameterSet const& pset);
46  virtual ~CosmicsGen();
47 
48  void produce(art::Event& evt);
49  virtual void beginRun(art::Run& run);
50  virtual void beginSubRun(art::SubRun& run);
51  virtual void endSubRun(art::SubRun& run);
52 
53  private:
54 
55  /// Method to cut the events not going through the DetectorBigBox,
56  /// which is defined as the box surrounding the DetectorEncosureBox
57  /// by adding an additional length in each of the dimensions.
58  /// The additional length is defined in detectrorbigbox.xml
59  void DetectorBigBoxCut(simb::MCTruth& truth);
60 
61 
62  /// Method to project the cosmic particles back to the surface for the near detector.
63  /// CRYHelper is setup such that it projects CRY particles, which generated
64  /// at the surface up to the end of the World volume.
65  /// This is because ndos, for instance, is above the surface.
66  /// CRYHelper is standard (in the EventGeneratorBase), used by other experiments.
67  /// So, we don't want to modify it.
68  /// But we do project cosmics back to the surface for the near detector
69  /// (maybe should do the same for the far???)
71 
72  /// Method to project muons generated at the near detector to the DetectorBigBox.
73  /// This is to get rid of the muon propagation through the rock, since
74  /// the secondary particles would probably not reach the detector.
75  /// Also we account for the muon's energy loss.
76  /// After the muon reaches the DetectorBigBox, we start propagating it through Geant.
78 
79  /// Method to take the MCTruth object and beak it into pieces for each individual Cosmic Ray
80  void SeparateMCTruth(simb::MCTruth& truth, std::vector<simb::MCTruth>* truthcol);
81 
82  /// Method to produce and enhanced sample with a richer population of high energy particles
83  /// with PDG set in fEnhancePDG. This is done at this step (before Geant) to minimize running
84  /// time. It only decides to keep or throw out complete events based on their most interesting
85  /// particle so as to not interfere with the CRYHelper process.
86  /// Currently set for neutrons, photons or protons independently.
87  void Enhance(simb::MCTruth& truth, bool* skipThisSet, double random_number);
88 
89  bool isIntersectTheBox(const double xyz[],
90  const double dxyz[],
91  double xlo, double xhi,
92  double ylo, double yhi,
93  double zlo, double zhi);
94 
95  void ProjectToSurface (const double xyz[],
96  const double dxyz[],
97  double xlo, double xhi,
98  double ylo, double yhi,
99  double zlo, double zhi,
100  const double surface_y,
101  double xyzout[]);
102 
103  void ProjectToBox (const double xyz[],
104  const double dxyz[],
105  double xlo, double xhi,
106  double ylo, double yhi,
107  double zlo, double zhi,
108  double xyzout[], double& dlout);
109 
110  private:
111 
112  evgb::CRYHelper* fCRYHelp; ///< CRY generator object
113  bool fIsBigBoxUsed; ///< Do we need to use the BigBox cut?
114  ///< The cosmic rays must go through the DetectorBigBox.
115  ///< Otherwise, don't store them
117  double fExposure; ///< Livetime this subrun, in seconds
118  bool fGenSingleEvents; ///< generate one cosmic per record
119  double fProjectedZOffset; ///< How far above the BigBox ND muons should be started
120  bool fEnhance; ///< Flag to enhace high energy particles
121  int fEnhancePDG; ///< PDG of enhanced particles (neutron, proton or photon)
122  int fCycle; ///< MC cycle generation number.
123  };
124 }
125 
126 ////////////////////////////////////////////////////////////////////////
127 
128 namespace evgen{
129 
130  //____________________________________________________________________________
132  : fIsBigBoxUsed (pset.get< bool >("IsBigBoxUsed"))
133  , fBreakMCTruth (pset.get< bool >("BreakMCTruth"))
134  , fExposure (0)
135  , fGenSingleEvents (pset.get< bool >("GenSingleEvents"))
136  , fProjectedZOffset (pset.get< double >("ProjectedZOffset"))
137  , fEnhance (pset.get< bool >("Enhance"))
138  , fEnhancePDG (pset.get< int >("EnhancePDG"))
139  , fCycle (pset.get< int >("Cycle", 0))
140  {
141 
142  // Create an ART Random Number engine
143  int seed = pset.get< int >("Seed", evgb::GetRandomNumberSeed());
144  createEngine(seed);
145 
146  // get the random number generator service and make some CLHEP generators
148  CLHEP::HepRandomEngine& engine = rng->getEngine();
149 
150  // Create a new CRYHelper
151  fCRYHelp = new evgb::CRYHelper(pset, engine);
152 
153  // Producers
154  produces< std::vector<simb::MCTruth> >();
155  produces< sumdata::SubRunData, art::InSubRun >();
156  produces< sumdata::RunData, art::InRun >();
157 
158  // This information is misleading when we're generating in normal mode, with
159  // many cosmics per spill. CosmicExposureInfo gets it right. But for single
160  // event mode we need it to do the accounting for overlays.
161  // We also need accounting of exposure for hadron enhanced samples.
162  if(fGenSingleEvents || fEnhance){
163  produces< sumdata::CosmicExposure, art::InSubRun >();
164  }
165  }
166 
167  //____________________________________________________________________________
169  {
170  if(fCRYHelp) delete fCRYHelp;
171  }
172 
173  //____________________________________________________________________________
175  {
176 
177  // grab the geometry object to see what geometry we are using
179 
180  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
181  geo->FileBaseName(),
182  geo->ExtractGDML()));
183 
184  run.put(std::move(runcol));
185 
186  return;
187  }
188 
189  //____________________________________________________________________________
191  {
192  fExposure = 0;
193  }
194 
195  //____________________________________________________________________________
197  {
198  if(fGenSingleEvents || fEnhance){
199  std::unique_ptr<sumdata::CosmicExposure> ce(new sumdata::CosmicExposure);
200  ce->totlivetime = fExposure;
201  subrun.put(std::move(ce));
202  }
203 
204  // store the cycle information
206  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
207  subrun.put(std::move(sd));
208 
209  }
210 
211  //____________________________________________________________________________
213  {
214  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
215 
218  double rantime = rng->getEngine().flat();
219 
220  simb::MCTruth truth;
222 
223  double sampletime;
224  double total_sampletime = 0.;
225 
226  // While skippedSet was implemented to regenerate the Truth object multiple
227  // times for enhanced samples but is skipped for normal CRY.
228  bool skippedSet = true;
229  int total_skipped = -1;
230  while ( skippedSet ){
231  total_skipped++;
232 
233  rantime = rng->getEngine().flat();
234  sampletime = fCRYHelp->Sample(truth, geo->SurfaceY(), geo->DetLength(), 0, rantime);
235 
236  // 1. Project the cosmics rays down to the surface (ND only)
237  this->ProjectCosmicsToSurface(truth);
238 
239  // 2. DetectorBigBox cut
240  if(fIsBigBoxUsed) this->DetectorBigBoxCut(truth);
241 
242  // 3. ND only: Rewind cosmics back to the surface after adjustment by 2.
243  this->ProjectCosmicsToSurface(truth);
244 
245  total_sampletime += sampletime;
246 
247  // * Option to produce an enhanced sample.
248  if(!fEnhance) break;
249  double rand_num = rng->getEngine().flat();
250  this->Enhance( truth, &skippedSet, rand_num);
251 
252  }
253 
254  // 4. ND only: Now project muons to the Big Box
255  this->ProjectMuonsToDetectorBigBox(truth);
256 
257  //5. Break MCTruth into individual truth objects for each cosmic
258  if(fBreakMCTruth){
259  this->SeparateMCTruth(truth, truthcol.get());
260  }
261  else {
262  truthcol->push_back(truth);
263  }
264 
265  evt.put(std::move(truthcol));
266 
267  fExposure += total_sampletime;
268 
269  return;
270 
271  }
272 
273  ////////////////////////////////////////////////////////////////////////
274  // Method to enhance high E Neutrons, Protons or Photons by deciding
275  // whether to keep or skip the current Truth based on a probability.
276  void CosmicsGen::Enhance(simb::MCTruth& truth, bool* skipThisSet, double random_number){
277 
278  // Parameters from the energy distribution fits docdb 13199
279  // TO DO: Make these fcl parameters
280 
281  double nNeutron = -5.18;
282  double E0Neutron = 4;
283 
284  double nPhoton = -2.63;
285  double E0Photon = 3;
286 
287  double nProton = -3.91;
288  double E0Proton = 4;
289 
290  // Enhance with weight w = k(E/E0)^n
291  double n, E0;
292 
293  if( fEnhancePDG == 2112 ){
294  n = nNeutron;
295  E0 = E0Neutron;
296  }
297  if( fEnhancePDG == 2212 ){
298  n = nProton;
299  E0 = E0Proton;
300  }
301  if( fEnhancePDG == 22 ){
302  n = nPhoton;
303  E0 = E0Photon;
304  }
305 
306  double event_score = 0.001; // probability to keep evt, better for high E fEnhancePDGs
307  double event_w = 1.; // assign weight to relate these events to normal CRY stats
308 
309  int nparticles = truth.NParticles();
310  for ( int part_Idx = 0; part_Idx < nparticles; ++part_Idx){
311  simb::MCParticle particle = truth.GetParticle(part_Idx);
312  int pdg = abs(particle.PdgCode());
313  // Calculate event score by looking at interesting particles
314 
315  if ( pdg == fEnhancePDG ){
316  const TLorentzVector& fourMomentum = particle.Momentum();
317  double energy = fourMomentum.E();
318  // Keep the largest particle score for the event
319  double score = pow((E0/energy),n);
320  if ( score > event_score ) event_score = score;
321  // Keep all events with score >= 1
322  if ( event_score > 1 ) break;
323  }
324 
325  }// All particles in truth
326  if ( event_score>1 ) event_score = 1;
327  event_w = 1.0/event_score;
328 
329  simb::MCTruth Enhancedtruth;
330  Enhancedtruth.SetOrigin(simb::kCosmicRay);
331  if ( random_number<event_score ) {
332  *skipThisSet = false;
333  int npart = truth.NParticles();
334  for( int ipart = 0; ipart < npart; ++ipart ){
335  simb::MCParticle particle = truth.GetParticle(ipart);
336 
337  simb::MCParticle part(particle.TrackId(),
338  particle.PdgCode(),
339  particle.Process(),
340  particle.Mother(),
341  particle.Mass(),
342  particle.StatusCode());
343  part.AddTrajectoryPoint(particle.Position(),
344  particle.Momentum());
345  part.SetWeight(event_w);
346 
347  Enhancedtruth.Add(part);
348  }
349  }
350  // Enhancedtruth remains empty unless we keep the event.
351  truth=Enhancedtruth;
352  }
353 
354  ////////////////////////////////////////////////////////////////////////
355  // Method to cut the events not going through the DetectorBigBox,
356  // which is defined as the box surrounding the DetectorEncosureBox
357  // by adding an additional length in each of the dimensions.
358  // The additional length is defined in detectrorbigbox.xml
360  simb::MCTruth mctruth_new;
361 
362  double x1=0;
363  double x2=0;
364  double y1=0;
365  double y2=0;
366  double z1=0;
367  double z2=0;
368 
370 
371  geo->DetectorBigBox(&x1, &x2, &y1, &y2, &z1, &z2);
372 
373  if ( x1==0 && x2==0 &&
374  y1==0 && y2==0 &&
375  z1==0 && z2==0 ) {
376  throw cet::exception("NoGeometryBoxCoords")
377  << "No Geometry Box Coordinates Set\n"
378  << __FILE__ << ":" << __LINE__ << "\n";
379  return;
380  }
381 
382  // Loop over the particles stored by CRY
383  int npart = truth.NParticles();
384  for(int ipart = 0; ipart < npart; ++ipart){
385  simb::MCParticle particle = truth.GetParticle(ipart);
386 
387  double px = particle.Px() / particle.P();
388  double py = particle.Py() / particle.P();
389  double pz = particle.Pz() / particle.P();
390 
391  double vx = particle.EndX();
392  double vy = particle.EndY();
393  double vz = particle.EndZ();
394 
395  // For the near detector, move all the cosmics from their position in
396  // the window, straight down, to 1cm above the detector big box. This
397  // ensures that a much larger fraction of them pass the IntersectTheBox
398  // cut below. Subsequently the ProjectCosmicsToSurface call will
399  // backtrack them up to the surface. The flux remains correct throughout
400  // this scheme due to the homogeneity of the flux in space and time.
401  // It does get the correlations wrong, but the chances of two correlated
402  // events making it to the ND are remote anyway.
403  if(geo->DetId() == novadaq::cnv::kNEARDET){
404  vy = y2+1;
405  }
406 
407  double xyz[3] = { vx, vy, vz};
408  double dxyz[3] = { px, py, pz};
409 
410 
411  // If intersecting the DetectorBigBox, add particle
412  //currently, using the methods in CosmicsGen code until the geometry methods can be validated
413  if(this->isIntersectTheBox(xyz, dxyz, x1, x2, y1, y2, z1, z2)){
414  simb::MCParticle part(particle.TrackId(),
415  particle.PdgCode(),
416  particle.Process(),
417  particle.Mother(),
418  particle.Mass(),
419  particle.StatusCode());
420  part.AddTrajectoryPoint(TLorentzVector(vx, vy, vz, particle.T()),
421  particle.Momentum());
422 
423  mctruth_new.Add(part);
424  }
425  }// end of loop over particles
426 
427  // reassign the MCTruth
428  truth = mctruth_new;
429  }
430 
431  ////////////////////////////////////////////////////////////////////////
432  // Method to project the cosmic particles back to the surface for the near detector.
433  // CRYHelper is setup such that it projects CRY particles, which generated
434  // at the surface up to the end of the World volume.
435  // This is because ndos, for instance, is above the surface.
436  // CRYHelper is standard (in the EventGeneratorBase), used by other experiments.
437  // So, we don't want to modify it.
438  // But we do project cosmics back to the surface for the near detector
439  // (maybe should do the same for the far???)
441 
443  if(geo->DetId() != novadaq::cnv::kNEARDET) return;
444 
445  simb::MCTruth truth_after_surface_projection;
446 
447  double x1=0;
448  double x2=0;
449  double y1=0;
450  double y2=0;
451  double z1=0;
452  double z2=0;
453 
454  // Get the World box
455  geo->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
456 
457  if ( x1==0 && x2==0 &&
458  y1==0 && y2==0 &&
459  z1==0 && z2==0 ) {
460  throw cet::exception("NoWorldBoxCoords")
461  << "No World Box Coordinates Set\n"
462  << __FILE__ << ":" << __LINE__ << "\n";
463  return;
464  }
465 
466  // Loop over the particles
467  int npart = truth.NParticles();
468  for(int ipart = 0; ipart < npart; ++ipart){
469  simb::MCParticle particle = truth.GetParticle(ipart);
470 
471  // make a new particle from this particle, where the first point
472  // is at the surface
473  simb::MCParticle patsurface(particle.TrackId(),
474  particle.PdgCode(),
475  particle.Process(),
476  particle.Mother(),
477  particle.Mass(),
478  particle.StatusCode());
479 
480  double px = particle.Px() / particle.P();
481  double py = particle.Py() / particle.P();
482  double pz = particle.Pz() / particle.P();
483 
484  double vx = particle.Vx();
485  double vy = particle.Vy();
486  double vz = particle.Vz();
487 
488  double xyz[3] = { vx, vy, vz};
489  double dxyz[3] = { px, py, pz};
490 
491  double xyzo[3];
492 
493  this->ProjectToSurface(xyz, dxyz, x1, x2, y1, y2, z1, z2, geo->SurfaceY(), xyzo);
494 
495  // reset vertex to be at the surface rather than the edge of the world
496  TLorentzVector pos(xyzo[0], xyzo[1], xyzo[2], particle.T());
497  TLorentzVector mom(particle.Px(), particle.Py(), particle.Pz(), particle.E());
498 
499  patsurface.AddTrajectoryPoint(pos, mom);
500 
501  truth_after_surface_projection.Add(patsurface);
502  }// end of loop over particles
503 
504  // redefine truth
505  truth = truth_after_surface_projection;
506  }
507 
508  ////////////////////////////////////////////////////////////////////////
509  // Method to project muons generated at the near detector to the DetectorBigBox.
510  // This is to get rid of the muon propagation through the rock, since
511  // the secondary particles would probably not reach the detector.
512  // Also we account for the muon's energy loss.
513  // After the muon reaches the DetectorBigBox, we start propagating it through Geant.
515 
517 
518  // Projecting muons only if BigBox cut was used
519  // Otherwise, the muon could miss the BigBox.
520  // We just don't want to deal with these kind of events.
521  if( !fIsBigBoxUsed || geo->DetId() != novadaq::cnv::kNEARDET) return;
522 
523  simb::MCTruth truth_after_bigbox_projection;
524 
525  double x1=0;
526  double x2=0;
527  double y1=0;
528  double y2=0;
529  double z1=0;
530  double z2=0;
531 
532  geo->DetectorBigBox(&x1, &x2, &y1, &y2, &z1, &z2);
533 
534  if ( x1==0 && x2==0 &&
535  y1==0 && y2==0 &&
536  z1==0 && z2==0 ) {
537  throw cet::exception("NoGeometryBoxCoords")
538  << "No Geometry Box Coordinates Set\n"
539  << __FILE__ << ":" << __LINE__ << "\n";
540  return;
541  }
542 
543  z2 += fProjectedZOffset;
544 
545  // Loop over the particles stored by CRY
546  int npart = truth.NParticles();
547  for(int ipart = 0; ipart < npart; ++ipart){
548 
549  simb::MCParticle particle = truth.GetParticle(ipart);
550 
551  //make a new particle from this particle, where the first point has the correct
552  //momentum after the rock if it is a muon
553  simb::MCParticle pafterrock(particle.TrackId(),
554  particle.PdgCode(),
555  particle.Process(),
556  particle.Mother(),
557  particle.Mass(),
558  particle.StatusCode());
559 
560  // Project only mu+ and mu-
561  if(abs(particle.PdgCode()) == 13 ){
562 
563  double px = particle.Px() / particle.P();
564  double py = particle.Py() / particle.P();
565  double pz = particle.Pz() / particle.P();
566 
567  double vx = particle.Vx();
568  double vy = particle.Vy();
569  double vz = particle.Vz();
570 
571  double xyz[3] = { vx, vy, vz};
572  double dxyz[3] = { px, py, pz};
573 
574  double xyzo[3];
575  double dlout=0;
576  this->ProjectToBox(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo, dlout);
577 
578  // In GeV
579  double m = particle.Mass();
580  double etot = particle.E();
581  double ke = etot - m;
582 
583  double X = 2.5 * dlout; //g/sm^2 these magic numbers come from Gaisser,
584  double ksi = 250000.; //g/cm^2 Cosmic Rays and Particle Physics equation (6.17)
585  double epsilon = 500.0; //GeV
586 
587  double KE_after_rock = (ke + epsilon)*exp(-X/ksi) - epsilon;
588  if(KE_after_rock < 0)continue;
589  double E_after_rock = KE_after_rock + m;
590 
591  double P_after_rock = sqrt(E_after_rock*E_after_rock - m*m);
592  double Px = particle.Px()/particle.P() * P_after_rock;
593  double Py = particle.Py()/particle.P() * P_after_rock;
594  double Pz = particle.Pz()/particle.P() * P_after_rock;
595 
596  TLorentzVector pos(xyzo[0], xyzo[1], xyzo[2], particle.T());
597  TLorentzVector mom(Px, Py, Pz, E_after_rock);
598 
599  pafterrock.AddTrajectoryPoint(pos, mom);
600 
601 
602  }// end of if that particle is mu- or mu+
603  else{
604  // Don't waste time tracking non-muons through geant when they likely
605  // won't reach the detector. As well as the waste of time, they'd break
606  // the empty spill detection.
607  continue;
608  }
609 
610  // Add the particle into the new stack
611  truth_after_bigbox_projection.Add(pafterrock);
612 
613  }// end of loop over particles
614 
615  // redefine the truth
616  truth = truth_after_bigbox_projection;
617  }
618 
619  /////////////////////////////////////////////////////
620  /////////////////////////////////////////////////////
621  //Method to separate MCTruth into separate objects for each particle
623  std::vector<simb::MCTruth>* truthcol){
624 
625  // Loop over the particles stored by CRY
626  int npart = truth.NParticles();
627  for(int ipart = 0; ipart < npart; ++ipart){
628  simb::MCParticle particle = truth.GetParticle(ipart);
629 
630  simb::MCTruth ptruth;
631 
632  simb::MCParticle part(particle.TrackId(),
633  particle.PdgCode(),
634  particle.Process(),
635  particle.Mother(),
636  particle.Mass(),
637  particle.StatusCode());
638  part.AddTrajectoryPoint(particle.Position(),
639  particle.Momentum());
640 
641  part.SetWeight(particle.Weight());
642 
643  ptruth.Add(part);
644 
645  truthcol->push_back(ptruth);
646  }
647 
648  }
649 
650  //------------------------------------------------------------------------------
651  // Suplementary methods
652  ////////////////////////////////////////////////////////////////////////
653  bool CosmicsGen::isIntersectTheBox(const double xyz[],
654  const double dxyz[],
655  double xlo, double xhi,
656  double ylo, double yhi,
657  double zlo, double zhi) {
658 
659  // If inside the box, obviously intesected
660  if(xyz[0]>=xlo && xyz[0]<=xhi &&
661  xyz[1]>=ylo && xyz[1]<=yhi &&
662  xyz[2]>=zlo && xyz[2]<=zhi) return true;
663 
664  // So, now we know that the particle is outside the box
665 
666  double dx = 0., dy = 0., dz = 0.;
667 
668  // Checking intersection with 6 planes
669 
670  // 1. Check intersection with the upper plane
671  dy = xyz[1] - yhi;
672  // Check whether the track going from above and down or from below and up.
673  // Otherwise not going to intersect this plane
674  if( (dy > 0 && dxyz[1] < 0) ||
675  (dy < 0 && dxyz[1] > 0) ){
676  double dl = fabs( dy / dxyz[1] );
677  double x = xyz[0] + dxyz[0] * dl;
678  double z = xyz[2] + dxyz[2] * dl;
679 
680  // is it inside the square?
681  if( x>=xlo && x<=xhi &&
682  z>=zlo && z<=zhi) return true;
683  }
684 
685  // 2. Check intersection with the lower plane
686  dy = xyz[1] - ylo;
687  // Check whether the track going from above and down or from below and up.
688  // Otherwise not going to intersect this plane
689  if( (dy > 0 && dxyz[1] < 0) ||
690  (dy < 0 && dxyz[1] > 0) ){
691  double dl = fabs( dy / dxyz[1] );
692  double x = xyz[0] + dxyz[0] * dl;
693  double z = xyz[2] + dxyz[2] * dl;
694 
695  // is it inside the square?
696  if( x>=xlo && x<=xhi &&
697  z>=zlo && z<=zhi) return true;
698  }
699 
700 
701  // 3. Check intersection with the right plane
702  dz = xyz[2] - zhi;
703  // Check whether the track going from right part to the left or from left part to right.
704  // Otherwise not going to intersect this plane
705  if( (dz > 0 && dxyz[2] < 0) ||
706  (dz < 0 && dxyz[2] > 0) ){
707  double dl = fabs( dz / dxyz[2] );
708  double x = xyz[0] + dxyz[0] * dl;
709  double y = xyz[1] + dxyz[1] * dl;
710 
711  // is it inside the square?
712  if( x>=xlo && x<=xhi &&
713  y>=ylo && y<=yhi) return true;
714  }
715 
716  // 4. Check intersection with the left plane
717  dz = xyz[2] - zlo;
718  // Check whether the track going from right part to the left or from left part to right.
719  // Otherwise not going to intersect this plane
720  if( (dz > 0 && dxyz[2] < 0) ||
721  (dz < 0 && dxyz[2] > 0) ){
722  double dl = fabs( dz / dxyz[2] );
723  double x = xyz[0] + dxyz[0] * dl;
724  double y = xyz[1] + dxyz[1] * dl;
725 
726  // is it inside the square?
727  if( x>=xlo && x<=xhi &&
728  y>=ylo && y<=yhi) return true;
729  }
730 
731  // 5. Check intersection with the far plane
732  dx = xyz[0] - xhi;
733  // Check whether the track going from far part toward us or from near part to right.
734  // Otherwise not going to intersect this plane
735  if( (dx > 0 && dxyz[0] < 0) ||
736  (dx < 0 && dxyz[0] > 0) ){
737  double dl = fabs( dx / dxyz[0] );
738  double y = xyz[1] + dxyz[1] * dl;
739  double z = xyz[2] + dxyz[2] * dl;
740 
741  // is it inside the square?
742  if( z>=zlo && z<=zhi &&
743  y>=ylo && y<=yhi) return true;
744  }
745 
746 
747  // 6. Check intersection with the near plane
748  dx = xyz[0] - xlo;
749  // Check whether the track going from far part toward us or from near part to right.
750  // Otherwise not going to intersect this plane
751  if( (dx > 0 && dxyz[0] < 0) ||
752  (dx < 0 && dxyz[0] > 0) ){
753  double dl = fabs( dx / dxyz[0] );
754  double y = xyz[1] + dxyz[1] * dl;
755  double z = xyz[2] + dxyz[2] * dl;
756 
757  // is it inside the square?
758  if( z>=zlo && z<=zhi &&
759  y>=ylo && y<=yhi) return true;
760  }
761 
762  return false;
763  }
764  //
765  ///----------------------------------------------------------------
766  /// Project Back to the surface
767  ///
768  /// \param xyz - The starting x,y,z location. Must be inside box.
769  /// \param dxyz - Direction vector
770  /// \param xlo - Low edge of box in x
771  /// \param xhi - High edge of box in x
772  /// \param ylo - Low edge of box in y
773  /// \param yhi - High edge of box in y
774  /// \param zlo - Low edge of box in z
775  /// \param zhi - High edge of box in z
776  /// \param surface_y - High edge of box in z
777  /// \param xyzout - On output, the position on the surface
778  ///
779  /// Note: It should be safe to use the same array for input and
780  /// output.
781  ///
782  void CosmicsGen::ProjectToSurface(const double xyz[],
783  const double dxyz[],
784  double xlo, double xhi,
785  double ylo, double yhi,
786  double zlo, double zhi,
787  const double surface_y,
788  double xyzout[]) {
790 
791  // Make sure that the cosmic muon is going down from above
792  // (except in ND where we're actually rewinding them to the surface).
793  assert(geom->DetId() == novadaq::cnv::kNEARDET ||
794  xyz[1] >= surface_y);
795  assert(dxyz[1] < 0);
796  // Also that the surface that we are projecting to is between ylo and yhi
797  assert(surface_y >= ylo && surface_y <= yhi);
798 
799  // Step length
800  double dl = - (xyz[1] - surface_y) / dxyz[1];
801 
802  // Make the step
803  for (int i=0; i<3; ++i)
804  xyzout[i] = xyz[i] + dxyz[i]*dl;
805 
806 
807  if(geom->DetId() != novadaq::cnv::kNEARDET){
808  // Make sure we're inside the box!
809  assert(xyzout[0]>=xlo && xyzout[0]<=xhi);
810  assert(xyzout[1]>=ylo && xyzout[1]<=yhi);
811  assert(xyzout[2]>=zlo && xyzout[2]<=zhi);
812  }
813  }
814 
815  ////////////////////////////////////////////////////////////////////////
816  ////////////////////////////////////////////////////////////////////////
817  void CosmicsGen::ProjectToBox (const double xyz[],
818  const double dxyz[],
819  double xlo, double xhi,
820  double ylo, double yhi,
821  double zlo, double zhi,
822  double xyzout[], double& dlout) {
823 
824 
825  // So, now we know that the particle is outside the box
826 
827  double dx, dy, dz;
828 
829  // Initialize to some big number
830  dlout = 999999999.;
831 
832  // Checking intersection with 6 planes
833  // 1. Check intersection with the upper plane
834  dy = xyz[1] - yhi;
835  // Check whether the track going from above and down or from below and up.
836  // Otherwise not going to intersect this plane
837  if( (dy > 0 && dxyz[1] < 0) ||
838  (dy < 0 && dxyz[1] > 0) ){
839  double dl = fabs( dy / dxyz[1] );
840  double x = xyz[0] + dxyz[0] * dl;
841  double z = xyz[2] + dxyz[2] * dl;
842 
843  // is it inside the square?
844  if( x>=xlo && x<=xhi &&
845  z>=zlo && z<=zhi &&
846  dl < dlout ) {
847  xyzout[0] = x;
848  xyzout[1] = yhi;
849  xyzout[2] = z;
850  dlout = dl;
851  }
852  }
853 
854  // 2. Check intersection with the lower plane
855  // Actually this cannot be the plane for the cosmics, because it would have to cross
856  // some other plane first. But what the heck, let's check
857  dy = xyz[1] - ylo;
858  // Check whether the track going from above and down or from below and up.
859  // Otherwise not going to intersect this plane
860  if( (dy > 0 && dxyz[1] < 0) ||
861  (dy < 0 && dxyz[1] > 0) ){
862  double dl = fabs( dy / dxyz[1] );
863  double x = xyz[0] + dxyz[0] * dl;
864  double z = xyz[2] + dxyz[2] * dl;
865 
866  // is it inside the square?
867  if( x>=xlo && x<=xhi &&
868  z>=zlo && z<=zhi &&
869  dl < dlout ) {
870  xyzout[0] = x;
871  xyzout[1] = ylo;
872  xyzout[2] = z;
873  dlout = dl;
874  }
875  }
876 
877 
878  // 3. Check intersection with the right plane
879  dz = xyz[2] - zhi;
880  // Check whether the track going from right part to the left or from left part to right.
881  // Otherwise not going to intersect this plane
882  if( (dz > 0 && dxyz[2] < 0) ||
883  (dz < 0 && dxyz[2] > 0) ){
884  double dl = fabs( dz / dxyz[2] );
885  double x = xyz[0] + dxyz[0] * dl;
886  double y = xyz[1] + dxyz[1] * dl;
887 
888  // is it inside the square?
889  if( x>=xlo && x<=xhi &&
890  y>=ylo && y<=yhi &&
891  dl < dlout ) {
892  xyzout[0] = x;
893  xyzout[1] = y;
894  xyzout[2] = zhi;
895  dlout = dl;
896  }
897  }
898 
899  // 4. Check intersection with the left plane
900  dz = xyz[2] - zlo;
901  // Check whether the track going from right part to the left or from left part to right.
902  // Otherwise not going to intersect this plane
903  if( (dz > 0 && dxyz[2] < 0) ||
904  (dz < 0 && dxyz[2] > 0) ){
905  double dl = fabs( dz / dxyz[2] );
906  double x = xyz[0] + dxyz[0] * dl;
907  double y = xyz[1] + dxyz[1] * dl;
908 
909  // is it inside the square?
910  if( x>=xlo && x<=xhi &&
911  y>=ylo && y<=yhi &&
912  dl < dlout ) {
913  xyzout[0] = x;
914  xyzout[1] = y;
915  xyzout[2] = zlo;
916  dlout = dl;
917  }
918  }
919 
920  // 5. Check intersection with the far plane
921  dx = xyz[0] - xhi;
922  // Check whether the track going from far part toward us or from near part to right.
923  // Otherwise not going to intersect this plane
924  if( (dx > 0 && dxyz[0] < 0) ||
925  (dx < 0 && dxyz[0] > 0) ){
926  double dl = fabs( dx / dxyz[0] );
927  double y = xyz[1] + dxyz[1] * dl;
928  double z = xyz[2] + dxyz[2] * dl;
929 
930  // is it inside the square?
931  if( z>=zlo && z<=zhi &&
932  y>=ylo && y<=yhi &&
933  dl < dlout ) {
934  xyzout[0] = xhi;
935  xyzout[1] = y;
936  xyzout[2] = z;
937  dlout = dl;
938  }
939  }
940 
941 
942  // 6. Check intersection with the near plane
943  dx = xyz[0] - xlo;
944  // Check whether the track going from far part toward us or from near part to right.
945  // Otherwise not going to intersect this plane
946  if( (dx > 0 && dxyz[0] < 0) ||
947  (dx < 0 && dxyz[0] > 0) ){
948  double dl = fabs( dx / dxyz[0] );
949  double y = xyz[1] + dxyz[1] * dl;
950  double z = xyz[2] + dxyz[2] * dl;
951 
952  // is it inside the square?
953  if( z>=zlo && z<=zhi &&
954  y>=ylo && y<=yhi &&
955  dl < dlout ) {
956  xyzout[0] = xlo;
957  xyzout[1] = y;
958  xyzout[2] = z;
959  dlout = dl;
960  }
961  }
962  }
963 
964 }// namespace
965 
966 /////////////////////////////////////////////////////
967 namespace evgen
968 {
970 }
double E(const int i=0) const
Definition: MCParticle.h:232
double fExposure
Livetime this subrun, in seconds.
Interface to the CRY cosmic-ray generator.
Definition: CRYHelper.h:26
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double EndZ() const
Definition: MCParticle.h:227
Int_t ipart
Definition: Style.C:10
double Weight() const
Definition: MCParticle.h:253
Float_t y1[n_points_granero]
Definition: compare.C:5
double fProjectedZOffset
How far above the BigBox ND muons should be started.
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
void ProjectToBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[], double &dlout)
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
T sqrt(T number)
Definition: d0nt_math.hpp:156
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
double Px(const int i=0) const
Definition: MCParticle.h:229
base_engine_t & createEngine(seed_t seed)
double DetLength() const
bool fEnhance
Flag to enhace high energy particles.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
Definition: MCParticle.h:214
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
double EndY() const
Definition: MCParticle.h:226
void ProjectMuonsToDetectorBigBox(simb::MCTruth &truth)
Definition: Run.h:31
int TrackId() const
Definition: MCParticle.h:209
int fCycle
MC cycle generation number.
CosmicsGen(fhicl::ParameterSet const &pset)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
evgb::CRYHelper * fCRYHelp
CRY generator object.
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Definition: CRYHelper.cxx:97
double dy[NP][NC]
double dx[NP][NC]
virtual double flat()=0
double dz[NP][NC]
unsigned int seed
Definition: runWimpSim.h:102
Interface to the CRY cosmic ray generator.
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
void DetectorBigBoxCut(simb::MCTruth &truth)
double P(const int i=0) const
Definition: MCParticle.h:233
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
int evt
double energy
Definition: plottest35.C:25
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
int fEnhancePDG
PDG of enhanced particles (neutron, proton or photon)
void ProjectToSurface(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, const double surface_y, double xyzout[])
z
Definition: test.py:28
Definition: run.py:1
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
double Vx(const int i=0) const
Definition: MCParticle.h:220
double SurfaceY() const
void produce(art::Event &evt)
bool fGenSingleEvents
generate one cosmic per record
ProductID put(std::unique_ptr< PROD > &&)
double epsilon
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void ProjectCosmicsToSurface(simb::MCTruth &truth)
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
virtual void endSubRun(art::SubRun &run)
double Vz(const int i=0) const
Definition: MCParticle.h:222
assert(nhit_max >=nhit_nbins)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
virtual void beginSubRun(art::SubRun &run)
A module to produce cosmic ray events.
bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
virtual void beginRun(art::Run &run)
void SeparateMCTruth(simb::MCTruth &truth, std::vector< simb::MCTruth > *truthcol)
Method to take the MCTruth object and beak it into pieces for each individual Cosmic Ray...
Event generator information.
Definition: MCTruth.h:32
double EndX() const
Definition: MCParticle.h:225
Float_t X
Definition: plot.C:38
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
void Enhance(simb::MCTruth &truth, bool *skipThisSet, double random_number)
Cosmic rays.
Definition: MCTruth.h:24
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
std::string FileBaseName() const