LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TextFileGen_module.cc
Go to the documentation of this file.
1 
58 #include <fstream>
59 #include <memory>
60 #include <string>
61 #include <utility>
62 
68 #include "cetlib_except/exception.h"
69 #include "fhiclcpp/ParameterSet.h"
71 
72 #include "TLorentzVector.h"
73 
78 
79 namespace evgen {
80  class TextFileGen;
81 }
82 
84 public:
85  explicit TextFileGen(fhicl::ParameterSet const& p);
86 
87  void produce(art::Event& e) override;
88  void beginJob() override;
89  void beginRun(art::Run& run) override;
90 
91 private:
92  std::pair<unsigned, unsigned> readEventInfo(std::istream& is);
94  unsigned long int fOffset;
95  std::ifstream* fInputFile;
96  std::string fInputFileName;
97  double fMoveY;
98 };
99 
100 //------------------------------------------------------------------------------
102  : EDProducer{p}
103  , fOffset{p.get<unsigned long int>("Offset")}
104  , fInputFile(0)
105  , fInputFileName{p.get<std::string>("InputFileName")}
106  , fMoveY{p.get<double>("MoveY", -1e9)}
107 {
108  if (fMoveY > -1e8) {
109  mf::LogWarning("TextFileGen") << "Particles will be moved to a new plane y = " << fMoveY
110  << " cm.\n";
111  }
112 
113  produces<std::vector<simb::MCTruth>>();
114  produces<sumdata::RunData, art::InRun>();
115 }
116 
117 //------------------------------------------------------------------------------
119 {
120  fInputFile = new std::ifstream(fInputFileName.c_str());
121 
122  // check that the file is a good one
123  if (!fInputFile->good())
124  throw cet::exception("TextFileGen")
125  << "input text file " << fInputFileName << " cannot be read.\n";
126 
127  for (unsigned i = 0; i != fOffset; ++i) {
128  auto const [eventNo, nparticles] = readEventInfo(*fInputFile);
129  for (unsigned p = 0; p != nparticles; ++p) {
130  constexpr auto all_chars_until = std::numeric_limits<unsigned>::max();
131  fInputFile->ignore(all_chars_until, '\n');
132  }
133  }
134 }
135 
136 //------------------------------------------------------------------------------
138 {
140  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
141 }
142 
143 //------------------------------------------------------------------------------
145 {
146  // check that the file is still good
147  if (!fInputFile->good())
148  throw cet::exception("TextFileGen")
149  << "input text file " << fInputFileName << " cannot be read in produce().\n";
150 
151  //Now, read the Event to be used.
152 
153  // check that the file is still good
154  if (!fInputFile->good())
155  throw cet::exception("TextFileGen")
156  << "input text file " << fInputFileName << " cannot be read in produce().\n";
157  auto truthcol = std::make_unique<std::vector<simb::MCTruth>>();
158  truthcol->push_back(readNextHepEvt());
159 
160  e.put(std::move(truthcol));
161 }
162 
164 {
165 
166  // declare the variables for reading in the event record
167  int status = 0;
168  int pdg = 0;
169  int firstMother = 0;
170  int secondMother = 0;
171  int firstDaughter = 0;
172  int secondDaughter = 0;
173  double xMomentum = 0.;
174  double yMomentum = 0.;
175  double zMomentum = 0.;
176  double energy = 0.;
177  double mass = 0.;
178  double xPosition = 0.;
179  double yPosition = 0.;
180  double zPosition = 0.;
181  double time = 0.;
182 
183  // read in line to get event number and number of particles
184  std::string oneLine;
185  std::istringstream inputLine;
186  simb::MCTruth nextEvent;
187  auto const [eventNo, nParticles] = readEventInfo(*fInputFile);
188 
189  // now read in all the lines for the particles
190  // in this interaction. only particles with
191  // status = 1 get tracked in Geant4.
192  for (unsigned short i = 0; i < nParticles; ++i) {
193  std::getline(*fInputFile, oneLine);
194  inputLine.clear();
195  inputLine.str(oneLine);
196 
197  inputLine >> status >> pdg >> firstMother >> secondMother >> firstDaughter >> secondDaughter >>
198  xMomentum >> yMomentum >> zMomentum >> energy >> mass >> xPosition >> yPosition >>
199  zPosition >> time;
200 
201  //Project the particle to a new y plane
202  if (fMoveY > -1e8) {
203  double totmom = sqrt(pow(xMomentum, 2) + pow(yMomentum, 2) + pow(zMomentum, 2));
204  double kx = xMomentum / totmom;
205  double ky = yMomentum / totmom;
206  double kz = zMomentum / totmom;
207  if (ky) {
208  double l = (fMoveY - yPosition) / ky;
209  xPosition += kx * l;
210  yPosition += ky * l;
211  zPosition += kz * l;
212  }
213  }
214 
215  TLorentzVector pos(xPosition, yPosition, zPosition, time);
216  TLorentzVector mom(xMomentum, yMomentum, zMomentum, energy);
217 
218  simb::MCParticle part(i, pdg, "primary", firstMother, mass, status);
219  part.AddTrajectoryPoint(pos, mom);
220 
221  nextEvent.Add(part);
222 
223  } // end loop on particles.
224 
225  return nextEvent;
226 }
227 
228 std::pair<unsigned, unsigned> evgen::TextFileGen::readEventInfo(std::istream& iss)
229 {
230  std::string line;
231  getline(iss, line);
232  std::istringstream buffer{line};
233 
234  // Parse read line for the event number and particles per event
235  unsigned event, nparticles;
236  buffer >> event >> nparticles;
237  return {event, nparticles};
238 }
239 
std::pair< unsigned, unsigned > readEventInfo(std::istream &is)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:258
constexpr auto fullRun()
unsigned long int fOffset
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
void produce(art::Event &e) override
Particle class.
Definition: Run.h:37
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
std::string fInputFileName
Name of text file containing events to simulate.
std::ifstream * fInputFile
TString part[npart]
Definition: Style.C:32
void beginJob() override
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
double energy
Definition: plottest35.C:25
void beginRun(art::Run &run) override
TextFileGen(fhicl::ParameterSet const &p)
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:203
Event generator information.
Definition: MCTruth.h:32
simb::MCTruth readNextHepEvt()
Float_t e
Definition: plot.C:35
Namespace collecting geometry-related classes utilities.
Event Generation using GENIE, cosmics or single particles.
double fMoveY
Project particles to a new y plane.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.