Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::CosmicsGen Class Reference

A module to produce cosmic ray events. More...

Inheritance diagram for evgen::CosmicsGen:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 CosmicsGen (fhicl::ParameterSet const &pset)
 
virtual ~CosmicsGen ()
 
void produce (art::Event &evt)
 
virtual void beginRun (art::Run &run)
 
virtual void beginSubRun (art::SubRun &run)
 
virtual void endSubRun (art::SubRun &run)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Member Functions

void DetectorBigBoxCut (simb::MCTruth &truth)
 
void ProjectCosmicsToSurface (simb::MCTruth &truth)
 
void ProjectMuonsToDetectorBigBox (simb::MCTruth &truth)
 
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. More...
 
void Enhance (simb::MCTruth &truth, bool *skipThisSet, double random_number)
 
bool isIntersectTheBox (const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
 
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[])
 
void ProjectToBox (const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[], double &dlout)
 

Private Attributes

evgb::CRYHelperfCRYHelp
 CRY generator object. More...
 
bool fIsBigBoxUsed
 
bool fBreakMCTruth
 
double fExposure
 Livetime this subrun, in seconds. More...
 
bool fGenSingleEvents
 generate one cosmic per record More...
 
double fProjectedZOffset
 How far above the BigBox ND muons should be started. More...
 
bool fEnhance
 Flag to enhace high energy particles. More...
 
int fEnhancePDG
 PDG of enhanced particles (neutron, proton or photon) More...
 
int fCycle
 MC cycle generation number. More...
 

Detailed Description

A module to produce cosmic ray events.

Definition at line 43 of file CosmicsGen_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

evgen::CosmicsGen::CosmicsGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 131 of file CosmicsGen_module.cc.

References art::EngineCreator::createEngine(), fCRYHelp, fEnhance, fGenSingleEvents, fhicl::ParameterSet::get(), art::RandomNumberGenerator::getEngine(), evgb::GetRandomNumberSeed(), art::EngineCreator::rng(), and seed.

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  }
double fExposure
Livetime this subrun, in seconds.
Interface to the CRY cosmic-ray generator.
Definition: CRYHelper.h:26
double fProjectedZOffset
How far above the BigBox ND muons should be started.
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
base_engine_t & createEngine(seed_t seed)
bool fEnhance
Flag to enhace high energy particles.
int fCycle
MC cycle generation number.
evgb::CRYHelper * fCRYHelp
CRY generator object.
base_engine_t & getEngine() const
unsigned int seed
Definition: runWimpSim.h:102
int fEnhancePDG
PDG of enhanced particles (neutron, proton or photon)
bool fGenSingleEvents
generate one cosmic per record
evgen::CosmicsGen::~CosmicsGen ( )
virtual

Definition at line 168 of file CosmicsGen_module.cc.

References fCRYHelp.

169  {
170  if(fCRYHelp) delete fCRYHelp;
171  }
evgb::CRYHelper * fCRYHelp
CRY generator object.

Member Function Documentation

void evgen::CosmicsGen::beginRun ( art::Run run)
virtual

Reimplemented from art::EDProducer.

Definition at line 174 of file CosmicsGen_module.cc.

References geo::GeometryBase::DetId(), geo::GeometryBase::ExtractGDML(), geo::GeometryBase::FileBaseName(), and art::Run::put().

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  }
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
Helper for AttenCurve.
Definition: Path.h:10
std::string FileBaseName() const
void evgen::CosmicsGen::beginSubRun ( art::SubRun run)
virtual

Reimplemented from art::EDProducer.

Definition at line 190 of file CosmicsGen_module.cc.

References fExposure.

191  {
192  fExposure = 0;
193  }
double fExposure
Livetime this subrun, in seconds.
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
void evgen::CosmicsGen::DetectorBigBoxCut ( simb::MCTruth truth)
private

Method to cut the events not going through the DetectorBigBox, which is defined as the box surrounding the DetectorEncosureBox by adding an additional length in each of the dimensions. The additional length is defined in detectrorbigbox.xml

Definition at line 359 of file CosmicsGen_module.cc.

