Nickel source implementation, including bug fix and example mac file

GNUmakefile_root

@@ -65,5 +65,4 @@ $(G4TMPDIR)/%.o: %.cxx
 
 clean :
 	@rm -f $(G4TMPDIR)/*.o
-	@rm -f ./WCSimRootDict*
 

TemaplateNiSpectrum.dat

@@ -0,0 +1,11 @@
+Number of gamma modes:
+4
+Cumulative branching ratios:
+0.1 0.5 0.7 1.0
+Gamma multiplicities:
+1 2 2 3
+Gamma energies (MeV):
+9.000 0.000 0.000 0.
+8.000 0.500 0.000 0.
+8.000 0.300 0.000 0.
+7.000 1.000 1.000 0.

include/WCSimGenerator_NiBall.hh

@@ -0,0 +1,50 @@
+#ifndef WCSimGenerator_NiBall_hh
+#define WCSimGenerator_NiBall_hh 1
+
+#include "G4ReactionProductVector.hh"
+#include "G4ThreeVector.hh"
+#include <vector>
+#include "WCSimDetectorConstruction.hh"
+#include "TF2.h"
+#include "TGraph.h"
+#include "TMath.h"
+
+
+using namespace std;
+
+class WCSimDetectorConstruction;
+
+class WCSimGenerator_NiBall
+{
+	public:
+		WCSimGenerator_NiBall(WCSimDetectorConstruction* myDC);
+		~WCSimGenerator_NiBall();
+		void Initialize(G4String file);
+
+		G4ThreeVector GetRandomVertex(G4int tSymNumber);
+		int GetNiGammaMode()			{ return fNiGammaMode; }
+		int GetNiGammaMultiplicity() 	 	{ return fNiGammaMultiplicity; }
+		double * GetNiGammaEnergy()		{ return fNiGammaEnergy; }
+		double * GetNiGammaPosition() 		{ return fNiGammaPosition;}
+
+		void SettingNiBall(double fNiBallPosition[3], double rn[6]);
+
+	private:
+
+		int fNiGammaMode;
+		int fNiGammaMultiplicity;
+	//        char fNiBallSpectrum[200];
+	        std::string fNiBallSpectrum;
+		double fNiGammaEnergy[4];
+		double fNiGammaPosition[3];
+		WCSimDetectorConstruction*      myDetector;
+
+		// Constants of gamma production as taken from SK
+//		static G4double NiBallBR[37];
+//		static G4int NiBallMulti[37];
+//		static G4double NiBallEnergy[4][37];
+		static std::vector<G4double> NiBallBR;
+		static std::vector<G4int> NiBallMulti;
+		static std::vector<std::vector<G4double>> NiBallEnergy;
+};
+#endif

include/WCSimPrimaryGeneratorAction.hh

@@ -10,6 +10,7 @@
 #include <fstream>
 
 #include "WCSimRootOptions.hh"
+#include "WCSimGenerator_NiBall.hh"
 #include "WCSimGenerator_Radioactivity.hh"
 
 #include "TFile.h"
@@ -99,15 +100,28 @@ class WCSimPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
         G4int fRnSymmetry;
 
         G4bool   usePoissonPMT;
+	G4bool   useNiBallEvt; // Pablo: Ni ball calibration
         G4double poissonPMTMean;
 
+        // Variables for Ni ball
+        WCSimGenerator_NiBall* myNiBallGenerator;
+        double fNiBallPosition[3];
+	      std::string fNiBallSpectrum;
+        double fNiGammaDirection[3];
+        double rn[4];
+        G4double niball_X;
+        G4double niball_Y;
+        G4double niball_Z;
+        G4String niball_spectrum;
+
         // For injector events
         G4int nPhotons;
         G4int injectorOnIdx;
         G4double twindow;
         G4double openangle;
         G4double wavelength;
 
+
         // These go with jhfNtuple
         G4int mode;
         //InteractionType_t mode;
@@ -198,9 +212,21 @@ class WCSimPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
         inline void SetRadonSymmetry(G4int choice) 		{ fRnSymmetry = choice; }
         inline G4int GetRadonSymmetry() 			{ return fRnSymmetry; }
 
+        inline void SetNiBallX(G4double choice) { niball_X = choice; }
+        inline G4double GetNiBallX()            { return niball_X; }
+        inline void SetNiBallY(G4double choice) { niball_Y = choice; }
+        inline G4double GetNiBallY()            { return niball_Y; }
+        inline void SetNiBallZ(G4double choice) { niball_Z = choice; }
+        inline G4double GetNiBallZ()            { return niball_Z; }
+        inline void SetNiBallSpectrum(G4String choice) { niball_spectrum = choice; }
+        inline G4String GetNiBallSpectrum()            { return niball_spectrum; }
+
         inline void SetPoissonPMT(G4bool choice) { usePoissonPMT = choice; }
         inline G4bool IsUsingPoissonPMT(){ return usePoissonPMT; }
 
+        inline void SetNiBallEvtGenerator(G4bool choice)         { useNiBallEvt = choice; }
+        inline G4bool IsUsingNiBallEvtGenerator()                { return useNiBallEvt; }
+
         inline void SetPoissonPMTMean(G4double val){ poissonPMTMean = val; }
         inline G4double GetPoissonPMTMean(){ return poissonPMTMean; }
 };

include/WCSimPrimaryGeneratorMessenger.hh

@@ -36,6 +36,10 @@ class WCSimPrimaryGeneratorMessenger: public G4UImessenger
   G4UIcmdWithAString* fileNameCmd;
   G4UIcmdWithABool*   poisCmd;
   G4UIcmdWithADouble* poisMeanCmd;
+  G4UIcmdWithADouble* niball_x_Cmd;
+  G4UIcmdWithADouble* niball_y_Cmd;
+  G4UIcmdWithADouble* niball_z_Cmd;
+  G4UIcmdWithAString* niball_spectrumCmd;
   G4UIcmdWithAString* isotopeCmd;
   G4UIcmdWithAString* radonScalingCmd;
   G4UIcmdWithADouble* radioactive_time_window_Cmd;

macros/mPMT_NiBall.mac

@@ -0,0 +1,253 @@
+# Sample setup macro with no visualization
+
+/run/verbose 1
+/tracking/verbose 0
+/hits/verbose 0
+
+## select the geometry
+# Default config if you do nothing is currently SuperK
+#
+
+#/WCSim/WCgeom nuPRISM
+
+#Select which PMT to use
+#/WCSim/nuPRISM/SetPMTType PMT8inch
+#/WCSim/nuPRISM/SetPMTPercentCoverage 40
+#Set height of nuPRISM inner detector
+#/WCSim/nuPRISM/SetDetectorHeight 6. m
+#Set vertical position of inner detector, in beam coordinates
+#/WCSim/nuPRISM/SetDetectorVerticalPosition 0. m
+#Set diameter of inner detector
+#/WCSim/nuPRISM/SetDetectorDiameter 8. m
+# Set Gadolinium doping (concentration is in percent)
+#/WCSim/DopingConcentration 0.1
+#/WCSim/DopedWater false
+#/WCSim/Construct
+#/WCSim/mPMT/PMTtype_inner PMT3inchR14374
+#/WCSim/mPMT/PMTtype_outer PMT3inchR14374
+
+#Use mPMTs settings (uncomment/delete the above)
+#/WCSim/WCgeom HyperK_HybridmPMT10PC
+#/WCSim/mPMT/PMTtype_inner PMT3inchR14374
+#/WCSim/mPMT/PMTtype_outer PMT3inchR14374
+#/WCSim/WCgeom HyperK_mPMT
+#/WCSim/WCgeom HyperK
+#/WCSim/WCgeom nuPRISM_mPMT
+#/WCSim/WCgeom nuPRISMShort_mPMT
+# Set Gadolinium doping (concentration is in percent)
+#/WCSim/DopingConcentration 0.1
+#/WCSim/DopedWater false
+/WCSim/Construct
+## OR for single mPMT mode or updating mPMT parameters:
+#/control/execute macros/mPMT_nuPrism1.mac         ## mPMT options: mPMT_nuPrism1.mac and 2.mac
+
+/WCSim/tuning/rayff 0.75
+/WCSim/tuning/bsrff 2.25
+/WCSim/tuning/abwff 1.30
+/WCSim/tuning/rgcff 0.32
+
+#/WCSim/mPMT/reflectorHeightID 7. mm
+#/WCSim/mPMT/reflectorHeightID 7. mm
+#/WCSim/mPMT/reflectorAngleID 48. deg
+
+#Added for the PMT QE option 08/17/10 (XQ)
+# 1. Stacking only mean when the photon is generated
+# the QE is applied to reduce the total number of photons
+# 2. Stacking and sensitivity detector
+# In the stacking part, the maximum QE is applied to reduce
+# the total number of photons
+# On the detector side, the rest of QE are applied according to QE/QE_max
+# distribution. This option is in particular important for the WLS
+# 3. The last option means all the QE are applied at the detector
+# Good for the low energy running.
+# 4. Switch off the QE, ie. set it at 100%
+/WCSim/PMTQEMethod     Stacking_Only
+#/WCSim/PMTQEMethod     Stacking_And_SensitiveDetector
+#/WCSim/PMTQEMethod     SensitiveDetector_Only
+#/WCSim/PMTQEMethod     DoNotApplyQE
+
+#turn on or off the collection efficiency
+/WCSim/PMTCollEff on
+
+# command to choose save or not save the pi0 info 07/03/10 (XQ)
+/WCSim/SavePi0 false
+
+#choose the Trigger & Digitizer type (and options)
+/DAQ/Digitizer SKI
+/DAQ/Trigger NDigits
+/DAQ/TriggerNDigits/Threshold 0
+/DAQ/TriggerNDigits/AdjustForNoise false
+
+#grab the other DAQ options (thresholds, timing windows, etc.)
+/control/execute macros/daq.mac
+
+
+# default dark noise frequency (and conversion factor) is PMT property (NEW), set in the code.
+# Below gives possibility to overwrite nominal values, eg. to switch OFF the Dark Noise.
+#/DarkRate/SetDarkRate 0 kHz   #Turn dark noise off
+#/DarkRate/SetDarkRate 4.2 kHz #This is the value for SKI set in SKDETSIM.
+#/DarkRate/SetDarkRate 8.4 kHz #For 20 inch HPDs and Box and Line PMTs, based on High QE 20in R3600 dark rate from EGADS nov 2014
+#/DarkRate/SetDarkRate 3.0 kHz #For 12 inch HPDs and Box and Line PMTs, based on High QE 20in R3600 dark rate from EGADS nov 2014
+
+# command to multiply the dark rate.
+# Convert dark noise frequency before digitization to after digitization by setting suitable factor
+# Again, this is now a PMT property and can be overridden here
+#/DarkRate/SetConvert 1.367  #For Normal PMT
+#/DarkRate/SetConvert 1.119 #For HPDs
+#/DarkRate/SetConvert 1.126 #For Box and Line PMTs
+
+# Select which time window(s) to add dark noise to
+#/DarkRate/SetDarkMode 0 to add dark noise to a time window starting at
+#/DarkRate/SetDarkLow to /DarkRate/SetDarkHigh [time in ns]
+#/DarkRate/SetDarkMode 1 adds dark noise hits to a window of
+#width /DarkRate/SetDarkWindow [time in ns] around each hit
+#i.e. hit time ± (/DarkRate/SetDarkWindow) / 2
+/DarkRate/SetDarkMode 1
+/DarkRate/SetDarkHigh 100000
+/DarkRate/SetDarkLow 0
+/DarkRate/SetDarkWindow 4000
+
+#Uncomment one of the lines below if you want to use the OGLSX or RayTracer visualizer
+#/control/execute macros/visOGLSX.mac
+#/control/execute macros/visRayTracer.mac
+#/control/execute macros/visOGLQT.mac             ## NEW
+
+## select the input nuance-formatted vector file
+## you can of course use your own
+#/mygen/generator muline
+#/mygen/vecfile inputvectorfile
+#/mygen/vecfile h2o.2km.001-009x3_G4.kin
+#/mygen/vecfile mu+.out
+
+# Or you can use the G4 Particle Gun
+# for a full list of /gun/ commands see:
+# http://geant4.web.cern.ch/geant4/G4UsersDocuments/UsersGuides/ForApplicationDeveloper/html/Control/UIcommands/_gun_.html
+#/mygen/generator gun
+#/gun/particle gamma
+#/gun/particle e-
+#/gun/particle mu-
+#/gun/particle pi0
+#/gun/energy 500 MeV
+#/gun/direction 1 0 0
+#/gun/position 0 0 0 m
+#/gun/number 100
+
+# Or you can use the G4 General Particle Source
+# you can do a lot more with this than a monoenergetic, monodirectional, single-particle gun
+# for a full list of /gps/ commands see:
+# https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch02s07.html
+/mygen/generator gps
+#/gps/particle gamma
+/gps/particle e-
+/gps/energy 10 MeV
+#/gps/energy 1 MeV
+/gps/direction 1 0 0
+#/gps/position 0. 0. 0. m
+#/gps/position 3.5 0. 0. m
+
+# Position: a surface
+#/gps/pos/type 
+#/gps/pos/shape Cylinder
+#/gps/pos/halfz 4.5 m
+#/gps/pos/radius 3.9 m
+
+# Position: a cylinder
+/gps/pos/type Volume
+/gps/pos/shape Cylinder
+/gps/pos/halfz 30.4 m
+/gps/pos/radius 38.4 m
+
+# Position: a sphere
+#/gps/pos/type Volume
+#/gps/pos/shape Sphere
+#/gps/pos/radius 200. cm
+#/gps/pos/centre 0. 0. 0. m
+
+# Direction: unidirection along x
+#/gps/ang/type beam1d
+#/gps/ang/rot1 0 0 1
+#/gps/ang/rot2 0 1 0
+
+# Direction: a sphere
+/gps/ang/type iso
+#/gps/ang/mintheta 0 deg
+#/gps/ang/maxtheta 180 deg
+#/gps/ang/minphi 0 deg
+#/gps/ang/maxphi 360 deg
+
+# spectrum: uniform between 200 and 1000 MeV
+#/gps/ene/type Lin
+#/gps/ene/gradient 0
+#/gps/ene/intercept 1
+#/gps/ene/type Mono
+#/gps/ene/mono 500 MeV
+#/gps/ene/min 200 MeV
+#/gps/ene/max 1000 MeV
+
+#/random/setSeeds 1225 1500
+
+# Or you can use the laser option
+# This is equivalent to the gps command, except that the gps particle energies are saved ignoring their mass
+# for a full list of /gps/ commands see:
+# https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch02s07.html
+# It is used for laser calibration simulation
+#/mygen/generator laser
+#/gps/particle opticalphoton
+#/gps/energy 3.5 keV
+#/gps/direction 1 0 0
+#/gps/position 3.5 0 0 m
+#/gps/number 1000
+#/gps/ang/type iso
+#/gps/ang/mintheta 0 deg
+#/gps/ang/maxtheta 0 deg
+#/gps/ang/minphi 0 deg
+#/gps/ang/maxphi 90 deg
+
+##### NEW
+/Tracking/fractionOpticalPhotonsToDraw 0.0
+
+# Radon option
+# IMPORTANT: You need to compile WCSim with cmake in order to use the radioactive decays
+# Generate Bi214 decay in the water following a Model extracted from SK-IV solar data
+#/mygen/generator radon
+# This window should match the one in biasprofile.dat if time bias is used (see below)
+#/mygen/radioactive_time_window 1.e-4
+# Needed if time bias is not used, in order to keep timing below one sec.
+#/DAQ/RelativeHitTime true
+# The following set which scenario is used to scale SK model to HK
+#/mygen/radon_scaling A # 	Relative scaling with respect to full ID volume (Pessimistic)
+#/mygen/radon_scaling B #	Absolute scaling with respect to ID border (Optimistic)
+#/mygen/radon_scaling 0 #	Uniform distribution
+#/mygen/radon_symmetry 1 #	Divide the detector in x slices using the theta symmetry
+#Clear the source and tell GPS to make multiple vertices
+#/gps/source/clear
+
+# Needed to prevent other decay than Bi-214
+#/grdm/nucleusLimits 214 214 83 83 
+# if you need to use time bias
+#/grdm/decayBiasProfile biasprofile.dat
+#/grdm/analogueMC 0
+#/grdm/BRbias 0
+
+# Nickel calibration
+# For the time being set the coordinates for the center position of the ball in cm.
+# In future modify to select the calibration port and the Z-position of the ball.
+/mygen/generator niball
+/mygen/niball_x 10.0
+/mygen/niball_y -300.0
+/mygen/niball_z 1.0
+/mygen/niball_spectrum TemaplateNiSpectrum.dat
+
+#Clear the source and tell GPS to make multiple vertices
+/gps/source/clear
+
+
+## change the name of the output root file, default = wcsim.root
+/WCSimIO/RootFile wcsim_output.root
+
+## Boolean to select whether to save the NEUT RooTracker vertices in the output file, provided you used
+## a NEUT vector file as input
+/WCSimIO/SaveRooTracker 0
+
+/run/beamOn 10
+#exit