All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
evgen::ldm::MeVPrtlGen Class Reference
Inheritance diagram for evgen::ldm::MeVPrtlGen:

Public Member Functions

 MeVPrtlGen (fhicl::ParameterSet const &p)
 
 MeVPrtlGen (MeVPrtlGen const &)=delete
 
 MeVPrtlGen (MeVPrtlGen &&)=delete
 
MeVPrtlGen & operator= (MeVPrtlGen const &)=delete
 
MeVPrtlGen & operator= (MeVPrtlGen &&)=delete
 
void produce (art::Event &e) override
 
void beginRun (art::Run &run) override
 
void endSubRun (art::SubRun &sr) override
 
bool Deweight (double &weight, double &max_weight)
 
 ~MeVPrtlGen () noexcept
 

Private Attributes

bool fProduce
 
bool fAnaOutput
 
bool fDoDeweight
 
double fSubRunPOT
 
std::unique_ptr
< evgen::ldm::IMesonGen > 		fGenTool
 
std::unique_ptr
< evgen::ldm::IMeVPrtlFlux > 		fFluxTool
 
std::unique_ptr
< evgen::ldm::IRayTrace > 		fRayTool
 
std::unique_ptr
< evgen::ldm::IMeVPrtlDecay > 		fDecayTool
 
double fGenMaxWeight
 
double fFluxMaxWeight
 
double fRayMaxWeight
 
double fDecayMaxWeight
 
double fRayDecayMaxWeight
 
TTree * fTree
 
MeVPrtlTruth * fMeVPrtl
 
CLHEP::HepJamesRandom * fEngine
 
std::array< uint64_t, 4 > fNCalls
 
std::array< double, 4 > fNTime
 

Detailed Description

Definition at line 58 of file MeVPrtlGen_module.cc.

Constructor & Destructor Documentation

evgen::ldm::MeVPrtlGen::MeVPrtlGen ( fhicl::ParameterSet const &  p)
explicit

Definition at line 116 of file MeVPrtlGen_module.cc.

