SnovaGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief An ART module for generating supernova events via GENIE
3 /// \author Andrey Sheshukov <_ash@list.ru>
4 ////////////////////////////////////////////////////////////////////////
5 
6 // C++ includes
7 #include <iostream>
8 #include <string>
9 #include <vector>
10 
11 // Framework includes
18 
19 #pragma GCC diagnostic push
20 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
21 // ROOT includes
22 #include "TVector3.h"
23 
24 // GENIE includes
25 #include "GENIE/Framework/Conventions/Constants.h"
26 #include "GENIE/Framework/EventGen/EventRecord.h"
27 #include "GENIE/Framework/EventGen/GMCJDriver.h"
28 #include "GENIE/Tools/Geometry/ROOTGeomAnalyzer.h"
29 #include "GENIE/Framework/GHEP/GHepUtils.h"
30 #include "GENIE/Framework/GHEP/GHepParticle.h"
31 #include "GENIE/Framework/Utils/AppInit.h"
32 #pragma GCC diagnostic pop
33 // NOvASoft includes
34 #include "Geometry/Geometry.h"
38 #include "SummaryData/RunData.h"
39 
40 
41 namespace supernova
42 {
43  class SnovaGen : public art::EDProducer
44  {
45  public:
46  explicit SnovaGen(const fhicl::ParameterSet& p);
47  ~SnovaGen();
48 
49  void produce(art::Event& evt);
50  void beginJob();
51  void beginRun(art::Run& run);
52 
53  double RecordTime() { return fRecord->Vertex()->T(); }
54  double EventStart() { return fTime; }
55  double EventEnd() { return fTime + fEventDuration; }
56  bool FluxEnd() { return fFluxDriver->End(); }
58 
60 
61  void PrepareEvent(art::Event& evt);
62  void SaveEvent(art::Event& evt);
63  void GenerateNextRecord();
64  void SaveRecordToEvent();
65 
66 
67  protected:
72 
75 
76  double DistanceKpc;
77  double Emin;
78 
79  TVector3 WindowPos;
80  double WindowSize;
81  double WindowX;
82  double WindowY;
83  double WindowZ;
84 
85  double fTime;
87 
89  std::unique_ptr<genie::geometry::ROOTGeomAnalyzer> fGeomDriver = nullptr;
90  std::unique_ptr<genie::supernova::GFluxSNova> fFluxDriver = nullptr;
91  std::unique_ptr<std::vector<simb::MCTruth>> fTruth = nullptr;
92 
93  TGeoManager* fGeoManager = nullptr;
94  TGeoVolume* fTopVolume = nullptr;
95 
96  size_t fRecordsCounter = 0;
97 
98  std::shared_ptr< genie::EventRecord > fRecord = nullptr; //< Generated interaction
99  };
100 }
101 
102 
103 ////////////////////////////////////////////////////////////////////////
104 namespace supernova
105 {
106  //.......................................................................
108  EDProducer(p),
109  MsgLevelFileName(p.get<std::string> ("MsgLevelFile")),
110  ModelFileName(p.get<std::string> ("ModelFile")),
111  MaxPathLenFileName(p.get<std::string> ("MaxPathLenFile")),
112  XsecFileName(p.get<std::string> ("XsecFile")),
113  TopVolume(p.get<std::string> ("TopVolumeName")),
114  EventGenList(p.get<std::string> ("EventGenList")),
115  DistanceKpc(p.get<double> ("DistanceKpc")),
116  Emin(p.get<double> ("Emin")),
117  WindowSize(p.get<double> ("Window.Size")),
118  WindowX(p.get<double> ("Window.X")),
119  WindowY(p.get<double> ("Window.Y")),
120  WindowZ(p.get<double> ("Window.Z")),
121  fTime(p.get<double> ("TimeStart")),
122  fEventDuration(p.get<double> ("EventDuration"))
123  {
124  cet::search_path sp0("GXMLPATH");
126 
127  cet::search_path sp("GSNOVA_DATA");
131 
132  std::cout << "Using files: "
133  << "\n ModelFile: \"" << ModelFile << "\""
134  << "\n MaxPathLenFile: \"" << MaxPathLenFile << "\""
135  << "\n XsecFile: \"" << XsecFile << "\""
136  << "\n MsgLevelFile: \"" << MsgLevelFile << "\""
137  << std::endl;
138 
139  WindowPos.SetX(WindowX);
140  WindowPos.SetY(WindowY);
141  WindowPos.SetZ(WindowZ);
142 
143  // Set messages level
145 
146  // setup geometry driver
148  fGeoManager = geo->ROOTGeoManager();
149  fTopVolume = fGeoManager->FindVolumeFast(TopVolume.c_str());
150  //fGeoManager->SetTopVolume(fTopVolume);
151 
153  fGeomDriver->SetTopVolName(TopVolume);
154  fGeomDriver->SetLengthUnits(genie::units::centimeter);
156  fGeomDriver->SetMixtureWeightsSum(1.);
157 
158  // Make model
159  auto snova_model = genie::supernova::GSNovaModel(ModelFile);
161 
162  snova_model.SetEmin(Emin);
163 
164  // Get flux driver
165  fFluxDriver.reset(new genie::supernova::GFluxSNova(snova_model, flux_window));
166 
167  // Read Xsec table
169 
170  // Create the GENIE monte carlo job driver
176 
177  // Charge the supernova!
179 
180  produces<std::vector<simb::MCTruth>> ();
181  produces<sumdata::RunData, art::InRun> ();
182  }
183 
184 
185  //......................................................................
187  {
188  }
189 
190 
191  //......................................................................
193  {
194  std::cout << "Begin job!" << std::endl;
195  }
196 
197 
198  //......................................................................
200  {
201  std::cout << "--- SupernovaGen::begin run\n";
202 
204  std::cout << "\tDetector ID: " << geo->DetId() << std::endl;
205 
206  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
207  geo->FileBaseName(),
208  geo->ExtractGDML()));
209  run.put(std::move(runcol));
210  return;
211  }
212 
213  //......................................................................
215  {
216  if(FluxEnd())return;
217  fRecord=std::unique_ptr<genie::EventRecord>(fMCJobDriver.GenerateEvent());
218  ++fRecordsCounter;
219  }
220 
221  //......................................................................
223  {
224  fTruth->push_back(MakeMCTruth(*fRecord));
225  }
226 
227  //......................................................................
229  {
230  fGeoManager->SetTopVolume(fTopVolume);
231  std::cout<<"Event Time=["<<EventStart()<<" - "<<EventEnd()<<"]"<<std::endl;
232  std::cout<<"~ starting gen ~"<<std::endl;
233  fTruth=std::make_unique<std::vector<simb::MCTruth>>();
234  }
235 
236  //......................................................................
238  {
239  std::cout<<"Event #"<<evt.id()<<" has "<<fTruth->size()<<" interactions"<<std::endl;
240  evt.put(std::move(fTruth));
241  fGeoManager->SetTopVolume(fGeoManager->GetMasterVolume());
242  EventTimeNext();
243  }
244 
245 
246  //......................................................................
248  {
249  ///This function was copied from GENIEHelper with minor modifications
250  simb::MCTruth truth;
251  TLorentzVector *vertex = record.Vertex();
252 
253  std::cout << "vtx=";
254  vertex->Print();
255  std::cout << std::endl;
256 
257  // Get the Interaction object from the record - this is the object
258  // that talks to the event information objects and is in m
259  genie::Interaction *inter = record.Summary();
260 
261  // Get the different components making up the interaction
262  const genie::InitialState &initState = inter->InitState();
263  const genie::ProcessInfo &procInfo = inter->ProcInfo();
264 
265  std::cout << "Neutrino E=" << initState.ProbeE(genie::kRfLab) << std::endl;
266 
267  // Add the particles from the interaction
268  TIter partitr(&record);
270 
271  // GHepParticles return units of GeV/c for p. The V_i are all in fermis
272  // and are relative to the center of the struck nucleus.
273  // Add the vertex X/Y/Z to the V_i for status codes 0 and 1
274  int trackid = 0;
275  std::string primary("primary");
276  while( (part = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ){
277 
278  simb::MCParticle tpart(trackid,
279  part->Pdg(),
280  primary,
281  part->FirstMother(),
282  part->Mass(),
283  part->Status());
284 
285  double vtx[4] = {part->Vx(), part->Vy(), part->Vz(), part->Vt()};
286  tpart.SetGvtx(vtx);
287  tpart.SetRescatter(part->RescatterCode());
288 
289  // Set the vertex location for the neutrino, nucleus and everything
290  // that is to be tracked. Vertex returns values in meters.
291  if(part->Status() == 0 || part->Status() == 1){
292  vtx[0] = 100*(part->Vx() * 1.e-15 + vertex->X());
293  vtx[1] = 100*(part->Vy() * 1.e-15 + vertex->Y());
294  vtx[2] = 100*(part->Vz() * 1.e-15 + vertex->Z());
295  //convert time to nanoseconds since event start
296  vtx[3] = part->Vt()+(vertex->T()-EventStart())*1e9;
297  }
298 
299  TLorentzVector pos(vtx[0], vtx[1], vtx[2], vtx[3]);
300  TLorentzVector mom(part->Px(), part->Py(), part->Pz(), part->E());
301 
302  tpart.AddTrajectoryPoint(pos, mom);
303 
304  if(part->PolzIsSet()) {
305  TVector3 polz;
306  part->GetPolarization(polz);
307  tpart.SetPolarization(polz);
308  }
309 
310  truth.Add(tpart);
311 
312  ++trackid;
313  } // End loop to convert GHepParticles to MCParticles
314 
315  // Get interaction category, mode, and type
316  int CCNC = procInfo.IsWeakNC() ? simb::kNC : simb::kCC;
318  int itype = simb::kUnknownInteraction;
319 
320  if (procInfo.IsNuElectronElastic() ) itype = simb::kNuElectronElastic;
321  else if (procInfo.IsInverseMuDecay() ) itype = simb::kInverseMuDecay;
322  else if (procInfo.IsQuasiElastic() ) itype = simb::kQE;
323  else if (procInfo.IsResonant() ) itype = simb::kRes;
324  else if (procInfo.IsDeepInelastic() ) itype = simb::kDIS;
325  else if (procInfo.IsCoherent() ) itype = simb::kCoh;
326  else if (procInfo.IsCoherentElas() ) itype = simb::kCohElastic;
327  else if (procInfo.IsElectronScattering() ) itype = simb::kElectronScattering;
328  else if (procInfo.IsIMDAnnihilation() ) itype = simb::kIMDAnnihilation;
329  else if (procInfo.IsInverseBetaDecay() ) itype = simb::kInverseBetaDecay;
330  else if (procInfo.IsGlashowResonance() ) itype = simb::kGlashowResonance;
331  else if (procInfo.IsAMNuGamma() ) itype = simb::kAMNuGamma;
332  else if (procInfo.IsMEC() ) itype = simb::kMEC;
333  else if (procInfo.IsDiffractive() ) itype = simb::kDiffractive;
334  else if (procInfo.IsEM() ) itype = simb::kEM;
335  else if (procInfo.IsWeakMix() ) itype = simb::kWeakMix;
336 
337  else itype = simb::kNuanceOffset + mode;
338 
339  // Set the neutrino information in MCTruth
341 
342  // The GENIE event kinematics are subtle different from the event
343  // kinematics that a experimentalist would calculate
344  // Instead of retriving the GENIE values for these kinematic variables
345  // calcuate them from the the final state particles
346  // while ingnoring the fermi momentum and the off-shellness of the bound nucleon.
347  genie::GHepParticle * hitnucl = record.HitNucleon();
348  TLorentzVector pdummy(0, 0, 0, 0);
349 
350  const TLorentzVector & k1 = *((record.Probe())->P4());
351  const TLorentzVector & k2 = *((record.FinalStatePrimaryLepton())->P4());
352  //const TLorentzVector & p1 = (hitnucl) ? *(hitnucl->P4()) : pdummy;
353 
355  TLorentzVector q = k1 - k2; // q = k1-k2, 4-p transfer
356 
357  double Q2 = -1 * q.M2(); // momemtum transfer
358  double v = (hitnucl) ? q.Energy() : -1; // v (E transfer to the nucleus)
359  double x = (hitnucl) ? 0.5 * Q2 / (M * v) : -1; // Bjorken x
360  double y = (hitnucl) ? v / k1.Energy() : -1; // Inelasticity, y = q * P1 / k1 * P1
361  double W2 = (hitnucl) ? M * M + 2 * M * v - Q2 : -1; // Hadronic invariant mass ^ 2
362  double W = (hitnucl) ? std::sqrt(W2) : -1; // Hardronic invariant mass
363 
364  truth.SetNeutrino(CCNC,
365  mode,
366  itype,
367  initState.Tgt().Pdg(),
368  initState.Tgt().HitNucPdg(),
369  initState.Tgt().HitQrkPdg(),
370  W,
371  x,
372  y,
373  Q2);
374 
375  return truth;
376  }
377 
378 
379  //.......................................................................
381  {
382  PrepareEvent(evt);
383  if(fRecord==nullptr)GenerateNextRecord();
384 
385  while(!FluxEnd() && RecordTime() < EventEnd()){
388  if(!FluxEnd())
389  std::cout<<"\rgenerated event #" << fRecordsCounter << " at time " << RecordTime() << std::endl;
390  }
391  SaveEvent(evt);
392  }
393 
394 } // End namespace supernova
395 
396 
397 ////////////////////////////////////////////////////////////////////////
398 namespace supernova
399 {
401 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:92
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
int RescatterCode(void) const
Definition: GHepParticle.h:66
std::shared_ptr< genie::EventRecord > fRecord
bool IsWeakMix(void) const
Definition: event.h:34
genie::GMCJDriver fMCJobDriver
std::string ModelFileName
void PrepareEvent(art::Event &evt)
std::string MsgLevelFile
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
neutrino electron elastic scatter
Definition: MCNeutrino.h:140
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
int HitNucPdg(void) const
Definition: Target.cxx:321
TGeoVolume * fTopVolume
int HitQrkPdg(void) const
Definition: Target.cxx:259
bool IsInverseMuDecay(void) const
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
const char * p
Definition: xmltok.h:285
int NuanceReactionCode(const GHepRecord *evrec)
Definition: GHepUtils.cxx:281
T sqrt(T number)
Definition: d0nt_math.hpp:156
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
bool IsElectronScattering(void) const
int Pdg(void) const
Definition: Target.h:72
void beginRun(art::Run &run)
TGeoManager * ROOTGeoManager() const
bool IsInverseBetaDecay(void) const
std::string MaxPathLenFile
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:95
bool IsDiffractive(void) const
double Mass(void) const
Mass that corresponds to the PDG code.
bool IsIMDAnnihilation(void) const
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
std::string EventGenList
std::string MsgLevelFileName
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
std::unique_ptr< genie::supernova::GFluxSNova > fFluxDriver
std::string find_file(std::string const &filename) const
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
Definition: Run.h:21
virtual TLorentzVector * Vertex(void) const
Definition: GHepRecord.h:141
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double Vt(void) const
Get production time.
Definition: GHepParticle.h:98
int Pdg(void) const
Definition: GHepParticle.h:64
SnovaGen(const fhicl::ParameterSet &p)
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
int FirstMother(void) const
Definition: GHepParticle.h:67
Summary information for an interaction.
Definition: Interaction.h:56
TString part[npart]
Definition: Style.C:32
A GENIE `MC Job Driver&#39;. Can be used for setting up complicated event generation cases involving deta...
Definition: GMCJDriver.h:47
double GlobProbScale(void) const
Definition: GMCJDriver.h:72
bool IsWeakNC(void) const
bool IsNuElectronElastic(void) const
std::unique_ptr< genie::geometry::ROOTGeomAnalyzer > fGeomDriver
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
bool IsAMNuGamma(void) const
std::string XsecFileName
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
static const double gram_centimeter3
Definition: Units.h:151
int evt
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:370
void SaveEvent(art::Event &evt)
std::string MaxPathLenFileName
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
void produce(art::Event &evt)
A ROOT/GEANT4 geometry driver.
TGeoManager * fGeoManager
bool PolzIsSet(void) const
Definition: run.py:1
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
OStream cout
Definition: OStream.cxx:6
bool IsMEC(void) const
bool IsEM(void) const
std::unique_ptr< std::vector< simb::MCTruth > > fTruth
void GetPolarization(TVector3 &polz)
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:82
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:80
double Vz(void) const
Get production z.
Definition: GHepParticle.h:97
inverse muon decay
Definition: MCNeutrino.h:141
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
An ART module for generating supernova events via GENIE.
Supernova neutrinos.
Definition: MCTruth.h:25
simb::MCTruth MakeMCTruth(const genie::EventRecord &record)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
static const double centimeter
Definition: Units.h:52
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
double Vy(void) const
Get production y.
Definition: GHepParticle.h:96
const Target & Tgt(void) const
Definition: InitialState.h:67
Event generator information.
Definition: MCTruth.h:32
bool IsCoherentElas(void) const
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
bool IsGlashowResonance(void) const
Helper for AttenCurve.
Definition: Path.h:10
#define W(x)
double ProbeE(RefFrame_t rf) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:179
A squre window, where neutrinos appear towards the origin (0,0,0)
Definition: GFluxWindow.h:12
EventID id() const
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string FileBaseName() const
double Vx(void) const
Get production x.
Definition: GHepParticle.h:95
Initial State information.
Definition: InitialState.h:49
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
enum BeamMode string