All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MUSUN_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file MUSUN_module.cc
3 /// \brief Generator for underground muon propagation
4 ///
5 /// Module designed to propagate muons underground
6 ///
7 /// For a description of how to use the module see DUNE DocDB
8 /// It is highly recommended that you read it before use.....
9 ///
10 /// Original MUSUN code was written by Vitaly A. Kudryavtsev (University of Sheffield).
11 /// Conversion from Fortran to C was done by Kareem Kazkaz (LLNL) with help from David Woodward (Sheffield)
12 /// Default slant depths and surface profile work done Martin Richardson (Sheffield)
13 /// Notable contributions to surface map profiling also made by Chao Zhang (USD) and Jeff de Jong (Oxford)
14 /// Interfaced with LArSoft by Karl Warburton (University of Sheffield) with necessary changes.
15 /// ------------START OF NECCESSARY CHANGES---------------
16 /// Changing variables to fcl parameters
17 /// Restructuring to fit LArSoft design
18 /// Co-ordinate transformation... y -> x, x -> z, z -> y ????SHOULD BE x -> -z AS PER VITALY????
19 /// Restore cavern angle convention to Vitaly's definition.
20 /// ------------END OF NECCESSARY CHANGES-----------------
21 /// \author k.warburton@sheffield.ac.uk
22 ///
23 /// Further Notes - Taken from Kareem Kazkaz;
24 ///
25 /// This C++ code is a port of the musun-surf.f and test-musun-surf.f code written
26 /// by Vitaly Kudryavtsev of the University of Sheffield. It generates muons with
27 /// energy, position, and direction specific to the Davis cavern at the Sanford
28 /// Underground Research Facility in Lead, South Dakota.
29 ///
30 /// This C++ code was ported by Kareem Kazkaz, kareem@llnl.gov, (925) 422-7208
31 ///
32 /// Here are the notes from Vitaly:
33 ///
34 ///c The code samples single atmospheric muons at the SURF
35 ///c underground laboratory (Davis' cavern)
36 ///c (taking into account the slant depth distribution)
37 ///c in the energy range E1-E2, zenith angle range theta1-theta2 (0-90 degrees)
38 ///c and azimuthal angle range phi1-phi2 (0-360 degrees).
39 ///c At present only the following ranges of parameters are supported:
40 ///c E1 = 1 GeV, E2 = 10^6 GeV, theta1 = 0, theta2 = 90 deg, phi1 = 0, phi2 = 360 deg.
41 ///c
42 ///c Program uses muon energy spectra at various depths and zenith
43 ///c angles obtained with MUSIC code for muon propagation and Gaisser's
44 ///c formula for muon spectrum at sea level
45 ///c (T.K.Gaisser, Cosmic Rays and Particle Physics, Cambridge
46 ///c University Press, 1990) modified for large zenith angles and
47 ///c prompt muon flux with normalisation and spectral index
48 ///c that fit LVD data: gamma = 2.77, A = 0.14.
49 ///c Density of rock is assumed to be 2.70 g/cm^3 but can be changed
50 ///c during initialisation (previous step, ask the author).
51 ///c
52 ///c Muon flux through a sphere (Chao's coordinates) = 6.33x10^(-9) cm^(-2) s^(-1) (gamma=2.77) - old
53 ///c Muon flux through a sphere (Martin's coordinates) = 6.16x10^(-9) cm^(-2) s^(-1) (gamma=2.77) - new
54 ///c Muon flux through the cuboid (30x22x24 m^3) = 0.0588 s^(-1) (gamma=2.77)
55 ///c
56 ///c Note: the muon spectrum at sea level does not take into account
57 ///c the change of the primary spectrum slope at and above the knee
58 ///c region (3*10^15-10^16 eV).
59 ///c
60 ///c Program uses the tables of muon energy spectra at various
61 ///c zenith and azimuthal angles at SURF
62 ///c calculated with the muon propagation code MUSIC and the
63 ///c angular distribution of muon intensities at SURF (4850 ft level).
64 ///c
65 ///c Coordinate system for the muon intensities
66 ///c is determined by the mountain profile provided
67 ///c by Chao Zhang (USD, South Dakota): x-axis is pointing to the East.
68 ///c Muons are sampled on a surface of a rectangular parallelepiped,
69 ///c the planes of which are parallel to the walls of the cavern.
70 ///c The long side of the cavern is pointing at 6.6 deg from the North
71 ///c to the East (or 90 - 6.6 deg from the East to the North).
72 ///c Muon coordinates and direction cosines are then given in the
73 ///c coordinate system related to the cavern with x-axis
74 ///c pointing along the long side of the cavern at 6.6 deg from the
75 ///c North to the East.
76 ///c The angle phi is measured from the positive x-axis towards
77 ///c positive y-axis (0-360 deg).
78 ///c Z-axis is pointing upwards.
79 ///
80 /// Further notes taken from Vitaly's original code, as written for LBNE
81 /// (note the differing definition for theta, which is used here,
82 /// theta is measured from East to South, not North to East.
83 ///c
84 ///c The code samples single atmospheric muons at the LBNE
85 ///c underground site (taking into account the slant depth
86 ///c distribution calculated by Martin Richardson, Sheffield)
87 ///c in the energy range E1-E2, zenith angle range theta1-theta2 (0-90 degrees)
88 ///c and azimuthal angle range phi1-phi2 (0-360 degrees).
89 ///c At present only the following ranges of parameters are supported:
90 ///c E1 = 1 GeV, E2 = 10^6 GeV, theta1 = 0, theta2 = 90 deg, phi1 = 0, phi2 = 360 deg.
91 ///c
92 ///c The LBNE far detector site has coordinates:
93 ///c Latitude = 44° 20' 45.21" N
94 ///c Longitude = 103° 45' 16.13" W
95 ///c Elevation = 355.8 ft (108.4 m)
96 ///c from e-mail from Virgil Bocean on 30 May 2013
97 ///c confirmed by Tracy Lundin on 16 December 2013
98 ///c
99 ///c For slant depth distribution and muon intensities the azimuthal
100 ///c angle phi is calculated from East to North.
101 ///c The long side of the cavern is assumed to be pointing to the
102 ///c beam (Fermilab) at an angle of 7 deg to the South from the
103 ///c East.
104 ///c Muon direction cosines are given in the coordinate system
105 ///c related to the cavern with positive x-axis pointing to the beam
106 ///c so direction cosine wrt x-axis is -1 if a muon is coming in the
107 ///c same direction (wrt to x-axis) as the neutrino beam.
108 ///c
109 ///c
110 ///c Program uses muon energy spectra at various depths and zenith
111 ///c angles obtained with MUSIC code for muon propagation and Gaisser's
112 ///c formula for muon spectrum at sea level
113 ///c (T.K.Gaisser, Cosmic Rays and Particle Physics, Cambridge
114 ///c University Press, 1990) modified for large zenith angles and
115 ///c prompt muon flux with parameters from the best fit to the LVD data.
116 ///c
117 ///c Muon flux through a sphere = 6.65x10^(-9) cm^(-2) s^(-1) (gamma=2.77) - old
118 ///c Muon flux through a sphere = 6.65x10^(-9) cm^(-2) s^(-1) (gamma=2.77) - new
119 ///c Muon flux through the cuboid (100x40x50 m^3) = 0.3074 s^(-1) (gamma=2.77)
120 ///c
121 ///c Note: the muon spectrum at sea level does not take into account
122 ///c the change of the primary spectrum slope at and above the knee
123 ///c region (3*10^15-10^16 eV).
124 ///c
125 ///c Program uses the tables of muon energy spectra at various
126 ///c zenith and azimuthal angles at DUSEL
127 ///c calculated with the muon propagation code MUSIC and the
128 ///c angular distribution of muon intensities at 4850 ft level
129 ///c (location of the LBNE far detector).
130 ///c
131 ///c Coordinate system is determined by the mountain profile provided
132 ///c by Jeff de Jong (Cambridge) and Chao Zhang (USD, South Dakota).
133 ///c This is done to allow the simulation of muons on a surface of
134 ///c a rectangular parallelepiped, the planes of which are parallel
135 ///c to the walls of the laboratory halls.
136 ///c The angle phi in the slant depth distribution file is measured
137 ///c from the positive x-axis towards positive y-axis.
138 ///c Z-axis is pointing upwards.
139 ///c Note that I use here the azimuthal angle range from 0 to 360 deg.
140 ///c
141 ////////////////////////////////////////////////////////////////////////
142 
143 // C++ includes.
144 #include <cmath>
145 #include <cstdlib>
146 #include <exception>
147 #include <fstream>
148 #include <iostream>
149 #include <memory>
150 #include <sstream>
151 #include <string>
152 
153 // Framework includes
154 #include "art/Framework/Core/EDProducer.h"
155 #include "art/Framework/Core/ModuleMacros.h"
156 #include "art/Framework/Principal/Event.h"
157 #include "art/Framework/Principal/Run.h"
158 #include "art/Framework/Services/Registry/ServiceHandle.h"
159 #include "art_root_io/TFileService.h"
160 #include "cetlib_except/exception.h"
161 #include "fhiclcpp/ParameterSet.h"
162 #include "messagefacility/MessageLogger/MessageLogger.h"
163 
164 // art extensions
165 #include "nurandom/RandomUtils/NuRandomService.h"
166 
167 // nusimdata includes
168 #include "nusimdata/SimulationBase/MCParticle.h"
169 #include "nusimdata/SimulationBase/MCTruth.h"
170 
171 // lar includes
174 
175 #include "TDatabasePDG.h"
176 #include "TTree.h"
177 
178 #include "CLHEP/Random/RandFlat.h"
179 #include "CLHEP/Random/RandGaussQ.h"
180 
181 // for c2: NTPCs is no longer used
182 //const int NTPCs = 300;
183 
184 namespace simb {
185  class MCTruth;
186 }
187 
188 namespace evgen {
189 
190  /// module to produce single or multiple specified particles in the detector
191  class MUSUN : public art::EDProducer {
192 
193  public:
194  explicit MUSUN(fhicl::ParameterSet const& pset);
195 
196  // This is called for each event.
197  void produce(art::Event& evt);
198  void beginJob();
199  void beginRun(art::Run& run);
200  void endRun(art::Run& run);
201 
202  private:
203  void SampleOne(unsigned int i, simb::MCTruth& mct, CLHEP::HepRandomEngine& engine);
204 
205  void initialization(double theta1,
206  double theta2,
207  double phi1,
208  double phi2,
209  int figflag,
210  double s_hor,
211  double s_ver1,
212  double s_ver2,
213  double& FI);
214 
215  void sampling(double& E,
216  double& theta,
217  double& phi,
218  double& dep,
219  CLHEP::HepRandomEngine& engine);
220 
221  static const int kGAUS = 1;
222 
223  CLHEP::HepRandomEngine& fEngine; ///< art-managed random-number engine
224 
225  int fPDG; ///< PDG code of particles to generate
226  double fChargeRatio; ///< Charge ratio of particle / anti-particle
227 
228  std::string fInputDir; ///< Input Directory
229  std::string fInputFile1; ///< Input File 1
230  std::string fInputFile2; ///< Input File 2
231  std::string fInputFile3; ///< Input File 3
232 
233  double fCavernAngle; ///< Angle of the detector from the North to the East.
234  double fRockDensity; ///< Default rock density is 2.70 g cm-3. If this is
235  ///< changed then the three input files need to be
236  ///< remade. If there is a desire for this contact
237  ///< Vitaly Kudryavtsev at V.Kudryavtsev@shef.ac.uk
238 
239  double fEmin; ///< Minimum Kinetic Energy (GeV)
240  double fEmax; ///< Maximum Kinetic Energy (GeV)
241 
242  double fThetamin; ///< Minimum theta
243  double fThetamax; ///< Maximum theta
244  double fPhimin; ///< Minimum phi
245  double fPhimax; ///< Maximum phi
246 
247  int figflag; ///< If want sampled from sphere or parallelepiped
248  double fXmin; ///< Minimum X position
249  double fYmin; ///< Minimum Y position
250  double fZmin; ///< Minimum Z position
251  double fXmax; ///< Maximum X position
252  double fYmax; ///< Maximum Y position
253  double fZmax; ///< Maximum Z position
254 
255  double fT0; ///< Central t position (ns) in world coordinates
256  double fSigmaT; ///< Variation in t position (ns)
257  int fTDist; ///< How to distribute t (gaus, or uniform)
258 
259  //Define TFS histograms.....
260  /*
261  TH1D* hPDGCode;
262  TH1D* hPositionX;
263  TH1D* hPositionY;
264  TH1D* hPositionZ;
265  TH1D* hTime;
266  TH1D* hMomentumHigh;
267  TH1D* hMomentum;
268  TH1D* hEnergyHigh;
269  TH1D* hEnergy;
270  TH1D* hDepth;
271  TH1D* hDirCosineX;
272  TH1D* hDirCosineY;
273  TH1D* hDirCosineZ;
274  TH1D* hTheta;
275  TH1D* hPhi;
276  */
277  int PdgCode;
278  double Energy, phi, theta, dep, Time;
279  double Momentum, px0, py0, pz0;
280  double x0, y0, z0, cx, cy, cz;
281 
282  //Define some variables....
283  double FI = 0.;
284  double s_hor = 0.;
285  double s_ver1 = 0.;
286  double s_ver2 = 0.;
287 
288  double spmu[121][62][51];
289  double fnmu[32401];
290  double depth[360][91];
291  double fmu[360][91];
292  // for c2: e1 and e2 are unused
293  //double e1, e2, the1, the2, ph1, ph2;
294  double the1, the2, ph1, ph2;
295  double se = 0.;
296  double st = 0.;
297  double sp = 0.;
298  double sd = 0.;
299 
300  unsigned int NEvents = 0;
301 
302  // TTree
303  TTree* fTree;
304  /*
305  // TTree for CryoPos
306  TTree* fCryos;
307  int NumTPCs;
308  double TPCMinX[NTPCs];
309  double TPCMaxX[NTPCs];
310  double TPCMinY[NTPCs];
311  double TPCMaxY[NTPCs];
312  double TPCMinZ[NTPCs];
313  double TPCMaxZ[NTPCs];
314  double CryoSize[6];
315  double DetHall[6];
316  */
317  };
318 }
319 
320 namespace evgen {
321 
322  //____________________________________________________________________________
323  MUSUN::MUSUN(fhicl::ParameterSet const& pset)
324  : art::EDProducer{pset}
325  // create a default random engine; obtain the random seed from NuRandomService,
326  // unless overridden in configuration with key "Seed"
327  , fEngine(art::ServiceHandle<rndm::NuRandomService> {}->createEngine(*this, pset, "Seed"))
328  , fPDG{pset.get<int>("PDG")}
329  , fChargeRatio{pset.get<double>("ChargeRatio")}
330  , fInputDir{pset.get<std::string>("InputDir")}
331  , fInputFile1{pset.get<std::string>("InputFile1")}
332  , fInputFile2{pset.get<std::string>("InputFile2")}
333  , fInputFile3{pset.get<std::string>("InputFile3")}
334  , fCavernAngle{pset.get<double>("CavernAngle")}
335  , fRockDensity{pset.get<double>("RockDensity")}
336  , fEmin{pset.get<double>("Emin")}
337  , fEmax{pset.get<double>("Emax")}
338  , fThetamin{pset.get<double>("Thetamin")}
339  , fThetamax{pset.get<double>("Thetamax")}
340  , fPhimin{pset.get<double>("Phimin")}
341  , fPhimax{pset.get<double>("Phimax")}
342  , figflag{pset.get<int>("igflag")}
343  , fXmin{pset.get<double>("Xmin")}
344  , fYmin{pset.get<double>("Ymin")}
345  , fZmin{pset.get<double>("Zmin")}
346  , fXmax{pset.get<double>("Xmax")}
347  , fYmax{pset.get<double>("Ymax")}
348  , fZmax{pset.get<double>("Zmax")}
349  , fT0{pset.get<double>("T0")}
350  , fSigmaT{pset.get<double>("SigmaT")}
351  , fTDist{pset.get<int>("TDist")}
352  {
353  produces<std::vector<simb::MCTruth>>();
354  produces<sumdata::RunData, art::InRun>();
355  }
356 
357  ////////////////////////////////////////////////////////////////////////////////
358  // Begin Job
359  ////////////////////////////////////////////////////////////////////////////////
360  void
362  {
363  // Make the Histograms....
364  art::ServiceHandle<art::TFileService const> tfs;
365  /*
366  hPDGCode = tfs->make<TH1D>("hPDGCode" ,"PDG Code" ,30 , -15 , 15 );
367  hPositionX = tfs->make<TH1D>("hPositionX" ,"Position (cm)" ,500, ( fXmin - 10 ), ( fXmax + 10 ) );
368  hPositionY = tfs->make<TH1D>("hPositionY" ,"Position (cm)" ,500, ( fYmin - 10 ), ( fYmax + 10 ) );
369  hPositionZ = tfs->make<TH1D>("hPositionZ" ,"Position (cm)" ,500, ( fZmin - 10 ), ( fZmax + 10 ) );
370  hTime = tfs->make<TH1D>("hTime" ,"Time (s)" ,500, 0 , 1e6 );
371  hMomentumHigh = tfs->make<TH1D>("hMomentumHigh","Momentum (GeV)",500, fEmin, fEmax );
372  hMomentum = tfs->make<TH1D>("hMomentum" ,"Momentum (GeV)",500, fEmin, fEmin+1e3 );
373  hEnergyHigh = tfs->make<TH1D>("hEnergyHigh" ,"Energy (GeV)" ,500, fEmin, fEmax );
374  hEnergy = tfs->make<TH1D>("hEnergy" ,"Energy (GeV)" ,500, fEmin, fEmin+1e3 );
375  hDepth = tfs->make<TH1D>("hDepth" ,"Depth (m)" ,800, 0 , 14000 );
376 
377  hDirCosineX = tfs->make<TH1D>("hDirCosineX","Normalised Direction cosine",500, -1, 1 );
378  hDirCosineY = tfs->make<TH1D>("hDirCosineY","Normalised Direction cosine",500, -1, 1 );
379  hDirCosineZ = tfs->make<TH1D>("hDirCosineZ","Normalised Direction cosine",500, -1, 1 );
380 
381  hTheta = tfs->make<TH1D>("hTheta" ,"Angle (degrees)",500, 0, 90 );
382  hPhi = tfs->make<TH1D>("hPhi" ,"Angle (degrees)",500, 0, 365 );
383  */
384  fTree = tfs->make<TTree>("Generator", "analysis tree");
385  fTree->Branch("particleID", &PdgCode, "particleID/I");
386  fTree->Branch("energy", &Energy, "energy/D");
387  fTree->Branch("time", &Time, "Time/D");
388  fTree->Branch("posX", &x0, "posX/D");
389  fTree->Branch("posY", &y0, "posY/D");
390  fTree->Branch("posZ", &z0, "posZ/D");
391  fTree->Branch("cosX", &cx, "cosX/D");
392  fTree->Branch("cosY", &cy, "cosY/D");
393  fTree->Branch("cosZ", &cz, "cosZ/D");
394  fTree->Branch("theta", &theta, "theta/D");
395  fTree->Branch("phi", &phi, "phi/D");
396  fTree->Branch("depth", &dep, "dep/D");
397  /*
398  fCryos = tfs->make<TTree>("CryoSizes","cryo tree");
399  fCryos->Branch("NumTPCs" , &NumTPCs , "NumTPCs/I" );
400  fCryos->Branch("TPCMinX" , &TPCMinX , "TPCMinX[NumTPCs]/D");
401  fCryos->Branch("TPCMaxX" , &TPCMaxX , "TPCMaxX[NumTPCs]/D");
402  fCryos->Branch("TPCMinY" , &TPCMinY , "TPCMinY[NumTPCs]/D");
403  fCryos->Branch("TPCMaxY" , &TPCMaxY , "TPCMaxY[NumTPCs]/D");
404  fCryos->Branch("TPCMinZ" , &TPCMinZ , "TPCMinZ[NumTPCs]/D");
405  fCryos->Branch("TPCMaxZ" , &TPCMaxZ , "TPCMaxZ[NumTPCs]/D");
406  fCryos->Branch("CryoSize", &CryoSize, "CryoSize[6]/D" );
407  fCryos->Branch("DetHall" , &DetHall , "DetHall[6]/D" );
408  */
409  }
410 
411  ////////////////////////////////////////////////////////////////////////////////
412  // Begin Run
413  ////////////////////////////////////////////////////////////////////////////////
414  void
415  MUSUN::beginRun(art::Run& run)
416  {
417  // Check fcl parameters were set correctly
418  if (fThetamax > 90.5)
419  throw cet::exception("MUSUNGen")
420  << "\nThetamax has to be less than " << M_PI / 2 << ", but was entered as " << fThetamax
421  << ", this causes an error so leaving program now...\n\n";
422  if (fThetamin < 0)
423  throw cet::exception("MUSUNGen")
424  << "\nThetamin has to be more than 0, but was entered as " << fThetamin
425  << ", this causes an error so leaving program now...\n\n";
426  if (fThetamax < fThetamin)
427  throw cet::exception("MUSUNGen") << "\nMinimum angle is bigger than maximum angle....causes "
428  "an error so leaving program now....\n\n";
429  if (fPhimax > 360.5)
430  throw cet::exception("MUSUNGen")
431  << "\nPhimax has to be less than " << 2 * M_PI << ", but was entered as " << fPhimax
432  << ", this cause an error so leaving program now...\n\n";
433  if (fPhimin < 0)
434  throw cet::exception("MUSUNGen")
435  << "\nPhimin has to be more than 0, but was entered as " << fPhimin
436  << ", this causes an error so leaving program now...\n\n";
437  if (fPhimax < fPhimin)
438  throw cet::exception("MUSUNGen") << "\nMinimum angle is bigger than maximum angle....causes "
439  "an error so leaving program now....\n\n";
440  if (fEmax < fEmin)
441  throw cet::exception("MUSUNGen")
442  << "\nMinimum energy is bigger than maximum energy....causes an error so leaving program "
443  "now....\n\n";
444 
445  // grab the geometry object to see what geometry we are using
446  art::ServiceHandle<geo::Geometry const> geo;
447  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
448 
449  // area of the horizontal plane of the parallelepiped
450  s_hor = (fZmax - fZmin) * (fXmax - fXmin);
451  // area of the vertical plane of the parallelepiped, perpendicular to z-axis
452  s_ver1 = (fXmax - fXmin) * (fYmax - fYmin);
453  // area of the vertical plane of the parallelepiped, perpendicular to x-axis
454  s_ver2 = (fZmax - fZmin) * (fYmax - fYmin);
455 
456  //std::cout << s_hor << " " << s_ver1 << " " << s_ver2 << std::endl;
457 
459 
460  std::cout << "Material - SURF rock" << std::endl;
461  std::cout << "Density = " << fRockDensity << " g/cm^3" << std::endl;
462  std::cout << "Parameters for muon spectrum are from LVD best fit" << std::endl;
463  std::cout << "Muon energy range = " << fEmin << " - " << fEmax << " GeV" << std::endl;
464  std::cout << "Zenith angle range = " << fThetamin << " - " << fThetamax << " degrees"
465  << std::endl;
466  std::cout << "Azimuthal angle range = " << fPhimin << " - " << fPhimax << " degrees"
467  << std::endl;
468  std::cout << "Global intensity = " << FI
469  << " (cm^2 s)^(-1) or s^(-1) (for muons on the surface)" << std::endl;
470  /*
471  NumTPCs = geo->NTPC(0);
472  std::cout << "There are " << NumTPCs << " in cryostat 0" << std::endl;
473  for (unsigned int c=0; c<geo->Ncryostats(); c++) {
474  const geo::CryostatGeo& cryostat=geo->Cryostat(c);
475  geo->CryostatBoundaries( CryoSize, 0 );
476  std::cout << "Cryo bounds " << CryoSize[0] << " "<< CryoSize[1] << " "<< CryoSize[2] << " "<< CryoSize[3] << " "<< CryoSize[4] << " "<< CryoSize[5] << std::endl;
477  for (unsigned int t=0; t<cryostat.NTPC(); t++) {
478  geo::TPCID id;
479  id.Cryostat=c;
480  id.TPC=t;
481  id.isValid=true;
482  const geo::TPCGeo& tpc=cryostat.TPC(id);
483  TPCMinX[t] = tpc.MinX();
484  TPCMaxX[t] = tpc.MaxX();
485  TPCMinY[t] = tpc.MinY();
486  TPCMaxY[t] = tpc.MaxY();
487  TPCMinZ[t] = tpc.MinZ();
488  TPCMaxZ[t] = tpc.MaxZ();
489  std::cout << t << "\t" << TPCMinX[t] << " " << TPCMaxX[t] << " " << TPCMinY[t] << " " << TPCMaxY[t] << " " << TPCMinZ[t] << " " << TPCMaxZ[t] << std::endl;
490  }
491  }
492  fCryos -> Fill();
493  */
494  return;
495  }
496 
497  ////////////////////////////////////////////////////////////////////////////////
498  // End Run
499  ////////////////////////////////////////////////////////////////////////////////
500  void
501  MUSUN::endRun(art::Run& run)
502  {
503  std::cout << "\n\nNumber of muons = " << NEvents << std::endl;
504  std::cout << "Mean muon energy = " << se / NEvents << " GeV" << std::endl;
505  std::cout << "Mean zenith angle (theta) = " << st / NEvents << " degrees" << std::endl;
506  std::cout << "Mean azimuthal angle (phi)= " << sp / NEvents << " degrees" << std::endl;
507  std::cout << "Mean slant depth = " << sd / NEvents << " m w.e." << std::endl;
508  }
509  ////////////////////////////////////////////////////////////////////////////////
510  // Produce
511  ////////////////////////////////////////////////////////////////////////////////
512  void
513  MUSUN::produce(art::Event& evt)
514  {
515  ///unique_ptr allows ownership to be transferred to the art::Event after the put statement
516  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
517 
518  simb::MCTruth truth;
519  truth.SetOrigin(simb::kSingleParticle);
520 
521  ++NEvents;
522 
523  SampleOne(NEvents, truth, fEngine);
524 
525  MF_LOG_DEBUG("MUSUN") << truth;
526 
527  truthcol->push_back(truth);
528 
529  evt.put(std::move(truthcol));
530  }
531  ////////////////////////////////////////////////////////////////////////////////
532  // Sample One
533  ////////////////////////////////////////////////////////////////////////////////
534  // Draw the type, momentum and position of a single particle from the
535  // FCIHL description
536  void
537  MUSUN::SampleOne(unsigned int i, simb::MCTruth& mct, CLHEP::HepRandomEngine& engine)
538  {
539  CLHEP::RandFlat flat(engine);
540  CLHEP::RandGaussQ gauss(engine);
541 
542  Energy = 0;
543  theta = 0;
544  phi = 0;
545  dep = 0;
546  Time = 0;
547 
548  sampling(Energy, theta, phi, dep, engine);
549 
550  theta = theta * M_PI / 180;
551 
552  // changing the angle phi so z-axis is positioned along the long side
553  // of the cavern pointing at 14 deg from the North to the East.
554  // phi += (90. - 14.0);
555  // Want our co-ord rotation going from East to South.
556  phi += fCavernAngle;
557  if (phi >= 360.) phi -= 360.;
558  if (phi < 0) phi += 360.;
559  phi *= M_PI / 180.;
560 
561  // set track id to -i as these are all primary particles and have id <= 0
562  int trackid = -1 * (i + 1);
563  std::string primary("primary");
564 
565  // Work out whether particle/antiparticle, and mass...
566  double m = 0.0;
567  PdgCode = fPDG;
568  double ChargeCheck = 1. / (1 + fChargeRatio);
569  double pdgfire = flat.fire();
570  if (pdgfire < ChargeCheck) PdgCode = -PdgCode;
571 
572  static TDatabasePDG pdgt;
573  TParticlePDG* pdgp = pdgt.GetParticle(PdgCode);
574  if (pdgp) m = pdgp->Mass();
575 
576  //std::cout << pdgfire << " " << ChargeCheck << " " << PdgCode << " " << m << std::endl;
577 
578  // Work out T0...
579  if (fTDist == kGAUS) { Time = gauss.fire(fT0, fSigmaT); }
580  else {
581  Time = fT0 + fSigmaT * (2.0 * flat.fire() - 1.0);
582  }
583 
584  // The minus sign above is for y-axis pointing up, so the y-momentum
585  // is always pointing down
586  cx = +sin(theta) * sin(phi);
587  cy = -cos(theta);
588  cz = +sin(theta) * cos(phi);
589  Momentum = std::sqrt(Energy * Energy - m * m); // Get momentum
590  px0 = Momentum * cx;
591  py0 = Momentum * cy;
592  pz0 = Momentum * cz;
593  TLorentzVector pvec(px0, py0, pz0, Energy);
594 
595  // Muon coordinates
596  double sh1 = s_hor * cos(theta);
597  double sv1 = s_ver1 * sin(theta) * fabs(cos(phi));
598  double sv2 = s_ver2 * sin(theta) * fabs(sin(phi));
599  double ss = sh1 + sv1 + sv2;
600  double xfl1 = flat.fire();
601  if (xfl1 <= sh1 / ss) {
602  x0 = (fXmax - fXmin) * flat.fire() + fXmin;
603  y0 = fYmax;
604  z0 = (fZmax - fZmin) * flat.fire() + fZmin;
605  }
606  else if (xfl1 <= (sh1 + sv1) / ss) {
607  x0 = (fXmax - fXmin) * flat.fire() + fXmin;
608  y0 = (fYmax - fYmin) * flat.fire() + fYmin;
609  if (cz >= 0)
610  z0 = fZmin;
611  else
612  z0 = fZmax;
613  }
614  else {
615  if (cx >= 0)
616  x0 = fXmin;
617  else
618  x0 = fXmax;
619  y0 = (fYmax - fYmin) * flat.fire() + fYmin;
620  z0 = (fZmax - fZmin) * flat.fire() + fZmin;
621  }
622  // Make Lorentz vector for x and t....
623  TLorentzVector pos(x0, y0, z0, Time);
624 
625  // Parameters written to the file muons_surf_v2_test*.dat
626  // nmu - muon sequential number
627  // id_part - muon charge (10 - positive, 11 - negative )
628  // Energy - total muon energy in GeV assuming ultrarelativistic muons
629  // x0, y0, z0 - muon coordinates on the surface of parallelepiped
630  // specified above; x-axis and y-axis are pointing in the
631  // directions such that the angle phi (from the slant depth
632  // distribution files) is measured from x to y. z-axis is
633  // pointing upwards.
634  // cx, cy, cz - direction cosines.
635 
636  simb::MCParticle part(trackid, PdgCode, primary);
637  part.AddTrajectoryPoint(pos, pvec);
638 
639  mct.Add(part);
640 
641  theta = theta * 180 / M_PI;
642  phi = phi * 180 / M_PI;
643 
644  // Sum energies, angles, depth for average outputs.
645  se += Energy;
646  st += theta;
647  sp += phi;
648  sd += dep;
649 
650  // Fill Histograms.....
651  /*
652  hPDGCode ->Fill (PdgCode);
653  hPositionX ->Fill (x0);
654  hPositionY ->Fill (y0);
655  hPositionZ ->Fill (z0);
656  hTime ->Fill (Time);
657  hMomentumHigh ->Fill (Momentum);
658  hMomentum ->Fill (Momentum);
659  hEnergyHigh ->Fill (Energy);
660  hEnergy ->Fill (Energy);
661  hDepth ->Fill (dep);
662  hDirCosineX ->Fill (cx);
663  hDirCosineY ->Fill (cy);
664  hDirCosineZ ->Fill (cz);
665  hTheta ->Fill (theta);
666  hPhi ->Fill (phi);
667  */
668  /*
669  // Write event by event outsputs.....
670  std::cout << "Theta: " << theta << " Phi: " << phi << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
671  std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
672  std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
673  std::cout << "Normalised..." << cx << " " << cy << " " << cz << std::endl;
674  */
675  fTree->Fill();
676  }
677 
678  ////////////////////////////////////////////////////////////////////////////////
679  // initialization
680  ////////////////////////////////////////////////////////////////////////////////
681  void
682  MUSUN::initialization(double theta1,
683  double theta2,
684  double phi1,
685  double phi2,
686  int figflag,
687  double s_hor,
688  double s_ver1,
689  double s_ver2,
690  double& FI)
691  {
692  //
693  // Read in the data files
694  //
695  int lineNumber = 0, index = 0;
696  char inputLine[10000];
697  std::string fROOTfile;
698 
699  for (int a = 0; a < 121; ++a)
700  for (int b = 0; b < 62; ++b)
701  for (int c = 0; c < 50; ++c)
702  spmu[a][b][c] = 0;
703  for (int a = 0; a < 23401; ++a)
704  fnmu[a] = 0;
705  for (int a = 0; a < 360; ++a)
706  for (int b = 0; b < 91; ++b)
707  depth[a][b] = 0;
708  for (int a = 0; a < 360; ++a)
709  for (int b = 0; b < 91; ++b)
710  fmu[a][b] = 0;
711 
712  std::ostringstream File1LocStream;
713  File1LocStream << fInputDir << fInputFile1;
714  std::string File1Loc = File1LocStream.str();
715  cet::search_path sp1("FW_SEARCH_PATH");
716  if (sp1.find_file(fInputFile1, fROOTfile)) File1Loc = fROOTfile;
717  std::ifstream file1(File1Loc.c_str(), std::ios::in);
718  if (!file1.good())
719  throw cet::exception("MUSUNGen")
720  << "\nFile1 " << fInputFile1 << " not found in FW_SEARCH_PATH or at " << fInputDir
721  << "\n\n";
722 
723  while (file1.good()) {
724  //std::cout << "Looking at file 1...." << std::endl;
725  file1.getline(inputLine, 9999);
726  char* token;
727  token = strtok(inputLine, " ");
728  while (token != NULL) {
729  //std::cout << "While loop file 1..." << std::endl;
730  fmu[index][lineNumber] = atof(token);
731  token = strtok(NULL, " ");
732  index++;
733  if (index == 360) {
734  //std::cout << "If statement file 1..." << std::endl;
735  index = 0;
736  lineNumber++;
737  }
738  }
739  }
740  file1.close();
741 
742  std::ostringstream File2LocStream;
743  File2LocStream << fInputDir << fInputFile2;
744  std::string File2Loc = File2LocStream.str();
745  cet::search_path sp2("FW_SEARCH_PATH");
746  if (sp2.find_file(fInputFile2, fROOTfile)) File2Loc = fROOTfile;
747  std::ifstream file2(File2Loc.c_str(), std::ios::binary | std::ios::in);
748  if (!file2.good())
749  throw cet::exception("MUSUNGen")
750  << "\nFile2 " << fInputFile2 << " not found in FW_SEARCH_PATH or at " << fInputDir
751  << "\n\n";
752 
753  int i1 = 0, i2 = 0, i3 = 0;
754  float readVal;
755  while (file2.good()) {
756  //std::cout << "Looking at file 2...." << std::endl;
757  file2.read((char*)(&readVal), sizeof(float));
758  spmu[i1][i2][i3] = readVal;
759  i1++;
760  if (i1 == 121) {
761  //std::cout << "First if statement file 2..." << std::endl;
762  i2++;
763  i1 = 0;
764  }
765  if (i2 == 62) {
766  //std::cout << "Second if statement file 2..." << std::endl;
767  i3++;
768  i2 = 0;
769  }
770  }
771  file2.close();
772  for (int i = 0; i < 120; i++)
773  for (int j = 0; j < 62; j++)
774  for (int k = 0; k < 51; k++)
775  spmu[i][j][k] = spmu[i + 1][j][k];
776  spmu[1][1][0] = 0.000853544;
777  //std::cout << "Set spmu to some value..." << std::endl;
778 
779  std::ostringstream File3LocStream;
780  File3LocStream << fInputDir << fInputFile3;
781  std::string File3Loc = File3LocStream.str();
782  cet::search_path sp3("FW_SEARCH_PATH");
783  if (sp3.find_file(fInputFile3, fROOTfile)) File3Loc = fROOTfile;
784  std::ifstream file3(File3Loc.c_str(), std::ios::in);
785  if (!file3.good())
786  throw cet::exception("MUSUNGen")
787  << "\nFile3 " << fInputFile3 << " not found in FW_SEARCH_PATH or at " << fInputDir
788  << "\n\n";
789 
790  lineNumber = index = 0;
791  while (file3.good()) {
792  //std::cout << "Looking at file 3...." << std::endl;
793  file3.getline(inputLine, 9999);
794  char* token;
795  token = strtok(inputLine, " ");
796  while (token != NULL) {
797  //std::cout << "While loop file 3..." << std::endl;
798  depth[index][lineNumber] = atof(token);
799  token = strtok(NULL, " ");
800  index++;
801  if (index == 360) {
802  //std::cout << "If statement file 3..." << std::endl;
803  index = 0;
804  lineNumber++;
805  }
806  }
807  }
808  file3.close();
809 
810  //
811  // Set up variables
812  //
813 
814  the1 = theta1;
815  the2 = theta2;
816  // for c2: c1 and c2 are unused
817  //double c1 = cos(M_PI/180.*theta1);
818  //double c2 = cos(M_PI/180.*theta2);
819  ph1 = M_PI / 180. * phi1;
820  ph2 = M_PI / 180. * phi2;
821  // for c2: dph is unused
822  //double dph = ph2-ph1;
823 
824  int ipc = 1;
825  double theta = theta1;
826  double dc = 1.;
827  double sc = 0.;
828  int iteration = 0;
829  while (theta < theta2 - dc / 2.) {
830  theta += dc / 2.;
831  double theta0 = M_PI / 180. * theta;
832  double cc = cos(theta0);
833  double ash = s_hor * cc;
834  double asv01 = s_ver1 * sqrt(1. - cc * cc);
835  double asv02 = s_ver2 * sqrt(1. - cc * cc);
836  int ic1 = (theta + 0.999);
837  int ic2 = ic1 + 1;
838  if (ic2 > 91) ic2 = 91;
839  if (ic1 < 1) ic1 = 1;
840  double phi = phi1;
841  double dp = 1.;
842 
843  while (phi < phi2 - dp / 2.) {
844  phi += dp / 2.;
845  // the long side of the cavern is pointing at 14 deg to the north:
846  // double phi0 = M_PI / 180. * (phi + 90. - 14);
847 
848  // Want our co-ord system going from East to South.
849  double phi0 = M_PI / 180. * (phi + fCavernAngle);
850 
851  double asv1 = asv01 * fabs(cos(phi0));
852  double asv2 = asv02 * fabs(sin(phi0));
853  double asv0 = ash + asv1 + asv2;
854  double fl = 1.;
855  if (figflag == 1) fl = asv0;
856  int ip1 = (phi + 0.999);
857  int ip2 = ip1 + 1;
858  if (ip2 > 360) ip2 = 1;
859  if (ip1 < 1) ip1 = 360;
860  double sp1 = 0.;
861 
862  for (int ii = 0; ii < 4; ii++) {
863  int iic = ii / 2;
864  int iip = ii % 2;
865  if (ip1 == 360 && (ii == 1 || ii == 3)) iip = -359;
866  if (fmu[ip1 + iip - 1][ic1 + iic - 1] < 0) {
867  if (pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 > 1e-6) {
868  std::cout << "Looking at fmu [ " << ip1 << " + " << iip << " - 1 (" << ip1 + iip - 1
869  << ") ] [ " << ic1 << " + " << iic << " - 1 (" << ic1 + iic - 1 << ") ] ."
870  << "\nChanging sp1 from " << sp1 << " to "
871  << sp1 + pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 << "..........."
872  << sp1 << " + 10 ^ (" << fmu[ip1 + iip - 1][ic1 + iic - 1] << ") / 4 "
873  << std::endl;
874  }
875  sp1 = sp1 + pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4;
876  }
877  }
878  /*
879  std::cout << iteration<< " time of new sc value! Theta " << theta << ", phi " << phi + dp / 2. << ", sc = " << sc + sp1 * fl * dp * M_PI / 180. * sin(theta0) * dc * M_PI / 180. << " = "
880  << sc << " + " << sp1 << " * " << fl << " * " << dp << " * " << M_PI/180 << " * sin(" << theta0 << ") * " << dc << " * " << M_PI/180 << ".....sin(theta)=" << sin(theta) << "\n"
881  << std::endl; */
882  sc = sc + sp1 * fl * dp * M_PI / 180. * sin(theta0) * dc * M_PI / 180.;
883  ++iteration;
884  ipc = ipc + 1;
885  fnmu[ipc - 1] = sc;
886  phi = phi + dp / 2.;
887  }
888 
889  theta = theta + dc / 2.;
890  }
891  //std::cout << *FI << " = " << sc << std::endl;
892  FI = sc;
893  for (int ipc1 = 0; ipc1 < ipc; ipc1++)
894  fnmu[ipc1] = fnmu[ipc1] / fnmu[ipc - 1];
895  }
896 
897  ////////////////////////////////////////////////////////////////////////////////
898  // sampling
899  ////////////////////////////////////////////////////////////////////////////////
900  void
901  MUSUN::sampling(double& E,
902  double& theta,
903  double& phi,
904  double& dep,
905  CLHEP::HepRandomEngine& engine)
906  {
907  CLHEP::RandFlat flat(engine);
908  CLHEP::RandGaussQ gauss(engine);
909 
910 #if 0 // this code is disabled for good
911  double xfl = flat.fire();
912  int loIndex = 0, hiIndex = 32400;
913  int i = (loIndex+hiIndex)/2;
914  bool foundIndex = false;
915  if( xfl < fnmu[loIndex] ) {
916  i = loIndex;
917  foundIndex = true;
918  } else if ( xfl > fnmu[hiIndex] ) {
919  i = hiIndex;
920  foundIndex = true;
921  } else if ( xfl > fnmu[i-1] && xfl <= fnmu[i] )
922  foundIndex = true;
923  while( !foundIndex ) {
924  if( xfl < fnmu[i] )
925  hiIndex = i;
926  else
927  loIndex = i;
928  i = (loIndex + hiIndex)/2;
929 
930  if( xfl > fnmu[i-1] && xfl <= fnmu[i] )
931  foundIndex = true;
932  }
933 #else
934  double xfl = flat.fire();
935  int i = 0;
936  while (xfl > fnmu[i])
937  ++i;
938 #endif
939  int ic = (i - 2) / 360;
940  int ip = i - 2 - ic * 360;
941 
942  xfl = flat.fire();
943  theta = the1 + ((double)ic + xfl);
944  xfl = flat.fire();
945  phi = ph1 + ((double)ip + xfl);
946  if (phi > 360) phi = phi - 360;
947  dep = depth[ip][ic] * fRockDensity;
948 
949  int ic1 = cos(M_PI / 180. * theta) * 50.;
950  if (ic1 < 0) ic1 = 0;
951  if (ic1 > 50) ic1 = 50;
952  int ip1 = dep / 200. - 16;
953  if (ip1 < 0) ip1 = 0;
954  if (ip1 > 61) ip1 = 61;
955 
956  xfl = flat.fire();
957 #if 0
958  loIndex = 0, hiIndex = 120;
959  int j = (loIndex+hiIndex)/2;
960  foundIndex = false;
961  if( xfl < spmu[loIndex][ip1][ic1] ) {
962  j = loIndex;
963  foundIndex = true;
964  } else if ( xfl > spmu[hiIndex][ip1][ic1] ) {
965  j = hiIndex;
966  foundIndex = true;
967  } else if ( xfl > spmu[j-1][ip1][ic1] && xfl <= spmu[j][ip1][ic1] )
968  foundIndex = true;
969  while( !foundIndex ) {
970  if( xfl < spmu[j][ip1][ic1] )
971  hiIndex = j;
972  else
973  loIndex = j;
974  j = (loIndex + hiIndex)/2;
975 
976  if( xfl > spmu[j-1][ip1][ic1] && xfl <= spmu[j][ip1][ic1] )
977  foundIndex = true;
978  }
979 #else
980  int j = 0;
981  while (xfl > spmu[j][ip1][ic1])
982  ++j;
983 #endif
984 
985  double En1 = 0.05 * (j - 1);
986  double En2 = 0.05 * (j);
987  E = pow(10., En1 + (En2 - En1) * flat.fire());
988 
989  return;
990  }
991 
992 } //end namespace evgen
993 
994 namespace evgen {
995 
996  DEFINE_ART_MODULE(MUSUN)
997 
998 } //end namespace evgen
MUSUN(fhicl::ParameterSet const &pset)
unsigned int NEvents
double fEmax
Maximum Kinetic Energy (GeV)
void beginRun(art::Run &run)
createEngine fInputDir
double fZmax
Maximum Z position.
double fYmax
Maximum Y position.
createEngine fYmax
double fCavernAngle
Angle of the detector from the North to the East.
createEngine fT0
double fEmin
Minimum Kinetic Energy (GeV)
createEngine fTDist
double fThetamin
Minimum theta.
std::string fInputDir
Input Directory.
createEngine fEmax
createEngine fEmin
process_name E
createEngine fPhimax
static const int kGAUS
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
double fT0
Central t position (ns) in world coordinates.
createEngine fPDG
process_name gaushit a
int fTDist
How to distribute t (gaus, or uniform)
double fRockDensity
createEngine fInputFile2
createEngine fXmin
int fPDG
PDG code of particles to generate.
void endRun(art::Run &run)
double fmu[360][91]
createEngine fXmax
createEngine fSigmaT
double fPhimax
Maximum phi.
std::string fInputFile3
Input File 3.
double fPhimin
Minimum phi.
double fChargeRatio
Charge ratio of particle / anti-particle.
fEngine(art::ServiceHandle< rndm::NuRandomService >() ->createEngine(*this, pset,"Seed"))
std::string fInputFile2
Input File 2.
createEngine fYmin
double depth[360][91]
int figflag
If want sampled from sphere or parallelepiped.
double spmu[121][62][51]
createEngine fThetamax
module to produce single or multiple specified particles in the detector
createEngine fThetamin
void initialization(double theta1, double theta2, double phi1, double phi2, int figflag, double s_hor, double s_ver1, double s_ver2, double &FI)
void SampleOne(unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
double fYmin
Minimum Y position.
double fZmin
Minimum Z position.
if &&[-z"$BASH_VERSION"] then echo Attempting to switch to bash bash shellSwitch exit fi &&["$1"= 'shellSwitch'] shift declare a IncludeDirectives for Dir in
createEngine fInputFile3
double fSigmaT
Variation in t position (ns)
createEngine fInputFile1
double fnmu[32401]
createEngine fZmin
createEngine fPhimin
do i e
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
double fThetamax
Maximum theta.
double fXmin
Minimum X position.
art::ServiceHandle< art::TFileService > tfs
createEngine fChargeRatio
void produce(art::Event &evt)
createEngine fZmax
TCEvent evt
Definition: DataStructs.cxx:8
pdgs k
Definition: selectors.fcl:22
void sampling(double &E, double &theta, double &phi, double &dep, CLHEP::HepRandomEngine &engine)
createEngine fRockDensity
std::string fInputFile1
Input File 1.
createEngine figflag
art framework interface to geometry description
BEGIN_PROLOG could also be cout
void beginJob()
double fXmax
Maximum X position.
createEngine fCavernAngle