117  : EDProducer{p}
118 {
119  // nullify pointers
120  fMeVPrtl = NULL;
121  fEngine = NULL;
122 
123  fProduce = p.get<bool>("Produce", true);
124  fAnaOutput = p.get<bool>("AnaOutput", false);
125 
126  fDoDeweight = p.get<bool>("Deweight", false);
127  fSubRunPOT = 0.;
128 
129  // Update constants
130  if (p.has_key("Constants")) Constants::Configure(p.get<fhicl::ParameterSet>("Constants"));
131 
132  // bring in the tools
133  fGenTool = art::make_tool<IMesonGen>(p.get<fhicl::ParameterSet>("MesonGen"));
134  fFluxTool = art::make_tool<IMeVPrtlFlux>(p.get<fhicl::ParameterSet>("Flux"));
135  fRayTool = art::make_tool<IRayTrace>(p.get<fhicl::ParameterSet>("RayTrace"));
136  fDecayTool = art::make_tool<IMeVPrtlDecay>(p.get<fhicl::ParameterSet>("Decay"));
137 
138  fGenMaxWeight = fGenTool->MaxWeight();
139  std::cout << "Gen max weight: " << fGenMaxWeight << std::endl;
140 
141  fFluxMaxWeight = fFluxTool->MaxWeight();
142  std::cout << "Flux max weight: " << fFluxMaxWeight << std::endl;
143 
144  fRayMaxWeight = fRayTool->MaxWeight();
145  std::cout << "Ray max weight: " << fRayMaxWeight << std::endl;
146 
147  fDecayMaxWeight = fDecayTool->MaxWeight();
148  std::cout << "Decay max weight: " << fDecayMaxWeight << std::endl;
149 
150  fRayDecayMaxWeight = fDecayMaxWeight*fRayMaxWeight;
151 
152  if (fProduce) {
153  // All the standard generator outputs
154  produces< std::vector<simb::MCTruth> >();
155  produces< std::vector<simb::MCFlux> >();
156  produces< sumdata::RunData, art::InRun >();
157  produces< sumdata::POTSummary, art::InSubRun >();
158  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
159  produces< std::vector<sim::BeamGateInfo> >();
160 
161  // also save info pertinent to the scalar
162  produces< std::vector<evgen::ldm::MeVPrtlTruth> >();
163  }
164 
165  // setup the random number engine
166  art::ServiceHandle<rndm::NuRandomService> seedSvc;
167  fEngine = new CLHEP::HepJamesRandom;
168  seedSvc->registerEngine(rndm::NuRandomService::CLHEPengineSeeder(fEngine), "MeVPrtlGen");
169 
170  if (fAnaOutput) {
171  art::ServiceHandle<art::TFileService> tfs;
172  fTree = tfs->make<TTree>("mevprtl_gen", "mevprtl_gen");
173  fMeVPrtl = new evgen::ldm::MeVPrtlTruth();
174  fTree->Branch("mevprtl", &fMeVPrtl);
175  }
176 
177  fNCalls = {0, 0, 0, 0};
178  fNTime = {0., 0., 0., 0.};
179 }
p
pdgs p
Definition: selectors.fcl:22
evgen::ldm::MeVPrtlGen::fEngine
CLHEP::HepJamesRandom * fEngine
Definition: MeVPrtlGen_module.cc:110
evgen::produces< sumdata::RunData, art::InRun >
produces< sumdata::RunData, art::InRun >()
evgen::ldm::MeVPrtlGen::fNTime
std::array< double, 4 > fNTime
Definition: MeVPrtlGen_module.cc:113
evgen::ldm::MeVPrtlGen::fDecayTool
std::unique_ptr< evgen::ldm::IMeVPrtlDecay > fDecayTool
Definition: MeVPrtlGen_module.cc:99
evgen::ldm::Constants::Configure
static void Configure(const fhicl::ParameterSet &p)
Definition: Constants.cpp:60
evgen::ldm::MeVPrtlGen::fRayTool
std::unique_ptr< evgen::ldm::IRayTrace > fRayTool
Definition: MeVPrtlGen_module.cc:98
evgen::ldm::MeVPrtlGen::fNCalls
std::array< uint64_t, 4 > fNCalls
Definition: MeVPrtlGen_module.cc:112
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
evgen::ldm::MeVPrtlGen::fGenTool
std::unique_ptr< evgen::ldm::IMesonGen > fGenTool
Definition: MeVPrtlGen_module.cc:96
evgen::ldm::MeVPrtlGen::fFluxTool
std::unique_ptr< evgen::ldm::IMeVPrtlFlux > fFluxTool
Definition: MeVPrtlGen_module.cc:97
evgen::ldm::MeVPrtlGen::MeVPrtlGen ( MeVPrtlGen const &  )
delete
evgen::ldm::MeVPrtlGen::MeVPrtlGen ( MeVPrtlGen &&  )
delete
evgen::ldm::MeVPrtlGen::~MeVPrtlGen ( )
inlinenoexcept

Definition at line 79 of file MeVPrtlGen_module.cc.

79  {
80  std::cout << "GenTool called (" << fNCalls[0] << ") times. Total duration (" << fNTime[0] << ") ms. Duration per call (" << (fNTime[0] / fNCalls[0]) << ") ms.\n";
81  std::cout << "FlxTool called (" << fNCalls[1] << ") times. Total duration (" << fNTime[1] << ") ms. Duration per call (" << (fNTime[1] / fNCalls[1]) << ") ms.\n";
82  std::cout << "RayTool called (" << fNCalls[2] << ") times. Total duration (" << fNTime[2] << ") ms. Duration per call (" << (fNTime[2] / fNCalls[2]) << ") ms.\n";
83  std::cout << "DcyTool called (" << fNCalls[3] << ") times. Total duration (" << fNTime[3] << ") ms. Duration per call (" << (fNTime[3] / fNCalls[3]) << ") ms.\n";
84 
85  if (fEngine) delete fEngine;
86  if (fMeVPrtl) delete fMeVPrtl;
87  }
evgen::ldm::MeVPrtlGen::fEngine
CLHEP::HepJamesRandom * fEngine
Definition: MeVPrtlGen_module.cc:110
evgen::ldm::MeVPrtlGen::fNTime
std::array< double, 4 > fNTime
Definition: MeVPrtlGen_module.cc:113
evgen::ldm::MeVPrtlGen::fNCalls
std::array< uint64_t, 4 > fNCalls
Definition: MeVPrtlGen_module.cc:112
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10