References simb::MCTruth::Add(), geo::GeometryBase::DetectorBigBox(), geo::GeometryBase::DetId(), simb::MCParticle::EndX(), simb::MCParticle::EndY(), simb::MCParticle::EndZ(), simb::MCTruth::GetParticle(), ipart, isIntersectTheBox(), novadaq::cnv::kNEARDET, simb::MCParticle::Mass(), simb::MCParticle::Momentum(), simb::MCParticle::Mother(), npart, simb::MCTruth::NParticles(), simb::MCParticle::P(), part, simb::MCParticle::PdgCode(), simb::MCParticle::Process(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::TrackId(), x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by produce().

359  {
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  }
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
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
double Px(const int i=0) const
Definition: MCParticle.h:229
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
double EndY() const
Definition: MCParticle.h:226
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
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
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
double Pz(const int i=0) const
Definition: MCParticle.h:231
bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Event generator information.
Definition: MCTruth.h:32
double EndX() const
Definition: MCParticle.h:225
Helper for AttenCurve.
Definition: Path.h:10
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
void evgen::CosmicsGen::endSubRun ( art::SubRun run)
virtual

Reimplemented from art::EDProducer.

Definition at line 196 of file CosmicsGen_module.cc.

References fCycle, fEnhance, fExposure, fGenSingleEvents, art::SubRun::put(), and sd().

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  }
double fExposure
Livetime this subrun, in seconds.
bool fEnhance
Flag to enhace high energy particles.
int fCycle
MC cycle generation number.
double sd(Eigen::VectorXd x)
bool fGenSingleEvents
generate one cosmic per record
Helper for AttenCurve.
Definition: Path.h:10
void evgen::CosmicsGen::Enhance ( simb::MCTruth truth,
bool *  skipThisSet,
double  random_number 
)
private

Method to produce and enhanced sample with a richer population of high energy particles with PDG set in fEnhancePDG. This is done at this step (before Geant) to minimize running time. It only decides to keep or throw out complete events based on their most interesting particle so as to not interfere with the CRYHelper process. Currently set for neutrons, photons or protons independently.

Definition at line 276 of file CosmicsGen_module.cc.

References abs(), simb::MCTruth::Add(), energy, fEnhancePDG, simb::MCTruth::GetParticle(), ipart, simb::kCosmicRay, simb::MCParticle::Mass(), simb::MCParticle::Momentum(), simb::MCParticle::Mother(), getGoodRuns4SAM::n, npart, simb::MCTruth::NParticles(), part, make_root_from_grid_output::pdg, simb::MCParticle::PdgCode(), simb::MCParticle::Position(), cet::pow(), simb::MCParticle::Process(), simb::MCTruth::SetOrigin(), simb::MCParticle::StatusCode(), and simb::MCParticle::TrackId().

Referenced by produce().

276  {
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  }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
Int_t ipart
Definition: Style.C:10
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
void abs(TH1 *hist)
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
double energy
Definition: plottest35.C:25
int fEnhancePDG
PDG of enhanced particles (neutron, proton or photon)
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
Event generator information.
Definition: MCTruth.h:32
Cosmic rays.
Definition: MCTruth.h:24
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
bool evgen::CosmicsGen::isIntersectTheBox ( const double  xyz[],
const double  dxyz[],
double  xlo,
double  xhi,
double  ylo,
double  yhi,
double  zlo,
double  zhi 
)
private

Definition at line 653 of file CosmicsGen_module.cc.

References dx, dy, dz, stan::math::fabs(), submit_syst::x, xhi, make_syst_table_plots::xlo, submit_syst::y, yhi, ylo, and test::z.

Referenced by DetectorBigBoxCut().

657  {
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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
z
Definition: test.py:28
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void evgen::CosmicsGen::produce ( art::Event evt)
virtual

Implements art::EDProducer.

Definition at line 212 of file CosmicsGen_module.cc.

References DetectorBigBoxCut(), geo::GeometryBase::DetLength(), Enhance(), fBreakMCTruth, fCRYHelp, fEnhance, fExposure, fIsBigBoxUsed, CLHEP::HepRandomEngine::flat(), art::RandomNumberGenerator::getEngine(), simb::kCosmicRay, ProjectCosmicsToSurface(), ProjectMuonsToDetectorBigBox(), art::Event::put(), art::EngineCreator::rng(), evgb::CRYHelper::Sample(), SeparateMCTruth(), simb::MCTruth::SetOrigin(), and geo::GeometryBase::SurfaceY().

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  }
double fExposure
Livetime this subrun, in seconds.
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
double DetLength() const
bool fEnhance
Flag to enhace high energy particles.
void ProjectMuonsToDetectorBigBox(simb::MCTruth &truth)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
evgb::CRYHelper * fCRYHelp
CRY generator object.
base_engine_t & getEngine() const
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Definition: CRYHelper.cxx:97
virtual double flat()=0
void DetectorBigBoxCut(simb::MCTruth &truth)
double SurfaceY() const
void ProjectCosmicsToSurface(simb::MCTruth &truth)
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
Helper for AttenCurve.
Definition: Path.h:10
void Enhance(simb::MCTruth &truth, bool *skipThisSet, double random_number)
Cosmic rays.
Definition: MCTruth.h:24
void evgen::CosmicsGen::ProjectCosmicsToSurface ( simb::MCTruth truth)
private

