IsoGen_module.cc
Go to the documentation of this file.
1 /// \author Zukai Wang - ricken@fnal.gov
2 /// \date Sept. 201
3 
4 // C++ includes.
5 #include <iostream>
6 #include <string>
7 #include <cmath>
8 #include <memory>
9 #include <vector>
10 
11 // ROOT includes.
12 #include "TDatabasePDG.h"
13 #include "TMath.h"
14 #include "TVector3.h"
15 
16 #include "Geometry/Geometry.h"
19 #include "nugen/EventGeneratorBase/evgenbase.h"
20 #include "SummaryData/RunData.h"
21 
22 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
32 
33 #include "CLHEP/Random/RandFlat.h"
35 
36 #include <fstream>
37 
38 #include "TTree.h"
39 
40 namespace evgen {
41  class SingleParticle;
42 
43  /// A module to check the results from the Monte Carlo generator
44  class IsoGen : public art::EDProducer {
45  public:
46  explicit IsoGen(fhicl::ParameterSet const& pset);
47  virtual ~IsoGen();
48 
49  void produce(art::Event& evt);
50  void reconfigure(const fhicl::ParameterSet& pset);
51  void beginRun(art::Run& run);
52 
53  private:
54 
55  void Sample(simb::MCTruth &truth);
56 
57  std::vector<int> fPDG; ///< PDG code of particles to generate
58 
59  std::vector<double> fP0; ///< Central momentum (GeV/c) to generate
60  std::vector<double> fSigmaP; ///< Variation in momenta (GeV/c)
61  std::vector<double> fPL; ///< momentum list (GeV/c) to generate
62  std::vector<int> fPDist; ///< How to distribute momenta (0=uniform, 1=gaussian)
63 
64  std::vector<double> fXhigh; ///< Up limit of x position (cm)
65  std::vector<double> fYhigh; ///< Up limit of y position (cm)
66  std::vector<double> fZhigh; ///< Up limit of z position (cm)
67  std::vector<double> fT0; ///< Central t position (ns)
68  std::vector<double> fXlow; ///< Low limit of x position (cm)
69  std::vector<double> fYlow; ///< Low Limit of y position (cm)
70  std::vector<double> fZlow; ///< Low Limit of z position (cm)
71  std::vector<double> fSigmaT; ///< Variation in t position (ns)
72 
73  TTree * isoGen_Tree;
76  void setInSurfID(CLHEP::RandFlat &flat, int const& i);
77  void setDirection(CLHEP::RandFlat &flat, int const& i);
78  bool getOutSurfPosition(int i);
79  bool projectToSurface(int axis, double surfaceLoc);
80  };
81 }
82 
83 
84 ////////////////////////////////////////////////////////////////////////
85 namespace evgen{
86 
87  //____________________________________________________________________________
89  {
90  this->reconfigure(pset);
91 
92  // get the random number seed, use a random default if not specified
93  // in the configuration file.
94  int seed = pset.get< unsigned int >("Seed", evgb::GetRandomNumberSeed());
95 
96  createEngine( seed );
97 
98  produces< std::vector<simb::MCTruth> >();
99  produces< sumdata::RunData, art::InRun >();
100  }
101 
102  //____________________________________________________________________________
104 
105  //____________________________________________________________________________
107  {
108  fXlow = pset.get< std::vector<double> >("Xlow");
109  fYlow = pset.get< std::vector<double> >("Ylow");
110  fZlow = pset.get< std::vector<double> >("Zlow");
111  fXhigh = pset.get< std::vector<double> >("Xhigh");
112  fYhigh = pset.get< std::vector<double> >("Yhigh");
113  fZhigh = pset.get< std::vector<double> >("Zhigh");
114  fPDG = pset.get< std::vector<int> >("PDG");
115  fP0 = pset.get< std::vector<double> >("P0");
116  fSigmaP = pset.get< std::vector<double> >("SigmaP");
117  fPL = pset.get< std::vector<double> >("PL");
118  fPDist = pset.get< std::vector<int> >("PDist");
119  fT0 = pset.get< std::vector<double> >("T0");
120  fSigmaT = pset.get< std::vector<double> >("SigmaT");
121  }
122 
123  //____________________________________________________________________________
125  {
126  // grab the geometry object to see what geometry we are using
128 
129  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
130  geo->FileBaseName(),
131  geo->ExtractGDML()));
132 
134  isoGen_Tree = tfs->make<TTree>("isoGen_Tree","Information from Generator");
135  isoGen_Tree -> Branch("inSurfID", &inSurfID);
136  isoGen_Tree -> Branch("outSurfID", &outSurfID);
137  isoGen_Tree -> Branch("cosx", &cosx);
138  isoGen_Tree -> Branch("cosy", &cosy);
139  isoGen_Tree -> Branch("cosz", &cosz);
140  isoGen_Tree -> Branch("ix", &init_x[0]);
141  isoGen_Tree -> Branch("iy", &init_x[1]);
142  isoGen_Tree -> Branch("iz", &init_x[2]);
143  isoGen_Tree -> Branch("it", &init_x[3]);
144  isoGen_Tree -> Branch("fx", &final_x[0]);
145  isoGen_Tree -> Branch("fy", &final_x[1]);
146  isoGen_Tree -> Branch("fz", &final_x[2]);
147  isoGen_Tree -> Branch("localPhi", &localPhi);
148  isoGen_Tree -> Branch("localTheta", &localTheta);
149 
150  run.put(std::move(runcol));
151 
152  return;
153  }
154 
155  //____________________________________________________________________________
157  {
158  ///unique_ptr allows ownership to be transferred to the art::Event after the ut statement
159  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
160 
161  simb::MCTruth truth;
163  Sample(truth);
164 
165  // std::cout << "put mctruth into the vector" << std::endl;
166  truthcol->push_back(truth);
167 
168  // std::cout << "add vector to the event " << truthcol->size() << std::endl;
169  evt.put(std::move(truthcol));
170 
171  return;
172  }
173 
174  //____________________________________________________________________________
176  {
177  // get the random number generator service and make some CLHEP generators
179 
180  CLHEP::HepRandomEngine &engine = rng->getEngine();
181  CLHEP::RandFlat flat(engine);
182  CLHEP::RandGaussQ gauss(engine);
183 
184  //std::ofstream ofs;
185  //ofs.open ("iso.txt", std::ofstream::out | std::ofstream::app);
186 
187  // every event will have one of each particle species in the fPDG array
188  for (unsigned int i=0; i<fPDG.size(); ++i) {
189  // Choose momentum
190  double p = 0.0;
191  double m = 0.0;
192 
193  if (fPL.size()==0){
194  if (fPDist[i]) {
195  p = gauss.fire(fP0[i], fSigmaP[i]);
196  }
197  else {
198  p = flat.fire(fP0[i]-fSigmaP[i], fP0[i]+fSigmaP[i]);
199  }
200  }
201  else {
202  int randomIndex = int(flat.fireInt(fPL.size()));
203  p = fPL[randomIndex];
204  }
205 
206  // ofs << "get the mass" << std::endl;
207  static TDatabasePDG* pdgt = TDatabasePDG::Instance();
208  TParticlePDG* pdgp = pdgt->GetParticle(fPDG[i]);
209  if (pdgp) m = pdgp->Mass();
210 
211  setInSurfID(flat, i);
212  // Choose position and time
213  init_x[3] = flat.fire(fT0[i]-fSigmaT[i], fT0[i]+fSigmaT[i]);
214  // Choose angles
215  setDirection(flat, i);
217 
218  isoGen_Tree -> Fill();
219 
220  TLorentzVector pvec(cosx*p,
221  cosy*p,
222  cosz*p,
223  sqrt(p*p+m*m));
224 
225  // ofs<<cosx<<" "<<cosy<<" "<<cosz<<std::endl;
226  TLorentzVector pos(init_x[0],
227  init_x[1],
228  init_x[2],
229  init_x[3]);
230 
231  // set track id to -i as these are all primary particles and have
232  // id <= 0
233  //
234  int trackid = -1*(i+1);
235 
236  std::string primary("primary");
237  simb::MCParticle part(trackid, fPDG[i], primary);
238  part.AddTrajectoryPoint(pos, pvec);
239 
240  mct.Add(part);
241  } //end loop over particles
242 
243  // ofs.close();
244  }
245 
247  //int inSurfID = 0;
248  //Choose an incident surface:
249  double lx,ly,lz;
250  lx = fXhigh[i]-fXlow[i];
251  ly = fYhigh[i]-fYlow[i];
252  lz = fZhigh[i]-fZlow[i];
253 
254  double Surf[6];
255  Surf[0] = lx*ly;
256  Surf[1] = lx*lz;
257  Surf[2] = ly*lz;
258  Surf[3] = lx*lz;
259  Surf[4] = ly*lz;
260  Surf[5] = lx*ly;
261 
262  double sumA =0;
263  for (int j=0; j!=6;++j) sumA += Surf[j];
264 
265  double dice = flat.fire(0,sumA);
266  if (dice<Surf[0]) {
267  inSurfID=0;
268  return;
269  }
270  //incident from end(5):
271  else if (dice<Surf[0]+Surf[5]) {
272  inSurfID=5;
273  return;
274  }
275  //incident from top(1):
276  else if (dice<Surf[0]+Surf[5]+Surf[1]) {
277  inSurfID=1;
278  return;
279  }
280  //incident from bottom(3):
281  else if (dice<sumA-Surf[4]-Surf[2]) {
282  inSurfID=3;
283  return;
284  }
285  //incident from left(4):
286  else if (dice<sumA-Surf[2]) {
287  inSurfID=4;
288  return;
289  }
290  //incident from right(2):
291  else {
292  inSurfID=2;
293  return;
294  }
295  //return inSurfID;
296  }
297 
299  //double localPhi;
300  localPhi = flat.fire(0,2.*TMath::Pi());
301  //double pass = flat.fire(0,1);
302  //incident from front(0):
303  if (inSurfID == 0) {
304  init_x[0] = flat.fire(fXlow[i],fXhigh[i]);
305  init_x[1] = flat.fire(fYlow[i],fYhigh[i]);
306  init_x[2] = fZhigh[i];
307 
308  cosz = -sqrt(flat.fire(0,1));
309  cosx = sqrt(1-cosz*cosz)*TMath::Cos(localPhi);
310  cosy = sqrt(1-cosz*cosz)*TMath::Sin(localPhi);
311  localTheta = TMath::ACos(-cosz);
312  return;
313  }
314 
315  //incident from end(5):
316  else if (inSurfID == 5) {
317  init_x[0] = flat.fire(fXlow[i],fXhigh[i]);
318  init_x[1] = flat.fire(fYlow[i],fYhigh[i]);
319  init_x[2] = fZlow[i];
320  cosz = sqrt(flat.fire(0,1));
321  cosx = sqrt(1-cosz*cosz)*TMath::Cos(localPhi);
322  cosy = sqrt(1-cosz*cosz)*TMath::Sin(localPhi);
323  localTheta = TMath::ACos(cosz);
324  return;
325  }
326 
327  //incident from top(1):
328  else if (inSurfID == 1) {
329  init_x[0] = flat.fire(fXlow[i],fXhigh[i]);
330  init_x[1] = fYhigh[i];
331  init_x[2] = flat.fire(fZlow[i],fZhigh[i]);
332  cosy = -sqrt(flat.fire(0,1));
333  cosz = sqrt(1-cosy*cosy)*TMath::Cos(localPhi);
334  cosx = sqrt(1-cosy*cosy)*TMath::Sin(localPhi);
335  localTheta = TMath::ACos(-cosy);
336  return;
337  }
338 
339  //incident from bottom(3):
340  else if (inSurfID == 3) {
341  init_x[0] = flat.fire(fXlow[i],fXhigh[i]);
342  init_x[1] = fYlow[i];
343  init_x[2] = flat.fire(fZlow[i],fZhigh[i]);
344  cosy = sqrt(flat.fire(0,1));
345  cosz = sqrt(1-cosy*cosy)*TMath::Cos(localPhi);
346  cosx = sqrt(1-cosy*cosy)*TMath::Sin(localPhi);
347  localTheta = TMath::ACos(cosy);
348  return;
349  }
350 
351  //incident from left(4):
352  else if (inSurfID == 4) {
353  init_x[0] = fXlow[i];
354  init_x[1] = flat.fire(fYlow[i],fYhigh[i]);
355  init_x[2] = flat.fire(fZlow[i],fZhigh[i]);
356  cosx = sqrt(flat.fire(0,1));
357  cosy = sqrt(1-cosx*cosx)*TMath::Cos(localPhi);
358  cosz = sqrt(1-cosx*cosx)*TMath::Sin(localPhi);
359  localTheta = TMath::ACos(cosx);
360  return;
361  }
362 
363  //incident from right(2):
364  else if (inSurfID == 2) {
365  init_x[0] = fXhigh[i];
366  init_x[1] = flat.fire(fYlow[i],fYhigh[i]);
367  init_x[2] = flat.fire(fZlow[i],fZhigh[i]);
368  cosx = -sqrt(flat.fire(0,1));
369  cosy = sqrt(1-cosx*cosx)*TMath::Cos(localPhi);
370  cosz = sqrt(1-cosx*cosx)*TMath::Sin(localPhi);
371  localTheta = TMath::ACos(-cosx);
372  return;
373  }
374  }
375 
377  //use http://nusoft.fnal.gov/nova/novasoft/doxygen/html/classevgen_1_1EventGeneratorTest.html#a1a77bf812bda6f7c03012ea11d1c2311
378  if(inSurfID!=2){
379  if(projectToSurface(0,fXhigh[i])){
380  if(fYlow[i] <=final_x[1] && final_x[1]<=fYhigh[i] && fZlow[i]<=final_x[2] && final_x[2]<=fZhigh[i]){
381  outSurfID = 2;
382  return true;
383  }
384  }
385  }
386  if(inSurfID!=4){
387  if( projectToSurface(0,fXlow[i])){
388  if(fYlow[i] <=final_x[1] && final_x[1]<=fYhigh[i] && fZlow[i]<=final_x[2] && final_x[2]<=fZhigh[i]){
389  outSurfID = 4;
390  return true;
391  }
392  }
393  }
394  if(inSurfID!=1){
395  if(projectToSurface(1,fYhigh[i])){
396  if(fXlow[i] <=final_x[0] && final_x[0]<=fXhigh[i] && fZlow[i]<=final_x[2] && final_x[2]<=fZhigh[i]){
397  outSurfID = 1;
398  return true;
399  }
400  }
401  }
402  if(inSurfID!=3){
403  if(projectToSurface(1,fYlow[i])){
404  if(fXlow[i] <=final_x[0] && final_x[0]<=fXhigh[i] && fZlow[i]<=final_x[2] && final_x[2]<=fZhigh[i]){
405  outSurfID = 3;
406  return true;
407  }
408  }
409  }
410  if(inSurfID!=0){
411  if(projectToSurface(2,fZhigh[i])){
412  if(fXlow[i] <=final_x[0] && final_x[0]<=fXhigh[i] && fYlow[i]<=final_x[1] && final_x[1]<=fYhigh[i]){
413  outSurfID = 0;
414  return true;
415  }
416  }
417  }
418  if(inSurfID!=5){
419  if(projectToSurface(2,fZlow[i])){
420  if(fXlow[i] <=final_x[0] && final_x[0]<=fXhigh[i] && fYlow[i]<=final_x[1] && final_x[1]<=fYhigh[i]){
421  outSurfID = 5;
422  return true;
423  }
424  }
425  }
426  return false;
427  }
428  bool IsoGen::projectToSurface(int axis, double surfaceLoc){
429  double momDir = 0, posDir = 0;
430  if(axis == 0){
431  momDir = cosx;
432  }
433  else if(axis == 1){
434  momDir = cosy;
435  }
436  else if(axis == 2){
437  momDir = cosz;
438  }
439 
440  posDir = init_x[axis];
441  double length1Dim = (surfaceLoc - posDir);
442  if(TMath::Abs(momDir)>0.){
443  length1Dim/=momDir;
444  final_x[0] = init_x[0] + length1Dim*cosx;
445  final_x[1] = init_x[1] + length1Dim*cosy;
446  final_x[2] = init_x[2] + length1Dim*cosz;
447  return true;
448  }
449  else{
450  return false;
451  }
452  }
453 
454 
455 } //end namespace
456 
457 ////////////////////////////////////////////////////////////////////////
458 namespace evgen
459 {
461 }
462 ////////////////////////////////////////////////////////////////////////
virtual ~IsoGen()
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
std::vector< double > fZhigh
Up limit of z position (cm)
double init_x[4]
const char * p
Definition: xmltok.h:285
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
T sqrt(T number)
Definition: d0nt_math.hpp:156
A module to check the results from the Monte Carlo generator.
void reconfigure(const fhicl::ParameterSet &pset)
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
std::vector< double > fYlow
Low Limit of y position (cm)
base_engine_t & createEngine(seed_t seed)
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
void setInSurfID(CLHEP::RandFlat &flat, int const &i)
void beginRun(art::Run &run)
std::vector< double > fSigmaT
Variation in t position (ns)
Definition: Run.h:31
void produce(art::Event &evt)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
unsigned int seed
Definition: runWimpSim.h:102
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
single particles thrown at the detector
Definition: MCTruth.h:26
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
std::vector< double > fP0
Central momentum (GeV/c) to generate.
const double j
Definition: BetheBloch.cxx:29
bool getOutSurfPosition(int i)
std::vector< double > fYhigh
Up limit of y position (cm)
void setDirection(CLHEP::RandFlat &flat, int const &i)
std::vector< double > fXhigh
Up limit of x position (cm)
static constexpr Double_t gauss
Definition: Munits.h:360
std::vector< double > fSigmaP
Variation in momenta (GeV/c)
Definition: run.py:1
std::vector< double > fPL
momentum list (GeV/c) to generate
T * make(ARGS...args) const
double final_x[3]
void Sample(simb::MCTruth &truth)
TTree * isoGen_Tree
bool projectToSurface(int axis, double surfaceLoc)
std::vector< int > fPDist
How to distribute momenta (0=uniform, 1=gaussian)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
std::vector< double > fXlow
Low limit of x position (cm)
Event generator information.
Definition: MCTruth.h:32
IsoGen(fhicl::ParameterSet const &pset)
std::vector< double > fT0
Central t position (ns)
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
std::vector< int > fPDG
PDG code of particles to generate.
Encapsulate the geometry of one entire detector (near, far, ndos)
long fireInt(long n)
std::string FileBaseName() const
std::vector< double > fZlow
Low Limit of z position (cm)
enum BeamMode string