Member Function Documentation

void evgen::ldm::MeVPrtlGen::beginRun ( art::Run &  run)
override

Definition at line 181 of file MeVPrtlGen_module.cc.

181  {
182  art::ServiceHandle<geo::Geometry const> geo;
183  if (fProduce) run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
184 }
bool evgen::ldm::MeVPrtlGen::Deweight ( double &  weight,
double &  max_weight 
)

Definition at line 196 of file MeVPrtlGen_module.cc.

196  {
197  if (!fDoDeweight || max_weight < 0) { // don't do deweighting procedure
198  return true;
199  }
200 
201  // Guard against bad max weight
202  //
203  // There is some question of whether to crash or handle this gracefully.
204  // Note that there is some literature that proposes that doing rejection
205  // sampling with an adaptive maximum weight (as done below) is valid.
206  // https://www.jstor.org/stable/4140534?seq=1
207  //
208  // However, this is really only true for a larger N than what is typically done
209  // in a larsoft job (10-150). Still, it's probably best not to crash.
210  if (max_weight < weight) {
211  std::cerr << "ERROR: weight (" << weight << ") with max weight (" << max_weight << "). Reconfiguration needed.\n";
212  std::cout << "ERROR: weight (" << weight << ") with max weight (" << max_weight << "). Reconfiguration needed.\n";
213  std::cout << "Updating max_weight to new value!\n";
214  max_weight = weight;
215  return true;
216  }
217 
218  // do deweighting
219  double rand = CLHEP::RandFlat::shoot(fEngine, 0, max_weight);
220  double test = weight;
221 
222  // update the weight value
223  weight = max_weight;
224  return rand <= test;
225 }
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
evgen::ldm::MeVPrtlGen::fEngine
CLHEP::HepJamesRandom * fEngine
Definition: MeVPrtlGen_module.cc:110
test
BEGIN_PROLOG hitmakerfive clustermakerfour pfparticlemakerthree showermakertwo END_PROLOG hitmakerfive clustermakerfour pfparticlemakerthree sequence::inline_paths sequence::inline_paths sequence::inline_paths showermakers test
Definition: assnschainutils_test.fcl:19
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
void evgen::ldm::MeVPrtlGen::endSubRun ( art::SubRun &  sr)
override

Definition at line 186 of file MeVPrtlGen_module.cc.

186  {
187  auto p = std::make_unique<sumdata::POTSummary>();
188  p->totpot = fSubRunPOT;
189  p->totgoodpot = fSubRunPOT;
190 
191  if (fProduce) sr.put(std::move(p));
192 
193  fSubRunPOT = 0.;
194 }
p
pdgs p
Definition: selectors.fcl:22
MeVPrtlGen& evgen::ldm::MeVPrtlGen::operator= ( MeVPrtlGen const &  )
delete
MeVPrtlGen& evgen::ldm::MeVPrtlGen::operator= ( MeVPrtlGen &&  )
delete
void evgen::ldm::MeVPrtlGen::produce ( art::Event &  e)
override

Definition at line 227 of file MeVPrtlGen_module.cc.