Method to project the cosmic particles back to the surface for the near detector. CRYHelper is setup such that it projects CRY particles, which generated at the surface up to the end of the World volume. This is because ndos, for instance, is above the surface. CRYHelper is standard (in the EventGeneratorBase), used by other experiments. So, we don't want to modify it. But we do project cosmics back to the surface for the near detector (maybe should do the same for the far???)

Definition at line 440 of file CosmicsGen_module.cc.

References simb::MCTruth::Add(), geo::GeometryBase::DetId(), simb::MCParticle::E(), simb::MCTruth::GetParticle(), ipart, novadaq::cnv::kNEARDET, simb::MCParticle::Mass(), simb::MCParticle::Mother(), npart, simb::MCTruth::NParticles(), simb::MCParticle::P(), simb::MCParticle::PdgCode(), elec2geo::pos, simb::MCParticle::Process(), ProjectToSurface(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), simb::MCParticle::StatusCode(), geo::GeometryBase::SurfaceY(), simb::MCParticle::T(), simb::MCParticle::TrackId(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), geo::GeometryBase::WorldBox(), x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by produce().

440  {
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  }
double E(const int i=0) const
Definition: MCParticle.h:232
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
Int_t ipart
Definition: Style.C:10
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
double Px(const int i=0) const
Definition: MCParticle.h:229
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
int TrackId() const
Definition: MCParticle.h:209
Int_t npart
Definition: Style.C:7
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
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
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
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[])
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
double Pz(const int i=0) const
Definition: MCParticle.h:231
double Vz(const int i=0) const
Definition: MCParticle.h:222
Event generator information.
Definition: MCTruth.h:32
Helper for AttenCurve.
Definition: Path.h:10
double Vy(const int i=0) const
Definition: MCParticle.h:221
void evgen::CosmicsGen::ProjectMuonsToDetectorBigBox ( simb::MCTruth truth)
private

Method to project muons generated at the near detector to the DetectorBigBox. This is to get rid of the muon propagation through the rock, since the secondary particles would probably not reach the detector. Also we account for the muon's energy loss. After the muon reaches the DetectorBigBox, we start propagating it through Geant.

Definition at line 514 of file CosmicsGen_module.cc.

