LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PhotonGen_module.cc
Go to the documentation of this file.
1 // Plugin Type: producer
3 // File: PhotonGen_module.cc
4 // Description:
5 // Produce photons at the vertex uniformly distributed in the active volume
6 // Oct. 20, 2020 by Mu Wei wmu@fnal.gov 2020
8 
9 // C++ includes.
10 #include <cmath>
11 #include <fstream>
12 #include <memory>
13 #include <string>
14 
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <time.h>
18 
19 // ART includes
25 #include "art_root_io/TFileService.h"
27 #include "cetlib_except/exception.h"
28 #include "fhiclcpp/ParameterSet.h"
30 
31 // art extensions
35 
36 // lar includes
39 
40 #include "TLorentzVector.h"
41 #include "TTree.h"
42 #include "TVector3.h"
43 
44 #include "CLHEP/Random/RandFlat.h"
45 #include "CLHEP/Random/RandGaussQ.h"
46 
47 namespace evgen {
48  class PhotonGen : public art::EDProducer {
49  public:
50  explicit PhotonGen(fhicl::ParameterSet const& pset);
51  void produce(art::Event& evt);
52  void beginRun(art::Run& run);
53 
54  private:
55  void Sample(simb::MCTruth& truth);
56 
57  //Flags
58  static const int kUNIF = 0;
59  static const int kGAUS = 1;
60 
61  //TTree to keep track of where particles have been shot
62  bool fFillTree; // Do we want to create a TTree of shot particles?
63  TTree* fPhotonGen;
64  TLorentzVector fShotPos;
65  TLorentzVector fShotMom;
66  Int_t fEvID;
67 
68  //Scan simulation
69  bool fScan; // Do we want to scan x, y, z or just simulate a random point?
70  double fPx;
71  double fPy; // Fixed coordinate
72  double fPz;
73 
74  //Parameters used to shoot in distributions
75  int fPosDist; //
76  int fTDist; // Random distributions to use : 1= gauss, 0= uniform
77  int fPDist; //
78  double fX; // central x position of source
79  double fY; // central y position of source
80  double fZ; // central z position of source
81  double fT; // central t position of source
82  double fSigmaT; // t width
83  double fP; // central momentm of photon
84  double fSigmaP; // mom width;
85 
86  // Number of photons per event
87  int fN; // number of photons per event
88 
89  CLHEP::HepRandomEngine& fEngine;
90 
91  //Boundaries of the detector
92  double fXmin;
93  double fXmax;
94  double fYmin;
95  double fYmax;
96  double fZmin;
97  double fZmax;
98  };
99 
100  //----------------------------------------------------------------
102  : art::EDProducer{pset}
103  , fFillTree{pset.get<bool>("FillTree")}
104  , fScan{pset.get<bool>("Scan")}
105  , fPosDist{pset.get<int>("PosDist")}
106  , fTDist{pset.get<int>("TDist")}
107  , fPDist{pset.get<int>("PDist")}
108  , fT{pset.get<double>("T0")}
109  , fSigmaT{pset.get<double>("SigmaT")}
110  , fP{pset.get<double>("P")}
111  , fSigmaP{pset.get<double>("SigmaP")}
112  , fN{pset.get<int>("N")}
113  , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
114  pset,
115  "Seed"))
116  {
117  produces<sumdata::RunData, art::InRun>();
118  produces<std::vector<simb::MCTruth>>();
119 
120  if (fFillTree) {
122  fPhotonGen = tfs->make<TTree>("PhGen", "PhGen");
123  fPhotonGen->Branch("X", &(fShotPos[0]), "X/D");
124  fPhotonGen->Branch("Y", &(fShotPos[1]), "Y/D");
125  fPhotonGen->Branch("Z", &(fShotPos[2]), "Z/D");
126  fPhotonGen->Branch("T", &(fShotPos[3]), "T/D");
127  fPhotonGen->Branch("PX", &(fShotMom[0]), "PX/D");
128  fPhotonGen->Branch("PY", &(fShotMom[1]), "PY/D");
129  fPhotonGen->Branch("PZ", &(fShotMom[2]), "PZ/D");
130  fPhotonGen->Branch("PT", &(fShotMom[3]), "PT/D");
131  fPhotonGen->Branch("EventID", &fEvID, "EventID/I");
132  }
133 
134  if (fScan) {
135  fPx = pset.get<double>("X");
136  fPy = pset.get<double>("Y");
137  fPz = pset.get<double>("Z");
138  std::cout << "Will generate photons from 3 points." << std::endl;
139  }
140  }
141 
142  //____________________________________________________________________________
144  {
146  std::cout << "Number of optical detector: " << int(geo->Cryostat().NOpDet()) << std::endl;
147 
148  auto const CryoBounds = geo->Cryostat().Boundaries();
149  fXmin = CryoBounds.MinX();
150  fXmax = CryoBounds.MaxX();
151  fYmin = CryoBounds.MinY();
152  fYmax = CryoBounds.MaxY();
153  fZmin = CryoBounds.MinZ();
154  fZmax = CryoBounds.MaxZ();
155  std::cout << "Cryo Boundaries:" << std::endl;
156  std::cout << "Xmin: " << fXmin << " Xmax: " << fXmax << " Ymin: " << fYmin << " Ymax: " << fYmax
157  << " Zmin: " << fZmin << " Zmax: " << fZmax << std::endl;
158  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
159  }
160 
161  //----------------------------------------------------------------
163  {
164  std::mt19937 rng;
165  rng.seed(std::random_device()());
166  std::uniform_real_distribution<double> distX(fXmin, fXmax);
167  std::uniform_real_distribution<double> distY(fYmin, fYmax);
168  std::uniform_real_distribution<double> distZ(fZmin, fZmax);
169  std::uniform_real_distribution<double> width(-2.0, 2.0);
170 
171  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
172  simb::MCTruth truth;
174 
175  if (fScan) {
176  fX = distX(rng);
177  fY = fPy + width(rng);
178  fZ = fPz + width(rng);
179  Sample(truth);
180  }
181  else {
182  fX = distX(rng);
183  fY = distY(rng);
184  fZ = distZ(rng);
185 
186  Sample(truth);
187  }
188 
189  truthcol->push_back(truth);
190  evt.put(std::move(truthcol));
191  }
192 
194  {
195  std::cout << "Photons Shooting at the Position: " << fX << " " << fY << " " << fZ << std::endl;
196 
197  CLHEP::RandFlat flat(fEngine);
198  CLHEP::RandGaussQ gauss(fEngine);
199 
200  for (int j = 0; j != fN; ++j) {
201  TVector3 pos;
202  pos[0] = fX;
203  pos[1] = fY;
204  pos[2] = fZ;
205 
206  double time;
207  if (fTDist == kGAUS) { time = gauss.fire(fT, fSigmaT); }
208  else {
209  time = fT + fSigmaT * (2.0 * flat.fire() - 1.0);
210  }
211  fShotPos = TLorentzVector(pos[0], pos[1], pos[2], time);
212 
213  //momentum (supplied in eV, convert to GeV)
214  double p = fP;
215  if (fPDist == kGAUS) { p = gauss.fire(fP, fSigmaP); }
216  else {
217  p = fP + fSigmaP * (2.0 * flat.fire() - 1.0);
218  }
219  p /= 1000000000.;
220 
221  //angles
222  double costh = 2 * flat.fire() - 1;
223  double sinth = std::sqrt(1 - costh * costh);
224  double phi = 2 * M_PI * flat.fire();
225 
226  //momentum 4-vector
227  fShotMom = TLorentzVector(p * sinth * cos(phi), p * sinth * sin(phi), p * costh, p);
228 
229  int trackid =
230  -1 *
231  (j + 1); // set track id to negative as these are all primary particles and have id <= 0
232  int PDG = 0; // optical photons have PDG 0
233  simb::MCParticle particle(trackid, PDG, "primary");
235 
236  if (fFillTree) { fPhotonGen->Fill(); }
237 
238  mct.Add(particle);
239  }
240  }
241 }
242 
base_engine_t & createEngine(seed_t seed)
PhotonGen(fhicl::ParameterSet const &pset)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:258
constexpr auto fullRun()
static const int kUNIF
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
TLorentzVector fShotMom
void beginRun(art::Run &run)
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
void Sample(simb::MCTruth &truth)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
CryostatGeo const & Cryostat(CryostatID const &cryoid=cryostat_zero) const
Returns the specified cryostat.
Particle class.
Definition: Run.h:37
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
single particles thrown at the detector
Definition: MCTruth.h:26
CLHEP::HepRandomEngine & fEngine
TLorentzVector fShotPos
static const int kGAUS
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
unsigned int NOpDet() const
Number of optical detectors in this TPC.
Definition: CryostatGeo.h:321
Definition: MVAAlg.h:12
geo::BoxBoundedGeo const & Boundaries() const
Returns boundaries of the cryostat (in centimetres).
Definition: CryostatGeo.h:114
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:203
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
Namespace collecting geometry-related classes utilities.
void produce(art::Event &evt)
Event Generation using GENIE, cosmics or single particles.
static ServiceHandle< RandomNumberGenerator > & rng()
art framework interface to geometry description