228 {
229  std::unique_ptr<std::vector<simb::MCFlux>> mcfluxColl(new std::vector<simb::MCFlux>);
230  std::unique_ptr<std::vector<simb::MCTruth>> mctruthColl(new std::vector<simb::MCTruth>);
231  std::unique_ptr<art::Assns<simb::MCTruth, simb::MCFlux>> truth2fluxAssn(new art::Assns<simb::MCTruth, simb::MCFlux>);
232  std::unique_ptr<std::vector<sim::BeamGateInfo>> beamgateColl(new std::vector<sim::BeamGateInfo>);
233 
234  std::unique_ptr<std::vector<evgen::ldm::MeVPrtlTruth>> mevprtlColl(new std::vector<evgen::ldm::MeVPrtlTruth>);
235 
236  // TODO: pileup? For now, don't worry
237 
238  // get the next MeVPrtl Truth
239  while (1) {
240 
241  std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
242  simb::MCFlux kaon = fGenTool->GetNext();
243  fNCalls[0] ++;
244  std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
245  std::chrono::duration<double, std::milli> duration = t2 - t1;
246  fNTime[0] = duration.count();
247 
248  evgen::ldm::KaonParent kaonp(kaon);
249  bool is_kaon = kaonp.kaon_pdg != 0;
250 
251  // (void) is_kaon;
252  if (is_kaon) {
253  std::cout << "Flux is kaon (" << is_kaon << "). Weight: " << kaonp.weight << ". Produced with energy: " << kaonp.mom.E()
254  << " M=" << kaonp.mom.M() << " P=(" << kaonp.mom.Px() << ", " << kaonp.mom.Py() << ", " << kaonp.mom.Pz() << ") At: ("
255  << kaonp.pos.X() << ", " << kaonp.pos.Y() << ", " << kaonp.pos.Z() << ")" << std::endl;
256  }
257 
258  bool success;
259 
260  evgen::ldm::MeVPrtlFlux flux;
261  double flux_weight;
262 
263  fNCalls[1] ++;
264  t1 = std::chrono::high_resolution_clock::now();
265  success = fFluxTool->MakeFlux(kaon, flux, flux_weight) && Deweight(flux_weight, fFluxMaxWeight);
266  t2 = std::chrono::high_resolution_clock::now();
267  duration = t2 - t1;
268  fNTime[1] += duration.count();
269 
270  if (!success) continue;
271 
272  std::cout << "New flux. E=" << flux.mom.E() << " At: (" << flux.pos.X() << ", " << flux.pos.Y() << ", " << flux.pos.Z() << ")" << std::endl;
273  std::cout << "P=(" << flux.mom.Px() << ", " << flux.mom.Py() << ", " << flux.mom.Pz() << ")" << std::endl;
274  std::cout << "Flux weight: " << flux_weight << std::endl;
275 
276  std::array<TVector3, 2> intersection;
277  double ray_weight;
278 
279  fNCalls[2] ++;
280  t1 = std::chrono::high_resolution_clock::now();
281  success = fRayTool->IntersectDetector(flux, intersection, ray_weight);
282  t2 = std::chrono::high_resolution_clock::now();
283  duration = t2 - t1;
284  fNTime[2] += duration.count();
285 
286  if (!success) continue;
287  std::cout << "Ray weight: " << ray_weight << std::endl;
288 
289  evgen::ldm::MeVPrtlDecay decay;
290  double decay_weight;
291 
292  fNCalls[3] ++;
293  t1 = std::chrono::high_resolution_clock::now();
294  success = fDecayTool->Decay(flux, intersection[0], intersection[1], decay, decay_weight);
295  t2 = std::chrono::high_resolution_clock::now();
296  duration = t2 - t1;
297  fNTime[3] += duration.count();
298 
299  if (!success) continue;
300 
301  std::cout << "Decay weight: " << decay_weight << std::endl;
302 
303  // Deweight the ray and decay weights together because they have some anti-correlation
304  double ray_decay_weight = ray_weight * decay_weight;
305  success = Deweight(ray_decay_weight, fRayDecayMaxWeight);
306 
307  if (!success) continue;
308 
309  std::cout << "RayDecay weight: " << ray_decay_weight << std::endl;
310 
311  std::cout << "PASSED!\n";
312 
313  // get the POT
314  double thisPOT = fGenTool->GetPOT();
315 
316  // if we are de-weighting, then the scaling all gets put into the POT variable
317  if (fDoDeweight) {
318  thisPOT = thisPOT / (flux_weight * ray_decay_weight);
319  flux_weight = 1.;
320  ray_weight = 1.;
321  decay_weight = 1.;
322  ray_decay_weight = 1.;
323  }
324 
325  fSubRunPOT += thisPOT;
326 
327  // build the output objects
328  evgen::ldm::MeVPrtlTruth mevprtl_truth(flux, decay,
329  intersection,
330  flux_weight,
331  1.,
332  ray_decay_weight,
333  thisPOT
334  );
335 
336  mevprtlColl->push_back(mevprtl_truth);
337 
338  simb::MCTruth mctruth;
339 
340  // Add the "Neutrino" as the 0th MCParticle
341  // This hopefully (???) won't do anything too bad and will give us
342  // the chance to use the neutrino energy in other things
343  simb::MCParticle fakenu(0, kaon.fntype, "primary", -1, 0, -1/* don't track */);
344  fakenu.AddTrajectoryPoint(mevprtl_truth.decay_pos, TLorentzVector(0, 0, flux.equiv_enu, flux.equiv_enu));
345  mctruth.Add(fakenu);
346  mctruth.SetNeutrino(-1, -1, -1, -1, -1, -1,
347  -1., -1., -1., -1.);
348 
349  for (unsigned i_d = 0; i_d < mevprtl_truth.daughter_mom.size(); i_d++) {
350  TLorentzVector daughter4p(mevprtl_truth.daughter_mom[i_d], mevprtl_truth.daughter_e[i_d]);
351  simb::MCParticle d(0, mevprtl_truth.daughter_pdg[i_d], "primary", -1, daughter4p.M());
352  d.AddTrajectoryPoint(mevprtl_truth.decay_pos, daughter4p);
353  mctruth.Add(d);
354  }
355 
356  // TODO:
357  //
358  // Flux systematic uncertainties are often evaluated on the neutrino energy.
359  // We need to figure out how to translate this for the case of heavy particles
360  // with different production kinematics. For now, we could save the neutrino energy
361  // so that the flux uncertainties "work" at some level.
362  //
363  // However, the existing MCTruth object has no way to "just" set a neutrino energy.
364  // This is __very__ very annoying.
365 
366  mctruthColl->push_back(mctruth);
367  mcfluxColl->push_back(kaon);
368 
369  // Make the associations only if we are producing stuff
370  // Otherwise this crashes
371  if (fProduce) {
372  art::PtrMaker<simb::MCFlux> MCFluxPtrMaker {evt};
373  art::PtrMaker<simb::MCTruth> MCTruthPtrMaker {evt};
374 
375  art::Ptr<simb::MCTruth> MCTruthPtr = MCTruthPtrMaker(mctruthColl->size() - 1);
376  art::Ptr<simb::MCFlux> MCFluxPtr = MCFluxPtrMaker(mcfluxColl->size() - 1);
377  truth2fluxAssn->addSingle(MCTruthPtr, MCFluxPtr);
378  }
379 
380  // TODO: implement for real
381  sim::BeamGateInfo gate;
382  beamgateColl->push_back(gate);
383 
384  if (fAnaOutput) {
385  *fMeVPrtl = mevprtl_truth;
386  fTree->Fill();
387  }
388 
389  break;
390  }
391 
392  if (fProduce) {
393  evt.put(std::move(mctruthColl));
394  evt.put(std::move(mcfluxColl));
395  evt.put(std::move(truth2fluxAssn));
396  evt.put(std::move(beamgateColl));
397 
398  evt.put(std::move(mevprtlColl));
399  }
400 
401 }
evgen::ldm::MeVPrtlFlux::mom
TLorentzVector mom
Definition: MeVPrtlFlux.h:14
evgen::ldm::MeVPrtlGen::fNTime
std::array< double, 4 > fNTime
Definition: MeVPrtlGen_module.cc:113
evgen::ldm::MeVPrtlGen::fDecayTool
std::unique_ptr< evgen::ldm::IMeVPrtlDecay > fDecayTool
Definition: MeVPrtlGen_module.cc:99
evgen::ldm::MeVPrtlGen::Deweight
bool Deweight(double &weight, double &max_weight)
Definition: MeVPrtlGen_module.cc:196
evgen::ldm::MeVPrtlGen::fRayTool
std::unique_ptr< evgen::ldm::IRayTrace > fRayTool
Definition: MeVPrtlGen_module.cc:98
evgen::ldm::MeVPrtlGen::fNCalls
std::array< uint64_t, 4 > fNCalls
Definition: MeVPrtlGen_module.cc:112
evgen::ldm::MeVPrtlFlux::pos
TLorentzVector pos
Definition: MeVPrtlFlux.h:11
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
evgen::ldm::MeVPrtlGen::fGenTool
std::unique_ptr< evgen::ldm::IMesonGen > fGenTool
Definition: MeVPrtlGen_module.cc:96
evgen::ldm::MeVPrtlGen::fFluxTool
std::unique_ptr< evgen::ldm::IMeVPrtlFlux > fFluxTool
Definition: MeVPrtlGen_module.cc:97

