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