diff --git a/CMakeLists.txt b/CMakeLists.txt
index ef68eed77aae64caa355c6b73efbe0d64bbb0f4f..1c3b76900a8e35eeae486219606a6752759c9e27 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -57,6 +57,8 @@ TARGET_LINK_LIBRARIES(plotHisto ${ROOT_LIBRARIES})
set(lightGuide_SCRIPTS
vis.mac
run1.mac
+ beam.mac
+ geom.mac
)
foreach(_script ${lightGuide_SCRIPTS})
diff --git a/beam.mac b/beam.mac
new file mode 100644
index 0000000000000000000000000000000000000000..638b8d57d450ec96da73d82837730f8b4fb317e9
--- /dev/null
+++ b/beam.mac
@@ -0,0 +1,41 @@
+#----- Set the beam energy profile -----
+/gps/particle opticalphoton
+/gps/ene/type Gauss
+/gps/ene/mono 3.4 eV
+/gps/ene/sigma 1 eV
+/gps/polarization 0 0 0
+
+
+
+#----- Set the beam geometry -----
+# Circular beam source. Simulates a single fiber optic
+/gps/pos/type Beam
+/gps/pos/shape Circle
+#/gps/pos/radius 0.75 mm
+/gps/pos/radius 102 mm
+#/gps/pos/sigma_r 0.002 mm
+
+# Rectangular plane for infinite resolution
+#/gps/pos/type Plane
+#/gps/pos/shape Rectangle
+#/gps/pos/halfx 40 mm
+#/gps/pos/halfy 68 mm
+
+# the incident surface is in the x-z plane just below the gas
+/gps/pos/rot1 1 0 0
+/gps/pos/rot2 0 0 1
+/gps/pos/centre 0. -0.1 0. mm
+
+# the photons are emitted around the y-axis
+/gps/ang/rot1 1 0 0
+/gps/ang/rot2 0 0 1
+/gps/ang/type beam1d
+
+#----- Set beam angular dispersion from histogram -----
+/gps/ang/surfnorm false
+/gps/ang/type user
+/control/execute hist.mac
+
+#----- Set beam angular dispersion for fiber optic -----
+# N.A. is 0.22 so acceptance cone is ~25.4 degrees
+#/gps/ang/sigma_r 25.40 deg
diff --git a/geom.mac b/geom.mac
new file mode 100644
index 0000000000000000000000000000000000000000..3dcb403ccde192e59f5917bd8a9d90de4c102806
--- /dev/null
+++ b/geom.mac
@@ -0,0 +1,58 @@
+#----- Set the light guide to be used and position it -----
+#If no model is selected, the program will default to a
+#simplified version of the 2018 testbeam light guide
+
+#/lightGuide/model/CADmodel ../models/LightGuide2018TB.stl
+#/lightGuide/model/rotate 0 90 0
+#/lightGuide/model/translate 0 -250 200 mm
+#/lightGuide/model/translatePMT 0 339 0 mm
+
+#/lightGuide/model/CADmodel ../models/Winston_cone.stl
+#/lightGuide/model/rotate 0 0 90
+#/lightGuide/model/translate 0 -111 0 mm
+#/lightGuide/model/translatePMT 0 339 0 mm
+
+#/lightGuide/model/CADmodel ../models/LightGuide2007BigPMT.stl
+#/lightGuide/model/rotate 0 -90 0
+#/lightGuide/model/translate 0 -250 -200 mm
+#/lightGuide/model/translatePMT 0 131 0 mm
+
+#/lightGuide/model/CADmodel ../models/LightGuide2007SmallPMT.stl
+#/lightGuide/model/rotate 0 -90 0
+#/lightGuide/model/translate 0 -250 -200 mm
+#/lightGuide/model/translatePMT 0 131 0 mm
+
+#/lightGuide/model/CADmodel ../models/run2.stl
+#/lightGuide/model/rotate 0 0 0
+#/lightGuide/model/translate 200 -250 0 mm
+#/lightGuide/model/translatePMT 0 -119 0 mm
+
+#/lightGuide/model/CADmodel ../models/run3wc.stl
+#/lightGuide/model/rotate 0 0 90
+#/lightGuide/model/translate 0 -319 0 mm
+#/lightGuide/model/translatePMT 0 -119 0 mm
+
+/lightGuide/model/CADmodel ../models/YuvalCone.stl
+/lightGuide/model/translate 200 -136.75 0 mm
+/lightGuide/model/translatePMT 0 -136.5 0 mm
+/lightGuide/model/PMTDiameter 170
+
+#/lightGuide/model/CADmodel ../models/YuvalHad.stl
+#/lightGuide/model/translate 0 -319 0 mm
+#/lightGuide/model/translatePMT 0 -136.5 0 mm
+#/lightGuide/model/PMTDiameter 170
+
+
+#----- Set surface properties of the light guide -----
+/lightGuide/surface/Model unified
+/lightGuide/surface/Type dielectric_metal
+/lightGuide/surface/Finish ground
+#/lightGuide/surface/SigmaAlpha .1
+#/lightGuide/surface/Property SPECULARLOBECONSTANT 0.000002 .4 0.000008 .4
+#/lightGuide/surface/Property SPECULARSPIKECONSTANT 0.000002 .1 0.000008 .1
+#/lightGuide/surface/Property BACKSCATTERCONSTANT 0.000002 .05 0.000008 .05
+#/lightGuide/surface/Property REFLECTIVITY 0.000002 .89 0.000008 .89
+
+
+#----- Set the refractive index of the gas the light guide is inside of -----
+#/lightGuide/gasProperty RINDEX 0.000002 1.0 0.000008 1.0
diff --git a/include/ASCIIPrimaryGenerator.hh b/include/ASCIIPrimaryGenerator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c4403ae4634f1516655660297ed14153869e4b87
--- /dev/null
+++ b/include/ASCIIPrimaryGenerator.hh
@@ -0,0 +1,70 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file ASCIIPrimaryGenerator.hh
+/// \brief Definition of the ASCIIPrimaryGenerator class
+//
+//
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#ifndef ASCIIPrimaryGenerator_h
+#define ASCIIPrimaryGenerator_h 1
+
+#include "G4VPrimaryGenerator.hh"
+#include "lgAnalysis.hh"
+
+#include <fstream>
+#include <vector>
+
+class G4Event;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class ASCIIPrimaryGenerator : public G4VPrimaryGenerator
+{
+ public:
+ ASCIIPrimaryGenerator();
+ ~ASCIIPrimaryGenerator();
+
+ virtual void SetInputFile(G4String _name);
+ virtual void GetNextEvent( );
+
+ public:
+ virtual void GeneratePrimaryVertex(G4Event*);
+
+ private:
+ std::vector< G4ThreeVector >* fPositionVec;
+ std::vector< G4ThreeVector >* fMomentumVec;
+ std::vector< G4double >* fEnergyVec;
+ G4int fEventNo;
+ G4int fnEvents;
+ std::ifstream fInputFile;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
diff --git a/include/EventAction.hh b/include/EventAction.hh
index cd2d9f38177a106a294afab6ca0cabd6850fe4f9..2f13f0f1c3e497a3a9a10ea17a722e3ef783dec9 100644
--- a/include/EventAction.hh
+++ b/include/EventAction.hh
@@ -29,6 +29,8 @@
#include "globals.hh"
#include "G4UserEventAction.hh"
+#include <vector>
+
class EventAction : public G4UserEventAction{
public:
@@ -36,11 +38,16 @@ class EventAction : public G4UserEventAction{
virtual ~EventAction();
virtual void BeginOfEventAction( const G4Event* event );
- virtual void EndOfEventAction( const G4Event* event );
+ virtual void EndOfEventAction( const G4Event* event );
+
+ inline std::vector< std::vector<double>* > GetVectors( ){return fPtrVec;}
private:
G4int hitsCollID;
G4int fEventNo;
+
+ std::vector<double> x,z,xHit,zHit,time,theta;
+ std::vector< std::vector<double>* > fPtrVec;
};
#endif
diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh
index 584a1c779e671d0f27ea4c9d5784bde34bb7796c..ca663a090eb39581c584264058fa728f010aad92 100644
--- a/include/PrimaryGeneratorAction.hh
+++ b/include/PrimaryGeneratorAction.hh
@@ -29,7 +29,11 @@
#define PrimaryGeneratorAction_h 1
#include "G4VUserPrimaryGeneratorAction.hh"
+#include "PrimaryGeneratorMessenger.hh"
+#include "ASCIIPrimaryGenerator.hh"
+
#include "globals.hh"
+#include "G4GeneralParticleSource.hh"
class G4GeneralParticleSource;
class G4Event;
@@ -43,9 +47,14 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
public:
virtual void GeneratePrimaries(G4Event*);
+ virtual void SetInputFile(G4String _name);
private:
- G4GeneralParticleSource* fParticleGun;
+ G4GeneralParticleSource* fParticleGun;
+ ASCIIPrimaryGenerator* fASCIIParticleGun;
+ PrimaryGeneratorMessenger* fMessenger;
+ G4bool fUseInput;
+
};
#endif /*PrimaryGeneratorAction_h*/
diff --git a/include/PrimaryGeneratorMessenger.hh b/include/PrimaryGeneratorMessenger.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4e5053811f1a97ac85cd4479b28a489964c007ae
--- /dev/null
+++ b/include/PrimaryGeneratorMessenger.hh
@@ -0,0 +1,66 @@
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file optical/OpNovice2/include/DetectorMessenger.hh
+/// \brief Definition of the DetectorMessenger class
+//
+//
+//
+
+#ifndef PrimaryGeneratorMessenger_h
+#define PrimaryGeneratorMessenger_h 1
+
+#include "globals.hh"
+#include "G4UImessenger.hh"
+
+class PrimaryGeneratorAction;
+class G4UIdirectory;
+class G4UIcommand;
+class G4UIcmdWithAString;
+class G4UIcmdWithAnInteger;
+class G4UIcmdWithADouble;
+class G4UIcmdWithADoubleAndUnit;
+class G4UIcmdWithoutParameter;
+class G4UIcmdWith3VectorAndUnit;
+
+
+class PrimaryGeneratorMessenger: public G4UImessenger{
+ public:
+
+ PrimaryGeneratorMessenger(PrimaryGeneratorAction* );
+ ~PrimaryGeneratorMessenger();
+
+ virtual void SetNewValue(G4UIcommand*, G4String);
+
+ private:
+
+ PrimaryGeneratorAction* fGenerator;
+
+ G4UIdirectory* fGeneratorDir;
+
+ G4UIcmdWithAString* fInputFileCmd;
+
+
+};
+#endif
diff --git a/run1.mac b/run1.mac
index bfb325ac4eace3f91e3a0c6b7921b183e7f23ac5..4675a885a6001ed97c85aa23149b08eeec4506be 100644
--- a/run1.mac
+++ b/run1.mac
@@ -1,8 +1,6 @@
-# Macro file for JZCaPA beam test 2018
+# Macro file for lightGuide
#
# To be run preferably in batch, without graphics:
-#
-# 31 TeV Pb ion
# Initialize kernel
/run/initialize
@@ -16,121 +14,23 @@
/hits/verbose 0
-
-#----- Set the light guide to be used and position it -----
-#If no model is selected, the program will default to a
-#simplified version of the 2018 testbeam light guide
-
-#/lightGuide/model/CADmodel ../models/LightGuide2018TB.stl
-#/lightGuide/model/rotate 0 90 0
-#/lightGuide/model/translate 0 -250 200 mm
-#/lightGuide/model/translatePMT 0 339 0 mm
-
-#/lightGuide/model/CADmodel ../models/Winston_cone.stl
-#/lightGuide/model/rotate 0 0 90
-#/lightGuide/model/translate 0 -111 0 mm
-#/lightGuide/model/translatePMT 0 339 0 mm
-
-#/lightGuide/model/CADmodel ../models/LightGuide2007BigPMT.stl
-#/lightGuide/model/rotate 0 -90 0
-#/lightGuide/model/translate 0 -250 -200 mm
-#/lightGuide/model/translatePMT 0 131 0 mm
-
-#/lightGuide/model/CADmodel ../models/LightGuide2007SmallPMT.stl
-#/lightGuide/model/rotate 0 -90 0
-#/lightGuide/model/translate 0 -250 -200 mm
-#/lightGuide/model/translatePMT 0 131 0 mm
-
-#/lightGuide/model/CADmodel ../models/run2.stl
-#/lightGuide/model/rotate 0 0 0
-#/lightGuide/model/translate 200 -250 0 mm
-#/lightGuide/model/translatePMT 0 -119 0 mm
-
-#/lightGuide/model/CADmodel ../models/run3wc.stl
-#/lightGuide/model/rotate 0 0 90
-#/lightGuide/model/translate 0 -319 0 mm
-#/lightGuide/model/translatePMT 0 -119 0 mm
-
-/lightGuide/model/CADmodel ../models/YuvalCone.stl
-/lightGuide/model/translate 200 -136.75 0 mm
-/lightGuide/model/translatePMT 0 -136.5 0 mm
-/lightGuide/model/PMTDiameter 170
-
-#/lightGuide/model/CADmodel ../models/YuvalHad.stl
-#/lightGuide/model/translate 0 -319 0 mm
-#/lightGuide/model/translatePMT 0 -136.5 0 mm
-#/lightGuide/model/PMTDiameter 170
-
-
-#----- Set surface properties of the light guide -----
-/lightGuide/surface/Model unified
-/lightGuide/surface/Type dielectric_metal
-/lightGuide/surface/Finish ground
-#/lightGuide/surface/SigmaAlpha .1
-#/lightGuide/surface/Property SPECULARLOBECONSTANT 0.000002 .4 0.000008 .4
-#/lightGuide/surface/Property SPECULARSPIKECONSTANT 0.000002 .1 0.000008 .1
-#/lightGuide/surface/Property BACKSCATTERCONSTANT 0.000002 .05 0.000008 .05
-#/lightGuide/surface/Property REFLECTIVITY 0.000002 .89 0.000008 .89
-
-
-
-#----- Set the refractive index of the gas the light guide is inside of -----
-#/lightGuide/gasProperty RINDEX 0.000002 1.0 0.000008 1.0
-
-
-
-#----- Set the beam energy profile -----
-/gps/particle opticalphoton
-/gps/ene/type Gauss
-/gps/ene/mono 3.4 eV
-/gps/ene/sigma 1 eV
-/gps/polarization 0 0 0
-
-
-
-#----- Set the beam geometry -----
-# Circular beam source. Simulates a single fiber optic
-/gps/pos/type Beam
-/gps/pos/shape Circle
-#/gps/pos/radius 0.75 mm
-/gps/pos/radius 102 mm
-#/gps/pos/sigma_r 0.002 mm
-
-# Rectangular plane for infinite resolution
-#/gps/pos/type Plane
-#/gps/pos/shape Rectangle
-#/gps/pos/halfx 40 mm
-#/gps/pos/halfy 68 mm
-
-# the incident surface is in the x-z plane just below the gas
-/gps/pos/rot1 1 0 0
-/gps/pos/rot2 0 0 1
-/gps/pos/centre 0. -0.1 0. mm
-
-# the photons are emitted around the y-axis
-/gps/ang/rot1 1 0 0
-/gps/ang/rot2 0 0 1
-/gps/ang/type beam1d
-
-#----- Set beam angular dispersion from histogram -----
-/gps/ang/surfnorm false
-/gps/ang/type user
-/control/execute hist.mac
-
-#----- Set beam angular dispersion for fiber optic -----
-# N.A. is 0.22 so acceptance cone is ~25.4 degrees
-#/gps/ang/sigma_r 25.40 deg
-
+#Lightguide selection and properties can be set with geom.mac
+/control/execute geom.mac
-#if using these commands in interactive mode to visualize, use this to allow more events to be displayed (do not use in batch mode)
-/vis/ogl/set/displayListLimit 10000
+#Uncomment to generate events from ASCII file, set the file path here.
+#make sure to set beamOn # to the number of events contained in the ASCII file (Have to figure out how to automatically handle this)
+/Input/FileName test.txt
+#Uncomment to use gps style of setting the beam profile
+#/control/execute beam.mac
+#set number of photons per event
+#/gps/number 10
##############################################################
##############################################################
# number of events
-/run/beamOn 10000000
+/run/beamOn 10
##############################################################
##############################################################
diff --git a/src/ASCIIPrimaryGenerator.cc b/src/ASCIIPrimaryGenerator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4311fcf5c2485d82dc69dca21e43456e57fc0904
--- /dev/null
+++ b/src/ASCIIPrimaryGenerator.cc
@@ -0,0 +1,177 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file ASCIIPrimaryGenerator.cc
+/// \brief Implementation of the PrimaryGenerator1 class
+//
+//
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "ASCIIPrimaryGenerator.hh"
+
+#include "G4Event.hh"
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4PrimaryParticle.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "Randomize.hh"
+
+/*
+ *
+ */
+ASCIIPrimaryGenerator::ASCIIPrimaryGenerator()
+: G4VPrimaryGenerator(),fEventNo(0)
+{
+ fPositionVec = new std::vector< G4ThreeVector >;
+ fMomentumVec = new std::vector< G4ThreeVector >;
+ fEnergyVec = new std::vector< G4double >;
+}
+
+
+/*
+ *
+ */
+ASCIIPrimaryGenerator::~ASCIIPrimaryGenerator()
+{ }
+
+
+/*
+ *
+ */
+void ASCIIPrimaryGenerator::GeneratePrimaryVertex(G4Event* event)
+{
+ G4ParticleDefinition* particleDefinition=
+ G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
+
+ std::vector< G4PrimaryVertex* > vertexVec;
+ std::vector< G4PrimaryParticle* > particleVec;
+
+ fPositionVec->clear();
+ fMomentumVec->clear();
+ fEnergyVec->clear();
+
+ GetNextEvent();
+
+ // Pass the original event number through to the ouput
+ auto analysisManager = G4AnalysisManager::Instance();
+ analysisManager->FillNtupleDColumn(6, fEventNo );
+
+ //For each photon read in, add a new particle and vertex to the event
+ int nPhotons = fPositionVec->size();
+ for(int i = 0; i < nPhotons; i++){
+ particleVec.push_back( new G4PrimaryParticle(particleDefinition) );
+ particleVec.back()->SetMomentumDirection( fMomentumVec->at(i) );
+ particleVec.back()->SetKineticEnergy( fEnergyVec->at(i)*eV );
+
+ vertexVec.push_back( new G4PrimaryVertex( fPositionVec->at(i) , 0*s ) );
+ vertexVec.back()->SetPrimary( particleVec.back() );
+ event->AddPrimaryVertex( vertexVec.back() );
+ }
+}
+
+
+/*
+ * Fill member vectors with position, momentum and energy of all
+ * photons in the next event stored in the input ascii file
+ */
+void ASCIIPrimaryGenerator::GetNextEvent(){
+ std::string str;
+ double x,y,z,px,py,pz,energy;
+
+ // The first line of each event should be "Event #"
+ // If it's not, there's an error in the file formatting
+ // or we've reached the end of the file
+ getline(fInputFile,str);
+ if(str[0] == 'E'){
+ fEventNo = std::atoi( str.substr( 6, str.size() ).c_str() );
+ }else{
+ G4cout << "Error reading input event. Searching for next event..." << G4endl;
+ while( !fInputFile.eof() ){
+ getline(fInputFile,str);
+ if(str[0] == 'E'){
+ fEventNo = std::atoi( str.substr( 6, str.size() ).c_str() );
+ break;
+ }
+ }
+ if( fInputFile.eof() ){
+ G4cout << "Reached end of input file" << G4endl;
+ fEventNo = -1;
+ return;
+ }
+ }
+
+ // Loop through the event and read all photons.
+ // Format is:
+ // (Vertex) V, x, y, z
+ // (Momentum) P, px, py, pz, energy
+ // Then push all of those to their respective member vectors
+ getline(fInputFile,str);
+ while( !fInputFile.eof() ){
+ //
+ //Get Vertex
+ if(str[0] != 'V'){
+ G4cout << "Error reading input vertex at position " << fInputFile.tellg() << G4endl;
+ return;
+ }
+ x = std::atof( str.substr(2,12).c_str() );
+ y = std::atof( str.substr(13,23).c_str() );
+ z = std::atof( str.substr(24,34).c_str() );
+ fPositionVec->push_back( G4ThreeVector(x*mm,y*mm,z*mm) );
+
+ //
+ //Get momentum and energy
+ getline(fInputFile,str);
+ if(str[0] != 'P'){
+ G4cout << "Error reading input momentum at position " << fInputFile.tellg() << G4endl;
+ return;
+ }
+ px = std::atof( str.substr(2,12).c_str() );
+ py = std::atof( str.substr(13,23).c_str() );
+ pz = std::atof( str.substr(24,34).c_str() );
+ fMomentumVec->push_back( G4ThreeVector(px,py,pz) );
+ energy = std::atof( str.substr(35,45).c_str() );
+ fEnergyVec->push_back( energy );
+
+ //
+ //The next line is either another Vertex, or the end of the event
+ getline(fInputFile,str);
+ if( str[0] == 'E' ) return;
+ }
+}
+
+/*
+ * Open the input file and get the total number of events from it
+ */
+void ASCIIPrimaryGenerator::SetInputFile(G4String _name){
+ std::string str;
+ fInputFile.open( _name.c_str() );
+
+ getline(fInputFile,str);
+ fnEvents = std::atoi( str.substr(15,str.size() ).c_str() );
+}
diff --git a/src/ActionInitialization.cc b/src/ActionInitialization.cc
index 04599319ab116fbfec2de04127d03a6bb2e53a30..06fada66d0c19f8f47240ab3c1c4dfad18aed6a7 100644
--- a/src/ActionInitialization.cc
+++ b/src/ActionInitialization.cc
@@ -57,7 +57,7 @@ void ActionInitialization::BuildForMaster() const{
*/
void ActionInitialization::Build() const{
SetUserAction(new PrimaryGeneratorAction());
- SetUserAction(new RunAction(ffileName));
SetUserAction(new SteppingAction());
SetUserAction(new EventAction());
+ SetUserAction(new RunAction(ffileName));
}
diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index af4d9ac1f67e0dc323f95bcbdf5373bd656499e7..afeb43eaf911682bc57d659d97b960a3a4724054 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -85,7 +85,7 @@ DetectorConstruction::DetectorConstruction()
*
*/
DetectorConstruction::~DetectorConstruction(){
-
+ delete m_DetectorMessenger;
}
/* Default geometry (2018 testbeam) is created here
diff --git a/src/EventAction.cc b/src/EventAction.cc
index 2581553e91ab6eaf99772355b28e8a2e7c606e00..d5277357f2eb9c7f6bbd1bae4be4de4cfffd5963 100644
--- a/src/EventAction.cc
+++ b/src/EventAction.cc
@@ -43,6 +43,12 @@
EventAction::EventAction()
: G4UserEventAction(),
hitsCollID(0){
+ fPtrVec.push_back(&x);
+ fPtrVec.push_back(&z);
+ fPtrVec.push_back(&xHit);
+ fPtrVec.push_back(&zHit);
+ fPtrVec.push_back(&time);
+ fPtrVec.push_back(&theta);
}
@@ -60,6 +66,10 @@ void EventAction::BeginOfEventAction(const G4Event* event){
fEventNo = event->GetEventID();
if(fEventNo%100000 == 0) G4cout << "Begin Event " << fEventNo << G4endl;
hitsCollID = 0;
+
+ for(uint i = 0; i < fPtrVec.size(); i++){
+ fPtrVec.at(i)->clear();
+ }
}
/*
@@ -71,16 +81,27 @@ void EventAction::EndOfEventAction(const G4Event* event){
if(HCE){
HitsCollection* HC = (HitsCollection*)(HCE->GetHC(hitsCollID));
if( HC->entries() ){
- PMTHit* aHit = (*HC)[0];
+ for(int i = 0; i < HC->entries(); i++ ){
+ PMTHit* aHit = (*HC)[i];
+ // fill vectors //
+ x.push_back ( aHit->getPos().x() );
+ z.push_back ( aHit->getPos().z() );
+ xHit.push_back ( aHit->getHit().x() );
+ zHit.push_back ( aHit->getHit().z() );
+ time.push_back ( aHit->getTime() );
+ G4double pTheta = atan(sqrt(pow(aHit->getMomentum().x(),2) + pow(aHit->getMomentum().z(),2) )/fabs(aHit->getMomentum().y() ));
+ theta.push_back( pTheta );
+ }
// fill ntuple //
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
- G4double pTheta = atan(sqrt(pow(aHit->getMomentum().x(),2) + pow(aHit->getMomentum().z(),2) )/fabs(aHit->getMomentum().y() ));
- analysisManager->FillNtupleDColumn(0, aHit->getPos().x() );
- analysisManager->FillNtupleDColumn(1, aHit->getPos().z() );
- analysisManager->FillNtupleDColumn(2, aHit->getHit().x() );
- analysisManager->FillNtupleDColumn(3, aHit->getHit().z() );
- analysisManager->FillNtupleDColumn(4, aHit->getTime() );
- analysisManager->FillNtupleDColumn(5, pTheta );
+ /* Not necessary if filling with vectors
+ analysisManager->FillNtupleDColumn(0, x );
+ analysisManager->FillNtupleDColumn(1, z );
+ analysisManager->FillNtupleDColumn(2, xHit );
+ analysisManager->FillNtupleDColumn(3, zHit );
+ analysisManager->FillNtupleDColumn(4, time );
+ analysisManager->FillNtupleDColumn(5, theta );
+ */
analysisManager->AddNtupleRow();
}
}
diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc
index af2cbe85de70203b6d3bace219729c289b2f4e9a..8933f22e77415d65776883d82afd2cf59ad472b4 100644
--- a/src/PrimaryGeneratorAction.cc
+++ b/src/PrimaryGeneratorAction.cc
@@ -25,23 +25,24 @@
/// \file /src/PrimaryGeneratorAction.cc
/// \brief Implementation of the PrimaryGeneratorAction class
#include "PrimaryGeneratorAction.hh"
+#include "lgAnalysis.hh"
#include "Randomize.hh"
#include "G4Event.hh"
-#include "G4GeneralParticleSource.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
-#include "G4SystemOfUnits.hh"
+
/*
*
*/
PrimaryGeneratorAction::PrimaryGeneratorAction()
- : G4VUserPrimaryGeneratorAction(),
- fParticleGun(0)
+ : G4VUserPrimaryGeneratorAction()
{
- fParticleGun = new G4GeneralParticleSource();
+ fParticleGun = new G4GeneralParticleSource();
+ fASCIIParticleGun = new ASCIIPrimaryGenerator();
+ fMessenger = new PrimaryGeneratorMessenger( this );
}
/*
@@ -50,6 +51,8 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
delete fParticleGun;
+ delete fASCIIParticleGun;
+ delete fMessenger;
}
/*
@@ -57,5 +60,14 @@ PrimaryGeneratorAction::~PrimaryGeneratorAction()
*/
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
- fParticleGun->GeneratePrimaryVertex(anEvent);
+ if(fUseInput){
+ fASCIIParticleGun->GeneratePrimaryVertex(anEvent);
+ }else{
+ fParticleGun->GeneratePrimaryVertex(anEvent);
+ }
+}
+
+void PrimaryGeneratorAction::SetInputFile(G4String _name){
+ fASCIIParticleGun->SetInputFile(_name);
+ fUseInput = true;
}
diff --git a/src/PrimaryGeneratorMessenger.cc b/src/PrimaryGeneratorMessenger.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e81b0c39068df91ac19ea61af1f58c18a3054f0c
--- /dev/null
+++ b/src/PrimaryGeneratorMessenger.cc
@@ -0,0 +1,73 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+/// \file /include/PrimaryGeneratorMessenger.cc
+/// \brief Implementation of the PrimaryGeneratorMessenger class
+//
+//
+
+#include "PrimaryGeneratorMessenger.hh"
+
+#include <sstream>
+#include <iostream>
+
+#include "G4OpticalSurface.hh"
+
+#include "PrimaryGeneratorAction.hh"
+#include "G4UIdirectory.hh"
+#include "G4UIcommand.hh"
+#include "G4UIparameter.hh"
+#include "G4UIcmdWithAString.hh"
+
+
+PrimaryGeneratorMessenger::PrimaryGeneratorMessenger(PrimaryGeneratorAction * Gen)
+:G4UImessenger(),fGenerator(Gen),fGeneratorDir(0),fInputFileCmd(0)
+{
+ fGeneratorDir = new G4UIdirectory("/Input/");
+ fGeneratorDir->SetGuidance("");
+
+ fInputFileCmd = new G4UIcmdWithAString("/Input/FileName", this);
+ fInputFileCmd->SetGuidance("Filepath of the previously generated events");
+ fInputFileCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
+ fInputFileCmd->SetToBeBroadcasted(false);
+
+}
+
+/*
+ *
+ */
+PrimaryGeneratorMessenger::~PrimaryGeneratorMessenger(){
+ delete fInputFileCmd;
+}
+
+/*
+ *
+ */
+void PrimaryGeneratorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
+{
+ if(command == fInputFileCmd){
+ fGenerator->SetInputFile(newValue);
+ }
+}
diff --git a/src/RunAction.cc b/src/RunAction.cc
index 3ac089f702ecc1e4ab3572782264b798e724a91b..e0f869057a6dfa3140c2738f403eed910bbb5f36 100644
--- a/src/RunAction.cc
+++ b/src/RunAction.cc
@@ -35,9 +35,11 @@
#include "G4Run.hh"
#include "RunAction.hh"
+#include "EventAction.hh"
+#include "lgAnalysis.hh"
+
#include "G4VAnalysisManager.hh"
#include "G4SystemOfUnits.hh"
-#include "lgAnalysis.hh"
/*
@@ -51,24 +53,35 @@ RunAction::RunAction(G4String fileName)
G4RunManager::GetRunManager()->SetPrintProgress(1);
// Create analysis manager. The choice of analysis
- // technology is done via selectin of a namespace
+ // technology is done via selection of a namespace
// in B4Analysis.hh
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
+ // Get event action
+ const EventAction* constEventAction
+ = static_cast<const EventAction*>(G4RunManager::GetRunManager()
+ ->GetUserEventAction());
+ EventAction* eventAction
+ = const_cast<EventAction*>(constEventAction);
+
+ std::vector< std::vector<double>* > pVec = eventAction->GetVectors( );
+
// Create directories
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
// Creating ntuple
analysisManager->CreateNtuple("lightGuide", "pos and momentum");
- analysisManager->CreateNtupleDColumn("X");
- analysisManager->CreateNtupleDColumn("Z");
- analysisManager->CreateNtupleDColumn("hitX");
- analysisManager->CreateNtupleDColumn("hitZ");
- analysisManager->CreateNtupleDColumn("time");
- analysisManager->CreateNtupleDColumn("theta");
+ analysisManager->CreateNtupleDColumn( "X", *pVec[0] );
+ analysisManager->CreateNtupleDColumn( "Z", *pVec[1] );
+ analysisManager->CreateNtupleDColumn( "hitX", *pVec[2] );
+ analysisManager->CreateNtupleDColumn( "hitZ", *pVec[3] );
+ analysisManager->CreateNtupleDColumn( "time", *pVec[4] );
+ analysisManager->CreateNtupleDColumn( "theta", *pVec[5] );
+ analysisManager->CreateNtupleIColumn( "EventNo" );
analysisManager->FinishNtuple();
+
}
/*