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  MsgLevelFileName(p.get<std::string> ("MsgLevelFile")),
109  ModelFileName(p.get<std::string> ("ModelFile")),
110  MaxPathLenFileName(p.get<std::string> ("MaxPathLenFile")),
111  XsecFileName(p.get<std::string> ("XsecFile")),
112  TopVolume(p.get<std::string> ("TopVolumeName")),
113  EventGenList(p.get<std::string> ("EventGenList")),
114  DistanceKpc(p.get<double> ("DistanceKpc")),
115  Emin(p.get<double> ("Emin")),
116  WindowSize(p.get<double> ("Window.Size")),
117  WindowX(p.get<double> ("Window.X")),
118  WindowY(p.get<double> ("Window.Y")),
119  WindowZ(p.get<double> ("Window.Z")),
120  fTime(p.get<double> ("TimeStart")),
121  fEventDuration(p.get<double> ("EventDuration"))
122  {
123  cet::search_path sp0("GXMLPATH");
125 
126  cet::search_path sp("GSNOVA_DATA");
130 
131  std::cout << "Using files: "
132  << "\n ModelFile: \"" << ModelFile << "\""
133  << "\n MaxPathLenFile: \"" << MaxPathLenFile << "\""
134  << "\n XsecFile: \"" << XsecFile << "\""
135  << "\n MsgLevelFile: \"" << MsgLevelFile << "\""
136  << std::endl;
137 
138  WindowPos.SetX(WindowX);
139  WindowPos.SetY(WindowY);
140  WindowPos.SetZ(WindowZ);
141 
142  // Set messages level
144 
145  // setup geometry driver
147  fGeoManager = geo->ROOTGeoManager();
148  fTopVolume = fGeoManager->FindVolumeFast(TopVolume.c_str());
149  //fGeoManager->SetTopVolume(fTopVolume);
150 
152  fGeomDriver->SetTopVolName(TopVolume);
153  fGeomDriver->SetLengthUnits(genie::units::centimeter);
155  fGeomDriver->SetMixtureWeightsSum(1.);
156 
157  // Make model
158  auto snova_model = genie::supernova::GSNovaModel(ModelFile);
160 
161  snova_model.SetEmin(Emin);
162 
163  // Get flux driver
164  fFluxDriver.reset(new genie::supernova::GFluxSNova(snova_model, flux_window));
165 
166  // Read Xsec table
168 
169  // Create the GENIE monte carlo job driver
175 
176  // Charge the supernova!
178 
179  produces<std::vector<simb::MCTruth>> ();
180  produces<sumdata::RunData, art::InRun> ();
181  }
182 
183 
184  //......................................................................
186  {
187  }
188 
189 
190  //......................................................................
192  {
193  std::cout << "Begin job!" << std::endl;
194  }
195 
196 
197  //......................................................................
199  {
200  std::cout << "--- SupernovaGen::begin run\n";
201 
203  std::cout << "\tDetector ID: " << geo->DetId() << std::endl;
204 
205  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
206  geo->FileBaseName(),
207  geo->ExtractGDML()));
208  run.put(std::move(runcol));
209  return;
210  }
211 
212  //......................................................................
214  {
215  if(FluxEnd())return;
216  fRecord=std::unique_ptr<genie::EventRecord>(fMCJobDriver.GenerateEvent());
217  ++fRecordsCounter;
218  }
219 
220  //......................................................................
222  {
223  fTruth->push_back(MakeMCTruth(*fRecord));
224  }
225 
226  //......................................................................
228  {
229  fGeoManager->SetTopVolume(fTopVolume);
230  std::cout<<"Event Time=["<<EventStart()<<" - "<<EventEnd()<<"]"<<std::endl;
231  std::cout<<"~ starting gen ~"<<std::endl;
232  fTruth=std::make_unique<std::vector<simb::MCTruth>>();
233  }
234 
235  //......................................................................
237  {
238  std::cout<<"Event #"<<evt.id()<<" has "<<fTruth->size()<<" interactions"<<std::endl;
239  evt.put(std::move(fTruth));
240  fGeoManager->SetTopVolume(fGeoManager->GetMasterVolume());
241  EventTimeNext();
242  }
243 
244 
245  //......................................................................
247  {
248  ///This function was copied from GENIEHelper with minor modifications
249  simb::MCTruth truth;
250  TLorentzVector *vertex = record.Vertex();
251 
252  std::cout << "vtx=";
253  vertex->Print();
254  std::cout << std::endl;
255 
256  // Get the Interaction object from the record - this is the object
257  // that talks to the event information objects and is in m
258  genie::Interaction *inter = record.Summary();
259 
260  // Get the different components making up the interaction
261  const genie::InitialState &initState = inter->InitState();
262  const genie::ProcessInfo &procInfo = inter->ProcInfo();
263 
264  std::cout << "Neutrino E=" << initState.ProbeE(genie::kRfLab) << std::endl;
265 
266  // Add the particles from the interaction
267  TIter partitr(&record);
269 
270  // GHepParticles return units of GeV/c for p. The V_i are all in fermis
271  // and are relative to the center of the struck nucleus.
272  // Add the vertex X/Y/Z to the V_i for status codes 0 and 1
273  int trackid = 0;
274  std::string primary("primary");
275  while( (part = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ){
276 
277  simb::MCParticle tpart(trackid,
278  part->Pdg(),
279  primary,
280  part->FirstMother(),
281  part->Mass(),
282  part->Status());
283 
284  double vtx[4] = {part->Vx(), part->Vy(), part->Vz(), part->Vt()};
285  tpart.SetGvtx(vtx);
286  tpart.SetRescatter(part->RescatterCode());
287 
288  // Set the vertex location for the neutrino, nucleus and everything
289  // that is to be tracked. Vertex returns values in meters.
290  if(part->Status() == 0 || part->Status() == 1){
291  vtx[0] = 100*(part->Vx() * 1.e-15 + vertex->X());
292  vtx[1] = 100*(part->Vy() * 1.e-15 + vertex->Y());
293  vtx[2] = 100*(part->Vz() * 1.e-15 + vertex->Z());
294  //convert time to nanoseconds since event start
295  vtx[3] = part->Vt()+(vertex->T()-EventStart())*1e9;
296  }
297 
298  TLorentzVector pos(vtx[0], vtx[1], vtx[2], vtx[3]);
299  TLorentzVector mom(part->Px(), part->Py(), part->Pz(), part->E());
300 
301  tpart.AddTrajectoryPoint(pos, mom);
302 
303  if(part->PolzIsSet()) {
304  TVector3 polz;
305  part->GetPolarization(polz);
306  tpart.SetPolarization(polz);
307  }
308 
309  truth.Add(tpart);
310 
311  ++trackid;
312  } // End loop to convert GHepParticles to MCParticles
313 
314  // Get interaction category, mode, and type
315  int CCNC = procInfo.IsWeakNC() ? simb::kNC : simb::kCC;
317  int itype = simb::kUnknownInteraction;
318 
319  if (procInfo.IsNuElectronElastic() ) itype = simb::kNuElectronElastic;
320  else if (procInfo.IsInverseMuDecay() ) itype = simb::kInverseMuDecay;
321  else if (procInfo.IsQuasiElastic() ) itype = simb::kQE;
322  else if (procInfo.IsResonant() ) itype = simb::kRes;
323  else if (procInfo.IsDeepInelastic() ) itype = simb::kDIS;
324  else if (procInfo.IsCoherent() ) itype = simb::kCoh;
325  else if (procInfo.IsCoherentElas() ) itype = simb::kCohElastic;
326  else if (procInfo.IsElectronScattering() ) itype = simb::kElectronScattering;
327  else if (procInfo.IsIMDAnnihilation() ) itype = simb::kIMDAnnihilation;
328  else if (procInfo.IsInverseBetaDecay() ) itype = simb::kInverseBetaDecay;
329  else if (procInfo.IsGlashowResonance() ) itype = simb::kGlashowResonance;
330  else if (procInfo.IsAMNuGamma() ) itype = simb::kAMNuGamma;
331  else if (procInfo.IsMEC() ) itype = simb::kMEC;
332  else if (procInfo.IsDiffractive() ) itype = simb::kDiffractive;
333  else if (procInfo.IsEM() ) itype = simb::kEM;
334  else if (procInfo.IsWeakMix() ) itype = simb::kWeakMix;
335 
336  else itype = simb::kNuanceOffset + mode;
337 
338  // Set the neutrino information in MCTruth
340 
341  // The GENIE event kinematics are subtle different from the event
342  // kinematics that a experimentalist would calculate
343  // Instead of retriving the GENIE values for these kinematic variables
344  // calcuate them from the the final state particles
345  // while ingnoring the fermi momentum and the off-shellness of the bound nucleon.
346  genie::GHepParticle * hitnucl = record.HitNucleon();
347  TLorentzVector pdummy(0, 0, 0, 0);
348 
349  const TLorentzVector & k1 = *((record.Probe())->P4());
350  const TLorentzVector & k2 = *((record.FinalStatePrimaryLepton())->P4());
351  //const TLorentzVector & p1 = (hitnucl) ? *(hitnucl->P4()) : pdummy;
352 
354  TLorentzVector q = k1 - k2; // q = k1-k2, 4-p transfer
355 
356  double Q2 = -1 * q.M2(); // momemtum transfer
357  double v = (hitnucl) ? q.Energy() : -1; // v (E transfer to the nucleus)
358  double x = (hitnucl) ? 0.5 * Q2 / (M * v) : -1; // Bjorken x
359  double y = (hitnucl) ? v / k1.Energy() : -1; // Inelasticity, y = q * P1 / k1 * P1
360  double W2 = (hitnucl) ? M * M + 2 * M * v - Q2 : -1; // Hadronic invariant mass ^ 2
361  double W = (hitnucl) ? std::sqrt(W2) : -1; // Hardronic invariant mass
362 
363  truth.SetNeutrino(CCNC,
364  mode,
365  itype,
366  initState.Tgt().Pdg(),
367  initState.Tgt().HitNucPdg(),
368  initState.Tgt().HitQrkPdg(),
369  W,
370  x,
371  y,
372  Q2);
373 
374  return truth;
375  }
376 
377 
378  //.......................................................................
380  {
381  PrepareEvent(evt);
382  if(fRecord==nullptr)GenerateNextRecord();
383 
384  while(!FluxEnd() && RecordTime() < EventEnd()){
387  if(!FluxEnd())
388  std::cout<<"\rgenerated event #" << fRecordsCounter << " at time " << RecordTime() << std::endl;
389  }
390  SaveEvent(evt);
391  }
392 
393 } // End namespace supernova
394 
395 
396 ////////////////////////////////////////////////////////////////////////
397 namespace supernova
398 {
400 }
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:80
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
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
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
std::string MsgLevelFileName
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
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:31
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
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
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
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
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
Definition: Event.h:56
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