Member Data Documentation

bool evgen::ldm::MeVPrtlGen::fAnaOutput
private

Definition at line 91 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fDecayMaxWeight
private

Definition at line 104 of file MeVPrtlGen_module.cc.

std::unique_ptr<evgen::ldm::IMeVPrtlDecay> evgen::ldm::MeVPrtlGen::fDecayTool
private

Definition at line 99 of file MeVPrtlGen_module.cc.

bool evgen::ldm::MeVPrtlGen::fDoDeweight
private

Definition at line 93 of file MeVPrtlGen_module.cc.

CLHEP::HepJamesRandom* evgen::ldm::MeVPrtlGen::fEngine
private

Definition at line 110 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fFluxMaxWeight
private

Definition at line 102 of file MeVPrtlGen_module.cc.

std::unique_ptr<evgen::ldm::IMeVPrtlFlux> evgen::ldm::MeVPrtlGen::fFluxTool
private

Definition at line 97 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fGenMaxWeight
private

Definition at line 101 of file MeVPrtlGen_module.cc.

std::unique_ptr<evgen::ldm::IMesonGen> evgen::ldm::MeVPrtlGen::fGenTool
private

Definition at line 96 of file MeVPrtlGen_module.cc.

MeVPrtlTruth* evgen::ldm::MeVPrtlGen::fMeVPrtl
private

Definition at line 108 of file MeVPrtlGen_module.cc.

std::array<uint64_t, 4> evgen::ldm::MeVPrtlGen::fNCalls
private

Definition at line 112 of file MeVPrtlGen_module.cc.

std::array<double, 4> evgen::ldm::MeVPrtlGen::fNTime
private

Definition at line 113 of file MeVPrtlGen_module.cc.

bool evgen::ldm::MeVPrtlGen::fProduce
private

Definition at line 90 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fRayDecayMaxWeight
private

Definition at line 105 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fRayMaxWeight
private

Definition at line 103 of file MeVPrtlGen_module.cc.

std::unique_ptr<evgen::ldm::IRayTrace> evgen::ldm::MeVPrtlGen::fRayTool
private

Definition at line 98 of file MeVPrtlGen_module.cc.

double evgen::ldm::MeVPrtlGen::fSubRunPOT
private

Definition at line 94 of file MeVPrtlGen_module.cc.

TTree* evgen::ldm::MeVPrtlGen::fTree
private

Definition at line 107 of file MeVPrtlGen_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5