References abs(), simb::MCTruth::Add(), geo::GeometryBase::DetectorBigBox(), geo::GeometryBase::DetId(), simb::MCParticle::E(), epsilon, stan::math::exp(), fIsBigBoxUsed, fProjectedZOffset, simb::MCTruth::GetParticle(), ipart, novadaq::cnv::kNEARDET, m, simb::MCParticle::Mass(), simb::MCParticle::Mother(), npart, simb::MCTruth::NParticles(), simb::MCParticle::P(), simb::MCParticle::PdgCode(), elec2geo::pos, simb::MCParticle::Process(), ProjectToBox(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), std::sqrt(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::TrackId(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), X, x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by produce().

514  {
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  }
double E(const int i=0) const
Definition: MCParticle.h:232
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
Int_t ipart
Definition: Style.C:10
Float_t y1[n_points_granero]
Definition: compare.C:5
double fProjectedZOffset
How far above the BigBox ND muons should be started.
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
double Mass() const
Definition: MCParticle.h:238
double Px(const int i=0) const
Definition: MCParticle.h:229
::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
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
int TrackId() const
Definition: MCParticle.h:209
Int_t npart
Definition: Style.C:7
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
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
double Vx(const int i=0) const
Definition: MCParticle.h:220
double epsilon
double Pz(const int i=0) const
Definition: MCParticle.h:231
double Vz(const int i=0) const
Definition: MCParticle.h:222
Event generator information.
Definition: MCTruth.h:32
Float_t X
Definition: plot.C:38
Helper for AttenCurve.
Definition: Path.h:10
double Vy(const int i=0) const
Definition: MCParticle.h:221
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
void evgen::CosmicsGen::ProjectToBox ( const double  xyz[],
const double  dxyz[],
double  xlo,
double  xhi,
double  ylo,
double  yhi,
double  zlo,
double  zhi,
double  xyzout[],
double &  dlout 
)
private

Definition at line 817 of file CosmicsGen_module.cc.

References DEFINE_ART_MODULE(), dx, dy, dz, stan::math::fabs(), submit_syst::x, xhi, make_syst_table_plots::xlo, submit_syst::y, yhi, ylo, and test::z.

Referenced by ProjectMuonsToDetectorBigBox().

822  {
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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
z
Definition: test.py:28
void evgen::CosmicsGen::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[] 
)
private

Project Back to the surface

Parameters
xyz- The starting x,y,z location. Must be inside box.
dxyz- Direction vector
xlo- Low edge of box in x
xhi- High edge of box in x
ylo- Low edge of box in y
yhi- High edge of box in y
zlo- Low edge of box in z
zhi- High edge of box in z
surface_y- High edge of box in z
xyzout- On output, the position on the surface

Note: It should be safe to use the same array for input and output.

Definition at line 782 of file CosmicsGen_module.cc.

References ana::assert(), geo::GeometryBase::DetId(), geom(), MECModelEnuComparisons::i, and novadaq::cnv::kNEARDET.

Referenced by ProjectCosmicsToSurface().

788  {
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  }
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
void evgen::CosmicsGen::SeparateMCTruth ( simb::MCTruth truth,
std::vector< simb::MCTruth > *  truthcol 
)
private

Method to take the MCTruth object and beak it into pieces for each individual Cosmic Ray.

Definition at line 622 of file CosmicsGen_module.cc.

References simb::MCTruth::Add(), simb::MCTruth::GetParticle(), ipart, simb::MCParticle::Mass(), simb::MCParticle::Momentum(), simb::MCParticle::Mother(), npart, simb::MCTruth::NParticles(), part, simb::MCParticle::PdgCode(), simb::MCParticle::Position(), simb::MCParticle::Process(), simb::MCParticle::StatusCode(), simb::MCParticle::TrackId(), and simb::MCParticle::Weight().

Referenced by produce().

623  {
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  }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
Int_t ipart
Definition: Style.C:10
double Weight() const
Definition: MCParticle.h:253
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
Event generator information.
Definition: MCTruth.h:32
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

bool evgen::CosmicsGen::fBreakMCTruth
private

Definition at line 116 of file CosmicsGen_module.cc.

Referenced by produce().

evgb::CRYHelper* evgen::CosmicsGen::fCRYHelp
private

CRY generator object.

Definition at line 112 of file CosmicsGen_module.cc.

Referenced by CosmicsGen(), produce(), and ~CosmicsGen().

int evgen::CosmicsGen::fCycle
private

MC cycle generation number.

Definition at line 122 of file CosmicsGen_module.cc.

Referenced by endSubRun().

bool evgen::CosmicsGen::fEnhance
private

Flag to enhace high energy particles.

Definition at line 120 of file CosmicsGen_module.cc.

Referenced by CosmicsGen(), endSubRun(), and produce().

int evgen::CosmicsGen::fEnhancePDG
private

PDG of enhanced particles (neutron, proton or photon)

Definition at line 121 of file CosmicsGen_module.cc.

Referenced by Enhance().

double evgen::CosmicsGen::fExposure
private

Livetime this subrun, in seconds.

Definition at line 117 of file CosmicsGen_module.cc.

Referenced by beginSubRun(), endSubRun(), and produce().

bool evgen::CosmicsGen::fGenSingleEvents
private

generate one cosmic per record

Definition at line 118 of file CosmicsGen_module.cc.

Referenced by CosmicsGen(), and endSubRun().

bool evgen::CosmicsGen::fIsBigBoxUsed
private

Do we need to use the BigBox cut? The cosmic rays must go through the DetectorBigBox. Otherwise, don't store them

Definition at line 113 of file CosmicsGen_module.cc.

Referenced by produce(), and ProjectMuonsToDetectorBigBox().

double evgen::CosmicsGen::fProjectedZOffset
private

How far above the BigBox ND muons should be started.

Definition at line 119 of file CosmicsGen_module.cc.

Referenced by ProjectMuonsToDetectorBigBox().


The documentation for this class was generated from the following file: