LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
FileMuons_module.cc
Go to the documentation of this file.
1 
10 // C++ includes
11 #include <cmath>
12 #include <fstream>
13 #include <iostream>
14 #include <memory>
15 #include <string>
16 #include <utility>
17 #include <vector>
18 
19 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
28 
29 // nusimdata includes
32 
33 // LArSoft includes
36 
37 // ROOT includves
38 #include "TDatabasePDG.h"
39 #include "TFile.h"
40 #include "TTree.h"
41 #include "TVector3.h"
42 
43 namespace evgen {
44 
46  class FileMuons : public art::EDProducer {
47  public:
48  explicit FileMuons(fhicl::ParameterSet const& pset);
49 
50  private:
51  void produce(art::Event& evt) override;
52  void beginJob() override;
53  void beginRun(art::Run& run) override;
54  void endJob() override;
55 
56  void ReadEvents(simb::MCTruth& mct);
57 
58  int fEventNumberOffset; // Where in file to start.
59  std::vector<int> fPDG;
60  std::vector<double> fXYZ_Off;
61  std::string fFileName;
62  std::string fMuonsFileType;
63  std::string fTreeName;
64  std::vector<std::string> fBranchNames;
65 
66  std::ifstream* fMuonFile;
67  TFile* fMuonFileR;
68  TTree* TNtuple;
69  unsigned int countFile;
70 
71  Float_t xtmp, ytmp, ztmp;
72  Float_t pxtmp, pytmp, pztmp;
73  Float_t charge;
74  Float_t E;
75  Float_t costheta;
76  Float_t phi;
77  Float_t xdet;
78  Float_t ydet;
79  Float_t zdet;
80 
81  TBranch* b_x;
82  TBranch* b_y;
83  TBranch* b_z;
84  TBranch* b_E;
85  TBranch* b_costheta;
86  TBranch* b_phi;
87  TBranch* b_xdet;
88  TBranch* b_ydet;
89  TBranch* b_zdet;
90  TBranch* b_px;
91  TBranch* b_py;
92  TBranch* b_pz;
93  TBranch* b_charge;
94  };
95 } // namespace
96 
97 namespace evgen {
98 
99  //____________________________________________________________________________
101  : EDProducer{pset}
102  , fEventNumberOffset(pset.get<int>("EventNumberOffset"))
103  , fPDG(pset.get<std::vector<int>>("PDG"))
104  , fXYZ_Off(pset.get<std::vector<double>>("InitialXYZOffsets"))
105  , fFileName(pset.get<std::string>("FileName"))
106  , fMuonsFileType(pset.get<std::string>("MuonsFileType"))
107  , fTreeName(pset.get<std::string>("TreeName"))
108  , fBranchNames(pset.get<std::vector<std::string>>("BranchNames"))
109  {
110 
111  produces<std::vector<simb::MCTruth>>();
112  produces<sumdata::RunData, art::InRun>();
113  }
114 
115  //____________________________________________________________________________
117  {
119  mf::LogInfo("FileMuons : starting at event ") << countFile << std::endl;
120 
121  if (fMuonsFileType.compare("source") == 0) {
122  std::cout << "FileMuons: Not yet equipped to walk through muons with TFS mojo." << std::endl;
123  }
124  else if (fMuonsFileType.compare("root") == 0) {
125  std::cout << "FileMuons: You have chosen to read muons from Root File " << fFileName
126  << std::endl;
127  }
128  else if (fMuonsFileType.compare("text") == 0) {
129  std::cout << "FileMuons: You have chosen to read muons from " << fFileName << "."
130  << std::endl;
131  }
132  else {
133  std::cout << "FileMuons: You must specify one of source/text/root file to read for muons."
134  << std::endl;
135  }
136 
137  if (fMuonsFileType.compare("text") == 0) {
138  fMuonFile = new std::ifstream(fFileName.c_str());
139  long begin = fMuonFile->tellg();
140  fMuonFile->seekg(0, std::ios::end);
141  long end = fMuonFile->tellg();
142  std::cout << "FileMuons: " << fFileName << " size is: " << (end - begin) << " bytes.\n";
143  fMuonFile->seekg(0, std::ios::beg);
144 
145  for (unsigned int header = 0; header < 3 && fMuonFile->good(); ++header) {
146  std::string line;
147  getline(*fMuonFile, line);
148  }
149  if (!fMuonFile->good()) {
150  std::cout << "FileMuons: Problem reading muon file header." << std::endl;
151  }
152  } // fMuonsFileType is a text file.
153  else if (fMuonsFileType.compare("root") == 0) {
154  fMuonFileR = new TFile(fFileName.c_str(), "READ");
155  TNtuple = (TTree*)(fMuonFileR->Get(fTreeName.c_str()));
156 
157  TNtuple->SetBranchAddress("x", &xtmp, &b_x);
158  TNtuple->SetBranchAddress("y", &ytmp, &b_y);
159  TNtuple->SetBranchAddress("z", &ztmp, &b_z);
160  TNtuple->SetBranchAddress("E", &E, &b_E);
161  TNtuple->SetBranchAddress("costheta", &costheta, &b_costheta);
162  TNtuple->SetBranchAddress("phi", &phi, &b_phi);
163  TNtuple->SetBranchAddress("xdet", &xdet, &b_xdet);
164  TNtuple->SetBranchAddress("ydet", &ydet, &b_ydet);
165  TNtuple->SetBranchAddress("zdet", &zdet, &b_zdet);
166  TNtuple->SetBranchAddress("px", &pxtmp, &b_px);
167  TNtuple->SetBranchAddress("py", &pytmp, &b_py);
168  TNtuple->SetBranchAddress("pz", &pztmp, &b_pz);
169  TNtuple->SetBranchAddress("charge", &charge, &b_charge);
170 
171  } // fMuonsFileType is a root file.
172  }
173 
174  //____________________________________________________________________________
176  {
177  if (fMuonsFileType.compare("text") == 0) fMuonFile->close();
178  if (fMuonsFileType.compare("root") == 0) fMuonFileR->Close();
179  }
180 
181  //____________________________________________________________________________
183  {
185  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
186  }
187 
188  //____________________________________________________________________________
190  {
191 
193  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
194 
195  simb::MCTruth truth;
197  ReadEvents(truth);
198 
199  truthcol->push_back(truth);
200 
201  evt.put(std::move(truthcol));
202  }
203 
204  //____________________________________________________________________________
206  {
208  auto const& cryostat = art::ServiceHandle<geo::Geometry const>
209  {
210  } -> Cryostat();
211  for (unsigned int i = 0; i < fPDG.size(); ++i) {
212 
213  // Choose momentum
214  double m(0.108);
215 
216  TVector3 x;
217  TVector3 p;
218  Double_t q = 0.;
219  Int_t pdgLocal;
220 
221  if (fMuonsFileType.compare("text") == 0) {
222 
223  std::string line;
224  getline(*fMuonFile, line);
225  if (!fMuonFile->good()) {
227  << "FileMuons: Problem reading muon file line ...." << countFile
228  << ". Perhaps you've exhausted the events in " << fFileName << std::endl;
229  }
230  countFile++;
231 
232  MF_LOG_DEBUG("FileMuons: countFile is ") << countFile << std::endl;
233  char *cstr, *ptok;
234 
235  // Split this line into tokens
236  cstr = new char[line.size() + 1];
237  strcpy(cstr, line.c_str());
238  // cstr now contains a c-string copy of str
239  ptok = strtok(cstr, "*");
240  unsigned int fieldCount = 0;
241  unsigned int posIndex = 0;
242  unsigned int pIndex = 0;
243  while (ptok != NULL) {
244  ptok = strtok(NULL, "*");
245  if (fieldCount == 9 || fieldCount == 10 || fieldCount == 11) {
246  p[pIndex] = atof(ptok);
247  pIndex++;
248  }
249  if (fieldCount == 6 || fieldCount == 7 || fieldCount == 8) {
250  x[posIndex] = atof(ptok);
251  // make the z axis point up for x, as with p
252  if (posIndex == 2) { x[posIndex] = -1.0 * x[posIndex]; }
253  posIndex++;
254  }
255  if (fieldCount == 12) { q = atof(ptok); }
256  fieldCount++;
257  }
258 
259  delete[] cstr;
260  }
261  else if (fMuonsFileType.compare("root") == 0) // from root file
262  {
263  TNtuple->GetEntry(countFile);
264 
265  x.SetXYZ(xdet, ydet, -zdet); // as with txt file, make z point up.
266  // Watch for units change to mm in Modern JdJ files!!
267  // This is for pre Spring-2012 JdJ Ntuples.
268  p.SetXYZ(pxtmp, pytmp, pztmp);
269  q = charge;
270 
271  countFile++;
272 
273  } // End read.
274 
275  static TDatabasePDG pdgt;
276  pdgLocal = -q * fPDG[i];
277 
278  TParticlePDG* pdgp = pdgt.GetParticle(pdgLocal);
279  if (pdgp) m = pdgp->Mass();
280 
281  // This gives coordinates at the center of the 300mx300m plate that is 3m above top of
282  // cavern. Got these by histogramming deJong's xdet,ydet,zdet.
283  const double cryoGap = 15.0;
284  x[0] -= fXYZ_Off[0];
285  x[1] -= fXYZ_Off[1];
286  x[2] -= fXYZ_Off[2]; // 3 for plate height above top of cryostat.
287  // Now, must rotate to TPC coordinates. Let's orient TPC axis along z axis,
288  // Cosmics, mostly going along deJong's +z axis must be going along TPC -y axis.
289  x.RotateX(-M_PI / 2);
290  p.RotateX(-M_PI / 2);
291  //add vector of the position of the center of the point between Cryostats
292  // level with top. (To which I've added 3m - in above code - in height.)
293  // This is referenced from origin at center-right of first cryostat.
294  TVector3 off3(cryostat.HalfWidth() * 0.01,
295  cryostat.HalfHeight() * 0.01,
296  cryostat.Length() * 0.01 + cryoGap * 0.01 / 2.0);
297  x += off3;
298 
299  TLorentzVector pos(x[0] * 100.0, x[1] * 100.0, x[2] * 100.0, 0.0);
300  TLorentzVector pvec(
301  p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0, std::sqrt(p.Mag2() * 1000.0 * 1000.0 + m * m));
302  std::cout << "x[m] and p [TeV] are " << std::endl;
303  x.Print();
304  p.Print();
305 
306  int trackid =
307  -1 * (i + 1); // set track id to -i as these are all primary particles and have id <= 0
308  std::string primary("primary");
309  simb::MCParticle part(trackid, pdgLocal, primary);
310  part.AddTrajectoryPoint(pos, pvec);
311  mct.Add(part);
312 
313  } //end loop over particles
314  }
315 
317 
318 } //end namespace evgen
Float_t x
Definition: compare.C:6
std::vector< double > fXYZ_Off
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:258
constexpr auto fullRun()
void beginRun(art::Run &run) override
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void beginJob() override
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
std::string fMuonsFileType
void ReadEvents(simb::MCTruth &mct)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
std::ifstream * fMuonFile
Particle class.
Definition: Run.h:37
std::vector< std::string > fBranchNames
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
TString part[npart]
Definition: Style.C:32
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
single particles thrown at the detector
Definition: MCTruth.h:26
unsigned int countFile
std::string fFileName
module to produce single or multiple specified particles in the detector
FileMuons(fhicl::ParameterSet const &pset)
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
#define MF_LOG_DEBUG(id)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:140
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
ROOT libraries.
void endJob() override
Event Generation using GENIE, cosmics or single particles.
std::vector< int > fPDG
art framework interface to geometry description
std::string fTreeName
void produce(art::Event &evt) override