aboutsummaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/G4BeamTestDetectorConstruction.cxx193
-rw-r--r--src/G4BeamTestEMPhysics.cxx71
-rw-r--r--src/G4BeamTestGeneralPhysics.cxx62
-rw-r--r--src/G4BeamTestHadronPhysics.cxx425
-rw-r--r--src/G4BeamTestIonPhysics.cxx112
-rw-r--r--src/G4BeamTestMuonPhysics.cxx94
-rw-r--r--src/G4BeamTestPhysicsList.cxx81
-rw-r--r--src/G4BeamTestRunManager.cxx143
-rw-r--r--src/G4BeamTestTank.cxx382
-rw-r--r--src/G4BeamTestUserSteppingAction.cxx42
-rw-r--r--src/G4BeamTestUserTrackingAction.cxx50
-rw-r--r--src/G4Interface.cxx292
-rw-r--r--src/G4TankIceSD.cxx176
13 files changed, 2123 insertions, 0 deletions
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
new file mode 100644
index 0000000..6eeaa54
--- /dev/null
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -0,0 +1,193 @@
+#include <G4NistManager.hh>
+#include <G4Material.hh>
+
+#include <G4LogicalVolume.hh>
+#include <G4PVPlacement.hh>
+
+#include <G4Box.hh>
+
+#include <G4VisAttributes.hh>
+#include <G4UserLimits.hh>
+
+#include "G4BeamTestDetectorConstruction.h"
+#include "G4BeamTestTank.h"
+
+
+G4BeamTestDetectorConstruction::G4BeamTestDetectorConstruction():
+ origin_(NAN, NAN, NAN), verboseLevel_(0)/* , tankList_(0) */
+{
+}
+
+
+G4BeamTestDetectorConstruction::~G4BeamTestDetectorConstruction()
+{
+}
+
+
+G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
+{
+ /* if(tankList_.empty()) return NULL; */
+
+ CreateMaterials();
+
+ // Determine bottom z position of snow volume
+ G4double zSnowBottom(tankList_.at(0)->GetPos().z());
+ BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+ {
+ // z position of bottom of tank
+ G4double z = tank->GetPos().z() - 0.5*tank->GetTankHeight_G4();
+ zSnowBottom = std::min(z, zSnowBottom);
+ }
+
+ // Subtract safety margin
+ zSnowBottom -= 1.0*CLHEP::m;
+
+ // Triangulate snow surface
+ G4Delaunay delaunay;
+ BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+ {
+ // z position of snow surface
+ G4double z = tank->GetPos().z() + 0.5 * tank->GetTankHeight_G4() + tank->GetSnowHeight_G4();
+
+ delaunay.AddPoint(tank->GetDelaunayPoint1().x(),
+ tank->GetDelaunayPoint1().y(),
+ z - zSnowBottom);
+ delaunay.AddPoint(tank->GetDelaunayPoint2().x(),
+ tank->GetDelaunayPoint2().y(),
+ z - zSnowBottom);
+ delaunay.AddPoint(tank->GetDelaunayPoint3().x(),
+ tank->GetDelaunayPoint3().y(),
+ z - zSnowBottom);
+ }
+
+ // Create tesselated snow volume
+ G4TessellatedSolid* solidSnow = new G4TessellatedSolid("solid_snow");
+ delaunay.BuildSolid(solidSnow, 50.0*CLHEP::m, 100.0*CLHEP::m);
+
+ // Determine World dimensions
+ G4double xHalfLength = 0.5 * (delaunay.GetXmax() - delaunay.GetXmin());
+ G4double yHalfLength = 0.5 * (delaunay.GetYmax() - delaunay.GetYmin());
+ G4double zHalfLength = 0.5 * (delaunay.GetZmax() + 20.0 * CLHEP::m); // 20 m of atmosphere
+
+ // World origin in IceCube coordinates
+ origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength);
+
+ // Create world volume
+ G4Box* world_box = new G4Box("solid_world", xHalfLength, yHalfLength, zHalfLength);
+ G4LogicalVolume* worldLog =
+ new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0);
+ G4VPhysicalVolume* worldPhys =
+ new G4PVPlacement(0, G4ThreeVector(), worldLog, "world", 0, false, 0);
+
+ // Snow layer
+ G4LogicalVolume* snowLog =
+ new G4LogicalVolume(solidSnow, G4Material::GetMaterial("Snow"), "log_snow", 0, 0, 0);
+ G4VPhysicalVolume* snowPhys =
+ new G4PVPlacement(0, G4ThreeVector(0, 0, -zHalfLength), snowLog, "snow", worldLog, false, 0);
+
+ // Instantiation of a set of visualization attributes with cyan colour
+ G4VisAttributes * snowVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
+ // Assignment of the visualization attributes to the logical volume
+ snowLog->SetVisAttributes(snowVisAtt);
+
+ // Install tanks
+ BOOST_FOREACH(G4BeamTestTank* tank, tankList_)
+ {
+ tank->InstallTank(snowPhys, origin_);
+ }
+
+ // User limits (energy cutoffs)
+ // Do not create photons or electrons below cherenkov threshold
+ // See also corresponding UserSpecialCuts in Physicslist !!!!
+ G4UserLimits* energyLimit = new G4UserLimits();
+ energyLimit->SetUserMinEkine(280.0 * CLHEP::keV); // Cherenkov threshold of electrons in ice TODO(shivesh)
+ worldLog->SetUserLimits(energyLimit);
+ snowLog->SetUserLimits(energyLimit);
+
+ return worldPhys;
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateMaterials()
+{
+ CreateAir();
+ /* CreateIce(); */
+ /* CreateSnow(); */
+ CreateWater();
+ CreatePlastic();
+ /* CreatePerlite(); */
+ CreateGlassSphere();
+ CreateEffectiveDOMMaterial();
+
+ //if(verboseLevel_>0) G4cout << *G4Material::GetMaterialTable() << G4endl;
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateAir()
+{
+ G4NistManager* nistManager = G4NistManager::Instance();
+ nistManager->ConstructNewGasMaterial("Air","G4_AIR"/* , (273.15 - 40.0)*CLHEP::kelvin, 670.0E-3*CLHEP::bar */);
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateWater()
+{
+ G4NistManager* nistManager = G4NistManager::Instance();
+ // G4Material* ice = new G4Material("Water", 1.0 * CLHEP::g / CLHEP::cm3, 2, kStateLiquid);
+ // ice->AddElement(nistManager->FindOrBuildElement("H"), 2);
+ // ice->AddElement(nistManager->FindOrBuildElement("O"), 1);
+ nistManager->ConstructNewGasMaterial("Water","G4_WATER");
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreatePlastic()
+{
+ G4NistManager* nistManager = G4NistManager::Instance();
+ /*
+ G4Material* plastic = new G4Material("Plastic", 1.5 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
+ plastic->AddElement(nistManager->FindOrBuildElement("H"), 20);
+ plastic->AddElement(nistManager->FindOrBuildElement("C"), 10);
+ plastic->AddElement(nistManager->FindOrBuildElement("O"), 5);
+ */
+
+ // POM
+ G4Material* plastic = new G4Material("Plastic", 1.425 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
+ plastic->AddElement(nistManager->FindOrBuildElement("H"), 2);
+ plastic->AddElement(nistManager->FindOrBuildElement("C"), 1);
+ plastic->AddElement(nistManager->FindOrBuildElement("O"), 1);
+
+ //nistManager->FindOrBuildMaterial("G4_POLYOXYMETHYLENE");
+
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateGlassSphere()
+{
+ // Elemental composition not exact for detailed line-up look at
+ // docushare collection 690 -> XRF and ICP Analysis of Pressure Sphere Glass
+
+ // 20 lbs. weight = 9072 g
+ // 6.5" outer radius & 0.5" thickness = 4024 cm3
+ G4NistManager* nistManager = G4NistManager::Instance();
+ G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
+ glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
+ glass->AddElement(nistManager->FindOrBuildElement("O"), 2);
+}
+
+/*****************************************************************/
+
+void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial()
+{
+ // Mass of the a complete DOM: 12050 g
+ // Mass without glass sphere: 2978 g
+ // 6" inner glass radius: volume = 14830 cm3
+ G4NistManager* nistManager = G4NistManager::Instance();
+ G4Material* glass = new G4Material("effectiveDOM", 0.2 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
+ glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
+ glass->AddElement(nistManager->FindOrBuildElement("O"), 2);
+}
diff --git a/src/G4BeamTestEMPhysics.cxx b/src/G4BeamTestEMPhysics.cxx
new file mode 100644
index 0000000..d36207a
--- /dev/null
+++ b/src/G4BeamTestEMPhysics.cxx
@@ -0,0 +1,71 @@
+/*
+ * copyright (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestEMPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+
+#include <globals.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4Gamma.hh>
+#include <G4Electron.hh>
+#include <G4Positron.hh>
+#include <G4NeutrinoE.hh>
+#include <G4AntiNeutrinoE.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestEMPhysics.h"
+
+
+G4BeamTestEMPhysics::G4BeamTestEMPhysics()
+ : G4VPhysicsConstructor("standard")
+{}
+
+
+G4BeamTestEMPhysics::~G4BeamTestEMPhysics()
+{}
+
+
+void G4BeamTestEMPhysics::ConstructParticle()
+{
+ // gamma
+ G4Gamma::GammaDefinition();
+
+ // electron
+ G4Electron::ElectronDefinition();
+ G4Positron::PositronDefinition();
+ G4NeutrinoE::NeutrinoEDefinition();
+ G4AntiNeutrinoE::AntiNeutrinoEDefinition();
+}
+
+
+void G4BeamTestEMPhysics::ConstructProcess()
+{
+ G4ProcessManager *pManager = 0;
+
+ // Gamma Physics
+ pManager = G4Gamma::Gamma()->GetProcessManager();
+ pManager->AddDiscreteProcess(&photoEffect);
+ pManager->AddDiscreteProcess(&comptonEffect);
+ pManager->AddDiscreteProcess(&pairProduction);
+
+ // Electron Physics
+ pManager = G4Electron::Electron()->GetProcessManager();
+ pManager->AddProcess(&electronMultipleScattering, -1, 1, 1);
+ pManager->AddProcess(&electronIonisation, -1, 2, 2);
+ pManager->AddProcess(&electronBremsStrahlung, -1, 3, 3);
+
+ //Positron Physics
+ pManager = G4Positron::Positron()->GetProcessManager();
+ pManager->AddProcess(&positronMultipleScattering, -1, 1, 1);
+ pManager->AddProcess(&positronIonisation, -1, 2, 2);
+ pManager->AddProcess(&positronBremsStrahlung, -1, 3, 3);
+ pManager->AddProcess(&annihilation, 0,-1, 4);
+}
diff --git a/src/G4BeamTestGeneralPhysics.cxx b/src/G4BeamTestGeneralPhysics.cxx
new file mode 100644
index 0000000..2def9b4
--- /dev/null
+++ b/src/G4BeamTestGeneralPhysics.cxx
@@ -0,0 +1,62 @@
+/*
+ * copyright (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestGeneralPhysics.cxx 152849 2017-01-20 21:44:25Z jgonzalez $
+ *
+ * @version $Revision: 152849 $
+ * @date $LastChangedDate: 2017-01-20 21:44:25 +0000 (Fri, 20 Jan 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ */
+
+#include "G4BeamTestGeneralPhysics.h"
+
+#include <iomanip>
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ProcessManager.hh>
+// Bosons
+#include <G4ChargedGeantino.hh>
+#include <G4Geantino.hh>
+#include <G4Version.hh>
+
+G4BeamTestGeneralPhysics::G4BeamTestGeneralPhysics()
+ : G4VPhysicsConstructor("general")
+{}
+
+
+G4BeamTestGeneralPhysics::~G4BeamTestGeneralPhysics()
+{}
+
+
+void G4BeamTestGeneralPhysics::ConstructParticle()
+{
+ // pseudo-particles
+ G4Geantino::GeantinoDefinition();
+ G4ChargedGeantino::ChargedGeantinoDefinition();
+}
+
+
+void G4BeamTestGeneralPhysics::ConstructProcess()
+{
+ //AddTransportation();
+
+// Decay processes are set somewhere else now
+#if G4VERSION_NUMBER < 1000
+ // Add Decay Process
+ theParticleIterator->reset();
+ while ((*theParticleIterator)()) {
+ G4ParticleDefinition *particle = theParticleIterator->value();
+ G4ProcessManager *pmanager = particle->GetProcessManager();
+ if (decay.IsApplicable(*particle)) {
+ pmanager->AddProcess(&decay);
+ // set ordering for PostStepDoIt and AtRestDoIt
+ pmanager->SetProcessOrdering(&decay, idxPostStep);
+ pmanager->SetProcessOrdering(&decay, idxAtRest);
+ }
+ }
+#endif
+}
+
diff --git a/src/G4BeamTestHadronPhysics.cxx b/src/G4BeamTestHadronPhysics.cxx
new file mode 100644
index 0000000..5150b67
--- /dev/null
+++ b/src/G4BeamTestHadronPhysics.cxx
@@ -0,0 +1,425 @@
+/*
+ * copyright (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestHadronPhysics.cxx 152814 2017-01-19 21:34:52Z jgonzalez $
+ *
+ * @version $Revision: 152814 $
+ * @date $LastChangedDate: 2017-01-19 21:34:52 +0000 (Thu, 19 Jan 2017) $
+ * @author Fabian Kislat <fabian.kislat@desy.de> Last changed by: $LastChangedBy: jgonzalez $
+ *
+ * Copied from tanktop
+ */
+
+#include <iomanip>
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4MesonConstructor.hh>
+#include <G4BaryonConstructor.hh>
+#include <G4ShortLivedConstructor.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestHadronPhysics.h"
+
+G4BeamTestHadronPhysics::G4BeamTestHadronPhysics()
+ : G4VPhysicsConstructor("hadron") {
+}
+
+/********************************************************************/
+
+G4BeamTestHadronPhysics::~G4BeamTestHadronPhysics()
+{
+ delete stringDecay_;
+}
+
+/********************************************************************/
+
+void G4BeamTestHadronPhysics::ConstructParticle()
+{
+
+ // Construct all mesons
+ G4MesonConstructor mesonConstructor;
+ mesonConstructor.ConstructParticle();
+
+ // Construct all barions
+ G4BaryonConstructor baryonConstructor;
+ baryonConstructor.ConstructParticle();
+
+ // Construct resonaces and quarks
+ G4ShortLivedConstructor shortLivedConstructor;
+ shortLivedConstructor.ConstructParticle();
+
+}
+
+/********************************************************************/
+
+void G4BeamTestHadronPhysics::ConstructProcess()
+{
+ G4ProcessManager *pManager = 0;
+
+ // G4cout << "" << G4endl;
+ // G4cout << "You are using the G4BeamTestHadronPhysics" << G4endl;
+ // G4cout << " - Note that this hadronic physics list is not optimized for any particular usage" << G4endl;
+ // G4cout << " - If you wish to have a starting point tailored for a particular area of work," << G4endl;
+ // G4cout << " please use one of the available physics lists by use-case." << G4endl;
+ // G4cout << " More information can also be found from the Geant4 HyperNews." << G4endl;
+ // G4cout << "" << G4endl;
+
+ // Elastic Process
+ elasticModel_ = new G4LElastic();
+ elasticProcess_.RegisterMe(elasticModel_);
+
+ // pi+ and pi-
+ preEquilib_ = new G4PreCompoundModel(&excitationHandler_);
+ cascade_.SetDeExcitation(preEquilib_);
+ theoModel_.SetTransport(&cascade_);
+ theoModel_.SetHighEnergyGenerator(&stringModel_);
+ stringDecay_ = new G4ExcitedStringDecay(&fragmentation_);
+ stringModel_.SetFragmentationModel(stringDecay_);
+ theoModel_.SetMinEnergy(15*CLHEP::GeV);
+ theoModel_.SetMaxEnergy(100*CLHEP::TeV);
+
+ // PionPlus
+ pManager = G4PionPlus::PionPlus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ lePionPlusModel_ = new G4LEPionPlusInelastic();
+ pionPlusInelastic_.RegisterMe(lePionPlusModel_);
+ pionPlusInelastic_.RegisterMe(&theoModel_);
+ pManager->AddDiscreteProcess(&pionPlusInelastic_);
+
+ pManager->AddProcess(&pionPlusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&pionPlusMult_);
+ pManager->SetProcessOrdering(&pionPlusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&pionPlusMult_, idxPostStep, 1);
+
+ // PionMinus
+ pManager = G4PionMinus::PionMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ lePionMinusModel_ = new G4LEPionMinusInelastic();
+ pionMinusInelastic_.RegisterMe(lePionMinusModel_);
+ pionMinusInelastic_.RegisterMe(&theoModel_);
+ pManager->AddDiscreteProcess(&pionMinusInelastic_);
+
+ pManager->AddProcess(&pionMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&pionMinusMult_);
+ pManager->SetProcessOrdering(&pionMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&pionMinusMult_, idxPostStep, 1);
+
+ pManager->AddRestProcess(&pionMinusAbsorption_, ordDefault);
+
+ // KaonPlus
+ pManager = G4KaonPlus::KaonPlus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leKaonPlusModel_ = new G4LEKaonPlusInelastic();
+ heKaonPlusModel_ = new G4HEKaonPlusInelastic();
+ kaonPlusInelastic_.RegisterMe(leKaonPlusModel_);
+ kaonPlusInelastic_.RegisterMe(&theoModel_);
+ pManager->AddDiscreteProcess(&kaonPlusInelastic_);
+
+ pManager->AddProcess(&kaonPlusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&kaonPlusMult_);
+ pManager->SetProcessOrdering(&kaonPlusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&kaonPlusMult_, idxPostStep, 1);
+
+ // KaonMinus
+ pManager = G4KaonMinus::KaonMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leKaonMinusModel_ = new G4LEKaonMinusInelastic();
+ heKaonMinusModel_ = new G4HEKaonMinusInelastic();
+ kaonMinusInelastic_.RegisterMe(leKaonMinusModel_);
+ kaonMinusInelastic_.RegisterMe(heKaonMinusModel_);
+ pManager->AddDiscreteProcess(&kaonMinusInelastic_);
+
+ pManager->AddProcess(&kaonMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&kaonMinusMult_);
+ pManager->SetProcessOrdering(&kaonMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&kaonMinusMult_, idxPostStep, 1);
+
+ pManager->AddRestProcess(&kaonMinusAbsorption_, ordDefault);
+
+ // KaonZeroL
+ pManager = G4KaonZeroLong::KaonZeroLong()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leKaonZeroLModel_ = new G4LEKaonZeroLInelastic();
+ heKaonZeroLModel_ = new G4HEKaonZeroInelastic();
+ kaonZeroLInelastic_.RegisterMe(leKaonZeroLModel_);
+ kaonZeroLInelastic_.RegisterMe(heKaonZeroLModel_);
+ pManager->AddDiscreteProcess(&kaonZeroLInelastic_);
+
+ // KaonZeroS
+ pManager = G4KaonZeroShort::KaonZeroShort()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leKaonZeroSModel_ = new G4LEKaonZeroSInelastic();
+ heKaonZeroSModel_ = new G4HEKaonZeroInelastic();
+ kaonZeroSInelastic_.RegisterMe(leKaonZeroSModel_);
+ kaonZeroSInelastic_.RegisterMe(heKaonZeroSModel_);
+ pManager->AddDiscreteProcess(&kaonZeroSInelastic_);
+
+ // Proton
+ pManager = G4Proton::Proton()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leProtonModel_ = new G4LEProtonInelastic();
+ heProtonModel_ = new G4HEProtonInelastic();
+ protonInelastic_.RegisterMe(leProtonModel_);
+ protonInelastic_.RegisterMe(&theoModel_);
+ pManager->AddDiscreteProcess(&protonInelastic_);
+
+ pManager->AddProcess(&protonIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&protonMult_);
+ pManager->SetProcessOrdering(&protonMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&protonMult_, idxPostStep, 1);
+
+ // anti-Proton
+ pManager = G4AntiProton::AntiProton()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiProtonModel_ = new G4LEAntiProtonInelastic();
+ heAntiProtonModel_ = new G4HEAntiProtonInelastic();
+ antiProtonInelastic_.RegisterMe(leAntiProtonModel_);
+ antiProtonInelastic_.RegisterMe(heAntiProtonModel_);
+ pManager->AddDiscreteProcess(&antiProtonInelastic_);
+
+ pManager->AddProcess(&antiProtonIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&antiProtonMult_);
+ pManager->SetProcessOrdering(&antiProtonMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&antiProtonMult_, idxPostStep, 1);
+
+ pManager->AddRestProcess(&antiProtonAnnihilation_);
+
+ // Neutron
+ pManager = G4Neutron::Neutron()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leNeutronModel_ = new G4LENeutronInelastic();
+ heNeutronModel_ = new G4HENeutronInelastic();
+ neutronInelastic_.RegisterMe(leNeutronModel_);
+ neutronInelastic_.RegisterMe(&theoModel_);
+ pManager->AddDiscreteProcess(&neutronInelastic_);
+
+ neutronFissionModel_ = new G4LFission();
+ neutronFission_.RegisterMe(neutronFissionModel_);
+ pManager->AddDiscreteProcess(&neutronFission_);
+
+ neutronCaptureModel_ = new G4LCapture();
+ neutronCapture_.RegisterMe(neutronCaptureModel_);
+ pManager->AddDiscreteProcess(&neutronCapture_);
+
+ // AntiNeutron
+ pManager = G4AntiNeutron::AntiNeutron()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiNeutronModel_ = new G4LEAntiNeutronInelastic();
+ heAntiNeutronModel_ = new G4HEAntiNeutronInelastic();
+ antiNeutronInelastic_.RegisterMe(leAntiNeutronModel_);
+ antiNeutronInelastic_.RegisterMe(heAntiNeutronModel_);
+ pManager->AddDiscreteProcess(&antiNeutronInelastic_);
+
+ pManager->AddRestProcess(&antiNeutronAnnihilation_);
+
+ // Lambda
+ pManager = G4Lambda::Lambda()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leLambdaModel_ = new G4LELambdaInelastic();
+ heLambdaModel_ = new G4HELambdaInelastic();
+ lambdaInelastic_.RegisterMe(leLambdaModel_);
+ lambdaInelastic_.RegisterMe(heLambdaModel_);
+ pManager->AddDiscreteProcess(&lambdaInelastic_);
+
+ // AntiLambda
+ pManager = G4AntiLambda::AntiLambda()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiLambdaModel_ = new G4LEAntiLambdaInelastic();
+ heAntiLambdaModel_ = new G4HEAntiLambdaInelastic();
+ antiLambdaInelastic_.RegisterMe(leAntiLambdaModel_);
+ antiLambdaInelastic_.RegisterMe(heAntiLambdaModel_);
+ pManager->AddDiscreteProcess(&antiLambdaInelastic_);
+
+ // SigmaMinus
+ pManager = G4SigmaMinus::SigmaMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leSigmaMinusModel_ = new G4LESigmaMinusInelastic();
+ heSigmaMinusModel_ = new G4HESigmaMinusInelastic();
+ sigmaMinusInelastic_.RegisterMe(leSigmaMinusModel_);
+ sigmaMinusInelastic_.RegisterMe(heSigmaMinusModel_);
+ pManager->AddDiscreteProcess(&sigmaMinusInelastic_);
+
+ pManager->AddProcess(&sigmaMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&sigmaMinusMult_);
+ pManager->SetProcessOrdering(&sigmaMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&sigmaMinusMult_, idxPostStep, 1);
+
+ // anti-SigmaMinus
+ pManager = G4AntiSigmaMinus::AntiSigmaMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiSigmaMinusModel_ = new G4LEAntiSigmaMinusInelastic();
+ heAntiSigmaMinusModel_ = new G4HEAntiSigmaMinusInelastic();
+ antiSigmaMinusInelastic_.RegisterMe(leAntiSigmaMinusModel_);
+ antiSigmaMinusInelastic_.RegisterMe(heAntiSigmaMinusModel_);
+ pManager->AddDiscreteProcess(&antiSigmaMinusInelastic_);
+
+ pManager->AddProcess(&antiSigmaMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&antiSigmaMinusMult_);
+ pManager->SetProcessOrdering(&antiSigmaMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&antiSigmaMinusMult_, idxPostStep, 1);
+
+ // SigmaPlus
+ pManager = G4SigmaPlus::SigmaPlus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leSigmaPlusModel_ = new G4LESigmaPlusInelastic();
+ heSigmaPlusModel_ = new G4HESigmaPlusInelastic();
+ sigmaPlusInelastic_.RegisterMe(leSigmaPlusModel_);
+ sigmaPlusInelastic_.RegisterMe(heSigmaPlusModel_);
+ pManager->AddDiscreteProcess(&sigmaPlusInelastic_);
+
+ pManager->AddProcess(&sigmaPlusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&sigmaPlusMult_);
+ pManager->SetProcessOrdering(&sigmaPlusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&sigmaPlusMult_, idxPostStep, 1);
+
+ // anti-SigmaPlus
+ pManager = G4AntiSigmaPlus::AntiSigmaPlus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiSigmaPlusModel_ = new G4LEAntiSigmaPlusInelastic();
+ heAntiSigmaPlusModel_ = new G4HEAntiSigmaPlusInelastic();
+ antiSigmaPlusInelastic_.RegisterMe(leAntiSigmaPlusModel_);
+ antiSigmaPlusInelastic_.RegisterMe(heAntiSigmaPlusModel_);
+ pManager->AddDiscreteProcess(&antiSigmaPlusInelastic_);
+
+ pManager->AddProcess(&antiSigmaPlusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&antiSigmaPlusMult_);
+ pManager->SetProcessOrdering(&antiSigmaPlusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&antiSigmaPlusMult_, idxPostStep, 1);
+
+ // XiMinus
+ pManager = G4XiMinus::XiMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leXiMinusModel_ = new G4LEXiMinusInelastic();
+ heXiMinusModel_ = new G4HEXiMinusInelastic();
+ xiMinusInelastic_.RegisterMe(leXiMinusModel_);
+ xiMinusInelastic_.RegisterMe(heXiMinusModel_);
+ pManager->AddDiscreteProcess(&xiMinusInelastic_);
+
+ pManager->AddProcess(&xiMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&xiMinusMult_);
+ pManager->SetProcessOrdering(&xiMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&xiMinusMult_, idxPostStep, 1);
+
+ // anti-XiMinus
+ pManager = G4AntiXiMinus::AntiXiMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiXiMinusModel_ = new G4LEAntiXiMinusInelastic();
+ heAntiXiMinusModel_ = new G4HEAntiXiMinusInelastic();
+ antiXiMinusInelastic_.RegisterMe(leAntiXiMinusModel_);
+ antiXiMinusInelastic_.RegisterMe(heAntiXiMinusModel_);
+ pManager->AddDiscreteProcess(&antiXiMinusInelastic_);
+
+ pManager->AddProcess(&antiXiMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&antiXiMinusMult_);
+ pManager->SetProcessOrdering(&antiXiMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&antiXiMinusMult_, idxPostStep, 1);
+
+ // XiZero
+ pManager = G4XiZero::XiZero()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leXiZeroModel_ = new G4LEXiZeroInelastic();
+ heXiZeroModel_ = new G4HEXiZeroInelastic();
+ xiZeroInelastic_.RegisterMe(leXiZeroModel_);
+ xiZeroInelastic_.RegisterMe(heXiZeroModel_);
+ pManager->AddDiscreteProcess(&xiZeroInelastic_);
+
+ // anti-XiZero
+ pManager = G4AntiXiZero::AntiXiZero()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiXiZeroModel_ = new G4LEAntiXiZeroInelastic();
+ heAntiXiZeroModel_ = new G4HEAntiXiZeroInelastic();
+ antiXiZeroInelastic_.RegisterMe(leAntiXiZeroModel_);
+ antiXiZeroInelastic_.RegisterMe(heAntiXiZeroModel_);
+ pManager->AddDiscreteProcess(&antiXiZeroInelastic_);
+
+ // OmegaMinus
+ pManager = G4OmegaMinus::OmegaMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leOmegaMinusModel_ = new G4LEOmegaMinusInelastic();
+ heOmegaMinusModel_ = new G4HEOmegaMinusInelastic();
+ omegaMinusInelastic_.RegisterMe(leOmegaMinusModel_);
+ omegaMinusInelastic_.RegisterMe(heOmegaMinusModel_);
+ pManager->AddDiscreteProcess(&omegaMinusInelastic_);
+
+ pManager->AddProcess(&omegaMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&omegaMinusMult_);
+ pManager->SetProcessOrdering(&omegaMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&omegaMinusMult_, idxPostStep, 1);
+
+ // anti-OmegaMinus
+ pManager = G4AntiOmegaMinus::AntiOmegaMinus()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ leAntiOmegaMinusModel_ = new G4LEAntiOmegaMinusInelastic();
+ heAntiOmegaMinusModel_ = new G4HEAntiOmegaMinusInelastic();
+ antiOmegaMinusInelastic_.RegisterMe(leAntiOmegaMinusModel_);
+ antiOmegaMinusInelastic_.RegisterMe(heAntiOmegaMinusModel_);
+ pManager->AddDiscreteProcess(&antiOmegaMinusInelastic_);
+
+ pManager->AddProcess(&antiOmegaMinusIonisation_, ordInActive,2, 2);
+
+ pManager->AddProcess(&antiOmegaMinusMult_);
+ pManager->SetProcessOrdering(&antiOmegaMinusMult_, idxAlongStep, 1);
+ pManager->SetProcessOrdering(&antiOmegaMinusMult_, idxPostStep, 1);
+
+}
diff --git a/src/G4BeamTestIonPhysics.cxx b/src/G4BeamTestIonPhysics.cxx
new file mode 100644
index 0000000..3c7e553
--- /dev/null
+++ b/src/G4BeamTestIonPhysics.cxx
@@ -0,0 +1,112 @@
+/*
+ * copyright (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestIonPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+#include <iomanip>
+
+#include <globals.hh>
+#include <G4ios.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4IonConstructor.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestIonPhysics.h"
+
+G4BeamTestIonPhysics::G4BeamTestIonPhysics()
+ : G4VPhysicsConstructor ("ion")
+{}
+
+/********************************************************************/
+
+G4BeamTestIonPhysics::~G4BeamTestIonPhysics()
+{}
+
+
+/********************************************************************/
+
+void G4BeamTestIonPhysics::ConstructParticle()
+{
+ // Construct light ions
+ G4IonConstructor pConstructor;
+ pConstructor.ConstructParticle();
+}
+
+/********************************************************************/
+
+
+
+
+void G4BeamTestIonPhysics::ConstructProcess()
+{
+ G4ProcessManager *pManager = 0;
+
+ // Elastic Process
+ elasticModel_ = new G4LElastic;
+ elasticProcess_.RegisterMe(elasticModel_);
+
+ // Generic Ion
+ pManager = G4GenericIon::GenericIon()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ pManager->AddProcess(&ionMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&ionIonisation_, -1, 2, 2);
+
+ // Deuteron
+ pManager = G4Deuteron::Deuteron()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ deuteronModel_ = new G4LEDeuteronInelastic;
+ deuteronProcess_.RegisterMe(deuteronModel_);
+ pManager->AddDiscreteProcess(&deuteronProcess_);
+
+ pManager->AddProcess(&deuteronMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&deuteronIonisation_, -1, 2, 2);
+
+ // Triton
+ pManager = G4Triton::Triton()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ tritonModel_ = new G4LETritonInelastic;
+ tritonProcess_.RegisterMe(tritonModel_);
+ pManager->AddDiscreteProcess(&tritonProcess_);
+
+ pManager->AddProcess(&tritonMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&tritonIonisation_, -1, 2, 2);
+
+ // Alpha
+ pManager = G4Alpha::Alpha()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ alphaModel_ = new G4LEAlphaInelastic;
+ alphaProcess_.RegisterMe(alphaModel_);
+ pManager->AddDiscreteProcess(&alphaProcess_);
+
+ pManager->AddProcess(&alphaMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&alphaIonisation_, -1, 2, 2);
+
+ // He3
+ pManager = G4He3::He3()->GetProcessManager();
+ // add process
+ pManager->AddDiscreteProcess(&elasticProcess_);
+
+ pManager->AddProcess(&he3MultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&he3Ionisation_, -1, 2, 2);
+
+
+}
+
+
+
diff --git a/src/G4BeamTestMuonPhysics.cxx b/src/G4BeamTestMuonPhysics.cxx
new file mode 100644
index 0000000..92264a2
--- /dev/null
+++ b/src/G4BeamTestMuonPhysics.cxx
@@ -0,0 +1,94 @@
+/*
+ * copyright (C) 2010
+ * The Icecube Collaboration
+ *
+ * $Id: G4BeamTestMuonPhysics.cxx 86420 2012-03-20 16:00:37Z kislat $
+ *
+ * @version $Revision: 86420 $
+ * @date $LastChangedDate: 2012-03-20 16:00:37 +0000 (Tue, 20 Mar 2012) $
+ * @author Fabian Kislat <fabian.kislat@desy.de>, Peter Nießen (tanktop)
+ * Last changed by: $LastChangedBy: kislat $
+ */
+
+
+#include <globals.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4MuonPlus.hh>
+#include <G4MuonMinus.hh>
+#include <G4TauMinus.hh>
+#include <G4TauPlus.hh>
+#include <G4NeutrinoTau.hh>
+#include <G4AntiNeutrinoTau.hh>
+#include <G4NeutrinoMu.hh>
+#include <G4AntiNeutrinoMu.hh>
+#include <G4ProcessManager.hh>
+
+#include "G4BeamTestMuonPhysics.h"
+
+
+
+G4BeamTestMuonPhysics::G4BeamTestMuonPhysics()
+ : G4VPhysicsConstructor("muon")
+{}
+
+/********************************************************************/
+
+G4BeamTestMuonPhysics::~G4BeamTestMuonPhysics()
+{}
+
+/********************************************************************/
+
+void G4BeamTestMuonPhysics::ConstructParticle()
+{
+ // Mu
+ G4MuonPlus::MuonPlusDefinition();
+ G4MuonMinus::MuonMinusDefinition();
+ G4NeutrinoMu::NeutrinoMuDefinition();
+ G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
+
+ // Tau
+ G4TauMinus::TauMinusDefinition();
+ G4TauPlus::TauPlusDefinition();
+ G4NeutrinoTau::NeutrinoTauDefinition();
+ G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();
+}
+
+/********************************************************************/
+
+
+void G4BeamTestMuonPhysics::ConstructProcess()
+{
+
+ G4ProcessManager *pManager = 0;
+
+ // Muon Plus Physics
+ pManager = G4MuonPlus::MuonPlus()->GetProcessManager();
+
+ pManager->AddProcess(&muPlusMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&muPlusIonisation_, -1, 2, 2);
+ pManager->AddProcess(&muPlusBremsstrahlung_, -1, 3, 3);
+ pManager->AddProcess(&muPlusPairProduction_, -1, 4, 4);
+
+ // Muon Minus Physics
+ pManager = G4MuonMinus::MuonMinus()->GetProcessManager();
+
+ pManager->AddProcess(&muMinusMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&muMinusIonisation_, -1, 2, 2);
+ pManager->AddProcess(&muMinusBremsstrahlung_, -1, 3, 3);
+ pManager->AddProcess(&muMinusPairProduction_, -1, 4, 4);
+
+ pManager->AddRestProcess(&muMinusCaptureAtRest_);
+
+ // Tau Plus Physics
+ pManager = G4TauPlus::TauPlus()->GetProcessManager();
+
+ pManager->AddProcess(&tauPlusMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&tauPlusIonisation_, -1, 2, 2);
+
+ // Tau Minus Physics
+ pManager = G4TauMinus::TauMinus()->GetProcessManager();
+
+ pManager->AddProcess(&tauMinusMultipleScattering_, -1, 1, 1);
+ pManager->AddProcess(&tauMinusIonisation_, -1, 2, 2);
+}
diff --git a/src/G4BeamTestPhysicsList.cxx b/src/G4BeamTestPhysicsList.cxx
new file mode 100644
index 0000000..5a5befd
--- /dev/null
+++ b/src/G4BeamTestPhysicsList.cxx
@@ -0,0 +1,81 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestPhysicsList.cxx 152849 2017-01-20 21:44:25Z jgonzalez $
+ *
+ * @file G4BeamTestPhysicsList.cxx
+ * @version $Rev: 152849 $
+ * @date $Date: 2017-01-20 21:44:25 +0000 (Fri, 20 Jan 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig, Javier Gonzalez
+ */
+
+
+#include <globals.hh>
+#include <G4Version.hh>
+#include "G4BeamTestPhysicsList.h"
+#include "G4BeamTestGeneralPhysics.h"
+#if G4VERSION_NUMBER < 1000
+#include "G4BeamTestEMPhysics.h"
+#include "G4BeamTestMuonPhysics.h"
+#include "G4BeamTestHadronPhysics.h"
+#include "G4BeamTestIonPhysics.h"
+#else
+#include "G4DecayPhysics.hh"
+#include "G4EmStandardPhysics.hh"
+#include "G4EmExtraPhysics.hh"
+#include "G4IonPhysics.hh"
+#include "G4StoppingPhysics.hh"
+#include "G4HadronElasticPhysics.hh"
+//#include "G4NeutronTrackingCut.hh"
+#include "G4DataQuestionaire.hh"
+#include "G4HadronPhysicsFTFP_BERT.hh"
+#include <FTFP_BERT.hh>
+#endif
+
+#include <G4ProcessManager.hh>
+#include <G4ParticleTypes.hh>
+#include <G4UserSpecialCuts.hh>
+
+
+G4BeamTestPhysicsList::G4BeamTestPhysicsList()
+ : G4VUserPhysicsList()
+{
+ defaultCutValue = 0.7*CLHEP::mm;
+ SetVerboseLevel(1);
+
+ RegisterPhysics(new G4BeamTestGeneralPhysics);
+#if G4VERSION_NUMBER < 1000
+ RegisterPhysics(new G4BeamTestEMPhysics);
+ RegisterPhysics(new G4BeamTestMuonPhysics);
+ RegisterPhysics(new G4BeamTestHadronPhysics);
+ RegisterPhysics(new G4BeamTestIonPhysics);
+#else
+ // The following is basically Geant4's FTFP_BERT physics list
+ G4DataQuestionaire it(photon); // this checks that G4LEVELGAMMADATA is there
+
+ RegisterPhysics(new G4EmStandardPhysics());
+ RegisterPhysics(new G4EmExtraPhysics());
+ RegisterPhysics(new G4DecayPhysics());
+ RegisterPhysics(new G4HadronElasticPhysics());
+ RegisterPhysics(new G4HadronPhysicsFTFP_BERT());
+ RegisterPhysics(new G4StoppingPhysics());
+ RegisterPhysics(new G4IonPhysics());
+ //RegisterPhysics(new G4NeutronTrackingCut());
+#endif
+}
+
+
+G4BeamTestPhysicsList::~G4BeamTestPhysicsList()
+{
+}
+
+
+void G4BeamTestPhysicsList::SetCuts()
+{
+ //G4VUserPhysicsList::SetCutsWithDefault method sets
+ //the default cut value for all particle types
+ //
+ SetCutsWithDefault();
+ if (verboseLevel>0) DumpCutValuesTable();
+}
+
diff --git a/src/G4BeamTestRunManager.cxx b/src/G4BeamTestRunManager.cxx
new file mode 100644
index 0000000..753deb8
--- /dev/null
+++ b/src/G4BeamTestRunManager.cxx
@@ -0,0 +1,143 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestRunManager.cxx 162348 2018-04-25 20:09:46Z nega $
+ *
+ * @file G4BeamTestRunManager.cxx
+ * @version $Rev: 162348 $
+ * @date $Date: 2018-04-25 21:09:46 +0100 (Wed, 25 Apr 2018) $
+ * @author Tilo Waldenmaier
+ */
+
+
+// On Sun, to prevent conflict with ObjectSpace, G4Timer.hh has to be
+// loaded *before* globals.hh...
+#include "G4Timer.hh"
+
+#include <g4-tankresponse/g4classes/G4BeamTestRunManager.h>
+#include <G4ParticleGun.hh>
+#include <G4Run.hh>
+
+
+G4BeamTestRunManager::G4BeamTestRunManager(): G4RunManager()
+{
+
+}
+
+
+void G4BeamTestRunManager::BeamOn(G4int n_event,const char* macroFile,G4int n_select)
+{
+ G4String text = "BeamOn method is not supported in G4BeamTestRunManager!";
+ G4Exception("G4BeamTestRunManager::BeamOn()", "G4BeamTestRunManager001", FatalException, text);
+}
+
+
+void G4BeamTestRunManager::InitializeRun()
+{
+ G4bool cond = ConfirmBeamOnCondition();
+ if(cond)
+ {
+ // Reset the event counter
+ numberOfEventToBeProcessed = 0;
+ ConstructScoringWorlds();
+ RunInitialization();
+
+ if(verboseLevel>0) timer->Start();
+ }
+}
+
+
+void G4BeamTestRunManager::InjectParticle(G4ParticleGun* particleGun)
+{
+ if(!currentRun)
+ {
+ G4String text = "Run needs to be initialized before injecting a particle.";
+ G4Exception("G4BeamTestRunManager::InjectParticle()", "G4BeamTestRunManager002", FatalException, text);
+ }
+ assert(currentRun); // the G4Exception() above calls abort(). This assert() silences the clang static analyzer
+
+ numberOfEventToBeProcessed++;
+ currentRun->SetNumberOfEventToBeProcessed(numberOfEventToBeProcessed);
+
+ currentEvent = GenerateEvent(numberOfEventToBeProcessed);
+ particleGun->GeneratePrimaryVertex(currentEvent);
+
+ eventManager->ProcessOneEvent(currentEvent);
+ AnalyzeEvent(currentEvent);
+ Update_Scoring();
+ StackPreviousEvent(currentEvent);
+ currentEvent = 0;
+
+ if(runAborted) TerminateRun();
+}
+
+
+G4Event* G4BeamTestRunManager::GenerateEvent(G4int i_event)
+{
+ G4Event* anEvent = new G4Event(i_event);
+
+ if(storeRandomNumberStatusToG4Event==1 || storeRandomNumberStatusToG4Event==3)
+ {
+ std::ostringstream oss;
+ CLHEP::HepRandom::saveFullState(oss);
+ randomNumberStatusForThisEvent = oss.str();
+ anEvent->SetRandomNumberStatus(randomNumberStatusForThisEvent);
+ }
+
+ if(storeRandomNumberStatus)
+ {
+ G4String fileN = randomNumberStatusDir + "currentEvent.rndm";
+ CLHEP::HepRandom::saveEngineStatus(fileN);
+ }
+
+ return anEvent;
+}
+
+
+void G4BeamTestRunManager::TerminateRun()
+{
+ if(verboseLevel>0)
+ {
+ timer->Stop();
+ G4cout << "Run terminated." << G4endl;
+ G4cout << "Run Summary" << G4endl;
+ if(runAborted)
+ {
+ G4cout << " Run Aborted after " << numberOfEventToBeProcessed << " events processed." << G4endl;
+ }
+ else
+ {
+ G4cout << " Number of events processed : " << numberOfEventToBeProcessed << G4endl;
+ }
+ G4cout << " " << *timer << G4endl;
+ }
+
+ RunTermination();
+}
+
+
+//
+// The following method is an exact copy of
+// UpdateScoring which is private in the G4RunManager
+//
+
+#include <G4ScoringManager.hh>
+#include <G4HCofThisEvent.hh>
+#include <G4VHitsCollection.hh>
+
+void G4BeamTestRunManager::Update_Scoring()
+{
+ G4ScoringManager* ScM = G4ScoringManager::GetScoringManagerIfExist();
+ if(!ScM) return;
+ G4int nPar = ScM->GetNumberOfMesh();
+ if(nPar<1) return;
+
+ G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
+ if(!HCE) return;
+ G4int nColl = HCE->GetCapacity();
+ for(G4int i=0;i<nColl;i++)
+ {
+ G4VHitsCollection* HC = HCE->GetHC(i);
+ if(HC) ScM->Accumulate(HC);
+ }
+}
diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
new file mode 100644
index 0000000..4f12c7a
--- /dev/null
+++ b/src/G4BeamTestTank.cxx
@@ -0,0 +1,382 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4BeamTestTank.cxx 158930 2017-10-20 01:17:41Z cweaver $
+ *
+ * @file G4BeamTestTank.cxx
+ * @version $Rev: 158930 $
+ * @date $Date: 2017-10-20 02:17:41 +0100 (Fri, 20 Oct 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestTank.h>
+#include <g4-tankresponse/g4classes/G4TankIceSD.h>
+#include <icetray/I3Units.h>
+
+#include <icetray/OMKey.h>
+#include <dataclasses/TankKey.h>
+#include <dataclasses/geometry/I3TankGeo.h>
+#include <dataclasses/geometry/I3Geometry.h>
+#include <dataclasses/I3Position.h>
+
+#include <G4LogicalVolume.hh>
+#include <G4PVPlacement.hh>
+#include <G4Material.hh>
+#include <G4Tubs.hh>
+#include <G4Sphere.hh>
+#include <G4Box.hh>
+#include <G4SDManager.hh>
+#include <G4ThreeVector.hh>
+
+#include <G4VisAttributes.hh>
+
+#include <G4UserLimits.hh>
+
+#include <boost/foreach.hpp>
+
+//prevent gcc to make something stupid with pretended unused variables
+#ifdef __GNUC__
+#define ATTRIBUTE_UNUSED __attribute__((unused))
+#else
+#define ATTRIBUTE_UNUSED
+#endif
+
+
+G4BeamTestTank::G4BeamTestTank(const TankKey& tankKey, const I3Geometry& geometry):
+tankKey_(tankKey), geometry_(geometry)
+{
+ const I3StationGeoMap& stationMap = geometry.stationgeo;
+ unsigned int tankID = tankKey.tank==TankKey::TankA?0:1;
+ I3StationGeoMap::const_iterator station_iter = stationMap.find(tankKey.string);
+
+ if(station_iter==stationMap.end())
+ {
+ log_fatal("The requested station %d in not in the geometry!", tankKey.string);
+ return;
+ }
+
+ if(station_iter->second.size()<tankID)
+ {
+ log_fatal("The number of tanks in station %d is not correct!", tankKey.string);
+ return;
+ }
+ const I3TankGeo& tankGeo = station_iter->second.at(tankID);
+
+ // Get tank dimensions
+ tankThickness_ = 0.5*CLHEP::cm;
+ tankHeight_ = (tankGeo.tankheight / I3Units::m) * CLHEP::m + tankThickness_;
+ innerRadius_ = (tankGeo.tankradius / I3Units::m) * CLHEP::m;
+ outerRadius_ = innerRadius_ + tankThickness_;
+
+ // Get fill and snow heights
+ fillHeight_ = (tankGeo.fillheight / I3Units::m) * CLHEP::m;
+ snowHeight_ = (tankGeo.snowheight / I3Units::m) * CLHEP::m;
+ perliteHeight_ = tankHeight_ - tankThickness_ - fillHeight_;
+
+ // Set DOM dimensions
+ glassOuterRadius_ = 6.5 * 2.54 * CLHEP::cm; // 6.5" outer glass sphere radius
+ glassThickness_ = 0.5 * 2.54 * CLHEP::cm; // 0.5" glass sphere thickness
+
+ // Calculate tank position (tank center)
+ // tankGeo.position corresponds to the average position of the two DOMs in a tank
+ position_.set((tankGeo.position.GetX() / I3Units::m) * CLHEP::m,
+ (tankGeo.position.GetY() / I3Units::m) * CLHEP::m,
+ (tankGeo.position.GetZ() / I3Units::m) * CLHEP::m - fillHeight_ + 0.5*tankHeight_);
+
+ // Get positions of the doms relativ to tank center
+ BOOST_FOREACH(const OMKey& omKey, tankGeo.omKeyList_)
+ {
+ I3OMGeoMap::const_iterator omGeo_iter = geometry_.omgeo.find(omKey);
+ if(omGeo_iter==geometry_.omgeo.end())
+ {
+ log_error_stream(omKey << " is missing in Tank " << tankKey_);
+ continue;
+ }
+
+ G4ThreeVector relPos(omGeo_iter->second.position.GetX() - tankGeo.position.GetX(),
+ omGeo_iter->second.position.GetY() - tankGeo.position.GetY(),
+ omGeo_iter->second.position.GetZ() - tankGeo.position.GetZ());
+
+ relDomPositions_[omKey] = (relPos / I3Units::m) * CLHEP::m;
+ }
+
+ //
+ // Calculate Delaunay points
+ //
+ G4ThreeVector triangleDir(NAN, NAN, NAN);
+ switch(station_iter->second.size())
+ {
+ case 1: // Single tank
+ {
+ // Vector orthogonal to DOM positions
+ triangleDir.set(relDomPositions_.begin()->second.y(),
+ -relDomPositions_.begin()->second.x(),
+ 0.0);
+ break;
+ }
+ case 2: // Two tanks
+ {
+ const I3TankGeo& neighborGeo = station_iter->second.at(tankID==0?1:0);
+ G4ThreeVector neighborPos(neighborGeo.position.GetX(),
+ neighborGeo.position.GetY(),
+ neighborGeo.position.GetZ());
+
+ // Convert to G4 units
+ neighborPos *= CLHEP::m/I3Units::m;
+
+ // Same z position same as other tank
+ neighborPos.setZ(position_.z());
+
+ triangleDir = position_ - neighborPos;
+ break;
+ }
+ default:
+ {
+ log_fatal("Invalid number of tanks (%zu) in station %d!",
+ station_iter->second.size(), tankKey_.string);
+ break;
+ }
+ }
+
+ // side length
+ double triangleLength = 10.0 * CLHEP::m;
+
+ triangleDir.setMag(0.5*triangleLength/cos(CLHEP::pi/6.0));
+ delaunayPoint1_ = position_ + triangleDir;
+
+ // Rotate by 120 deg
+ triangleDir.rotateZ(CLHEP::pi/1.5);
+ delaunayPoint2_ = position_ + triangleDir;
+
+ // Rotate by 120 deg
+ triangleDir.rotateZ(CLHEP::pi/1.5);
+ delaunayPoint3_ = position_ + triangleDir;
+}
+
+
+G4BeamTestTank::~G4BeamTestTank()
+{
+
+}
+
+
+G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const G4ThreeVector& origin)
+{
+ // User limits (energy cutoffs)
+ // Do not create photons or electrons below cherenkov threshold
+ // See also corresponding UserSpecialCuts in Physicslist !!!!
+ // Maybe do all of this as stepping action ??????
+ G4UserLimits* energyLimit = new G4UserLimits();
+ energyLimit->SetUserMinEkine(280.0 * CLHEP::keV); // Cherenkov threshold of electrons in ice
+
+ std::string tankName=boost::lexical_cast<std::string>(tankKey_);
+
+ // Define plastic frame
+ G4Material* plastic = G4Material::GetMaterial("Plastic");
+ G4Tubs* solidTank = new G4Tubs(("solid_tank_" + tankName).c_str(),
+ 0.0 * CLHEP::m, outerRadius_, 0.5 * tankHeight_,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+ tankLog_ = new G4LogicalVolume(solidTank, plastic,
+ ("log_tank_" + tankName).c_str(), 0, 0, 0);
+
+ // Define ice volume
+ G4Material* ice = G4Material::GetMaterial("Ice");
+ G4Tubs* solidIce = new G4Tubs(("solid_ice_" + tankName).c_str(),
+ 0.0 * CLHEP::m, innerRadius_, 0.5 * fillHeight_,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+ G4LogicalVolume* logIce =
+ new G4LogicalVolume(solidIce, ice, ("log_ice_" + tankName).c_str(), 0, 0, 0);
+ G4ThreeVector physIcePosition(0, 0, -0.5*tankHeight_ + tankThickness_ + 0.5*fillHeight_);
+ G4VPhysicalVolume* physIce ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, physIcePosition, logIce,
+ ("ice_" + tankName).c_str(), tankLog_, false, 0);
+
+ // Define perlite volume
+ G4Material* perlite = G4Material::GetMaterial("Perlite");
+ G4Tubs* solidPerlite = new G4Tubs(("solid_perlite_" + tankName).c_str(),
+ 0.0 * CLHEP::m, innerRadius_, 0.5 * perliteHeight_,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg);
+ G4LogicalVolume* logPerlite =
+ new G4LogicalVolume(solidPerlite, perlite, ("log_perlite_" + tankName).c_str(), 0, 0, 0);
+ G4ThreeVector physPerlitePosition(0, 0, -0.5 * tankHeight_ + 0.5 * CLHEP::cm + fillHeight_ +
+ 0.5 * perliteHeight_);
+ G4VPhysicalVolume* physPerlite ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, physPerlitePosition, logPerlite,
+ ("perlite_" + tankName).c_str(), tankLog_, false, 0);
+
+ // Define glass sphere & effective DOM material splitted in upper and lower part
+ G4Material* glass = G4Material::GetMaterial("Glass");
+ G4Material* effectiveDOM = G4Material::GetMaterial("effectiveDOM");
+
+ std::map<OMKey, G4ThreeVector> domPosIce;
+ std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+ for(om_iter=relDomPositions_.begin(); om_iter!=relDomPositions_.end(); ++om_iter)
+ {
+ const OMKey& omKey = om_iter->first;
+
+ G4ThreeVector upperDOMpos(om_iter->second.x(), om_iter->second.y(), -0.5 * perliteHeight_);
+ G4ThreeVector lowerDOMpos(om_iter->second.x(), om_iter->second.y(), 0.5 * fillHeight_);
+
+ domPosIce[omKey] = lowerDOMpos;
+
+ std::string omName=boost::lexical_cast<std::string>(omKey);
+
+ G4Sphere *upperglasssphere = new G4Sphere (("solid_dom_up_" + omName).c_str(),
+ 0.0 * CLHEP::m, glassOuterRadius_,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+ 0.0 * CLHEP::deg, 90.0 * CLHEP::deg);
+ G4Sphere *lowerglasssphere = new G4Sphere (("solid_dom_lo_" + omName).c_str(),
+ 0.0 * CLHEP::m, glassOuterRadius_,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+ 90.0 * CLHEP::deg, 180.0 * CLHEP::deg);
+
+ G4double domInnerRadius = glassOuterRadius_ - glassThickness_;
+ G4Sphere *upperdomsphere = new G4Sphere (("solid_inside_dom_up_" + omName).c_str(),
+ 0.0 * CLHEP::m, domInnerRadius,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+ 0.0 * CLHEP::deg, 90.0 * CLHEP::deg);
+ G4Sphere *lowerdomsphere = new G4Sphere (("solid_inside_dom_lo_" + omName).c_str(),
+ 0.0 * CLHEP::m, domInnerRadius,
+ 0.0 * CLHEP::deg, 360.0 * CLHEP::deg,
+ 90.0 * CLHEP::deg, 180.0 * CLHEP::deg);
+
+ G4LogicalVolume* logUpperGlass =
+ new G4LogicalVolume(upperglasssphere, glass,
+ ("log_dom_up_" + omName).c_str(), 0, 0, 0);
+ G4LogicalVolume* logLowerGlass =
+ new G4LogicalVolume(lowerglasssphere, glass,
+ ("log_dom_lo_" + omName).c_str(), 0, 0, 0);
+ G4LogicalVolume* logUpperDOM =
+ new G4LogicalVolume(upperdomsphere, effectiveDOM,
+ ("log_inside_dom_up_" + omName).c_str(), 0, 0, 0);
+ G4LogicalVolume* logLowerDOM =
+ new G4LogicalVolume(lowerdomsphere, effectiveDOM,
+ ("log_inside_dom_lo_" + omName).c_str(), 0, 0, 0);
+ G4VPhysicalVolume* physUpperGlass ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, upperDOMpos, logUpperGlass,
+ ("dom_up_" + omName).c_str(), logPerlite, false, 0);
+ G4VPhysicalVolume* physLowerGlass ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, lowerDOMpos, logLowerGlass,
+ ("dom_lo_" + omName).c_str(), logIce, false, 0);
+ G4VPhysicalVolume* physUpperDOM ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, G4ThreeVector(0,0,0), logUpperDOM,
+ ("inside_dom_up_" + omName).c_str(), logUpperGlass, false, 0);
+ G4VPhysicalVolume* physLowerDOM ATTRIBUTE_UNUSED =
+ new G4PVPlacement(0, G4ThreeVector(0,0,0), logLowerDOM,
+ ("inside_dom_lo_" + omName).c_str(), logLowerGlass, false, 0);
+
+ // apply energy limits
+ logUpperGlass->SetUserLimits(energyLimit);
+ logLowerGlass->SetUserLimits(energyLimit);
+ logUpperDOM->SetUserLimits(energyLimit);
+ logLowerDOM->SetUserLimits(energyLimit);
+ }
+
+ // Define sensitive detector
+ G4SDManager* sdManager = G4SDManager::GetSDMpointer();
+ iceSD_ = new G4TankIceSD(("ice_SD_" + tankName).c_str(), domPosIce);
+ sdManager->AddNewDetector(iceSD_);
+ logIce->SetSensitiveDetector(iceSD_);
+
+ // Instantiation of a set of visualization attributes with red colour
+ G4VisAttributes * tankVisAtt = new G4VisAttributes(G4Colour(1,0,0));
+ // Set the forced wireframe style
+ //snowVisAtt->SetForceWireFrame(true);
+ // Assignment of the visualization attributes to the logical volume
+ tankLog_->SetVisAttributes(tankVisAtt);
+
+ G4ThreeVector tankPos = position_ - origin - mother->GetTranslation();
+
+ G4VPhysicalVolume* tankPhys = new G4PVPlacement(0, tankPos, tankLog_,
+ ("tank_" + tankName).c_str(),
+ mother->GetLogicalVolume(), false, 0);
+
+ // apply energy limits
+ tankLog_->SetUserLimits(energyLimit);
+ logPerlite->SetUserLimits(energyLimit);
+ logIce->SetUserLimits(energyLimit);
+
+ return tankPhys;
+}
+
+
+double G4BeamTestTank::GetNumCherenkov(const OMKey& omKey)
+{
+ return std::max(iceSD_->GetNumCherenkov(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetNumCherenkovWeight(const OMKey& omKey)
+{
+ return std::max(iceSD_->GetNumCherenkovWeight(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetEDep_G4(const OMKey& omKey)
+{
+ return std::max(iceSD_->GetEDep(omKey), 0.);
+}
+
+
+double G4BeamTestTank::GetTime_G4(const OMKey& omKey)
+{
+ return iceSD_->GetTime(omKey);
+}
+
+
+double G4BeamTestTank::GetEDep_I3(const OMKey& omKey)
+{
+ return std::max(iceSD_->GetEDep(omKey), 0.) / CLHEP::keV * I3Units::keV;
+}
+
+
+double G4BeamTestTank::GetTime_I3(const OMKey& omKey)
+{
+ return (iceSD_->GetTime(omKey) / CLHEP::s) * I3Units::s;
+}
+
+
+I3Position G4BeamTestTank::GetPos_I3()
+{
+ I3Position pos((position_.x() / CLHEP::m) * I3Units::m,
+ (position_.y() / CLHEP::m) * I3Units::m,
+ (position_.z() / CLHEP::m) * I3Units::m);
+ return pos;
+}
+
+
+double G4BeamTestTank::GetX_I3()
+{
+ return (position_.x() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetY_I3()
+{
+ return (position_.y() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetZ_I3()
+{
+ return (position_.z() / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetTankHeight_I3()
+{
+ return (tankHeight_ / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetTankRadius_I3()
+{
+ return (outerRadius_ / CLHEP::m) * I3Units::m;
+}
+
+
+double G4BeamTestTank::GetSnowHeight_I3()
+{
+ return (snowHeight_ / CLHEP::m) * I3Units::m;
+}
diff --git a/src/G4BeamTestUserSteppingAction.cxx b/src/G4BeamTestUserSteppingAction.cxx
new file mode 100644
index 0000000..b16098c
--- /dev/null
+++ b/src/G4BeamTestUserSteppingAction.cxx
@@ -0,0 +1,42 @@
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserSteppingAction.cxx
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestUserSteppingAction.h>
+
+#include "G4Step.hh"
+#include "G4Track.hh"
+#include "G4LogicalVolume.hh"
+#include "G4UserLimits.hh"
+#include "G4SteppingManager.hh"
+
+G4BeamTestUserSteppingAction::G4BeamTestUserSteppingAction(){}
+
+void G4BeamTestUserSteppingAction::UserSteppingAction(const G4Step*)
+{
+ G4Track* track = fpSteppingManager->GetTrack();
+ if(track)
+ {
+ const G4LogicalVolume *volume = track->GetLogicalVolumeAtVertex();
+ G4UserLimits *limit = volume->GetUserLimits();
+ if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
+ G4double threshold = limit->GetUserMinEkine(*track);
+ //check if particle is a gamma
+ if(track->GetDefinition()->GetParticleName() == "gamma")
+ {
+ //check if particle energy is below threshold; if true, kill the particle
+ G4double energy = track->GetTotalEnergy();
+ if(energy < threshold) track->SetTrackStatus(fStopAndKill);
+ }
+ }
+}
diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx
new file mode 100644
index 0000000..d554152
--- /dev/null
+++ b/src/G4BeamTestUserTrackingAction.cxx
@@ -0,0 +1,50 @@
+/**
+ * Copyright (C) 2011
+ * The IceCube collaboration
+ * ID: $Id$
+ *
+ * @file G4BeamTestUserTrackingAction.cxx
+ * @version $Revision$
+ * @date $Date$
+ * @author Thomas Melzig
+ *
+ * $LastChangedBy$
+ */
+
+
+#include <g4-tankresponse/g4classes/G4BeamTestUserTrackingAction.h>
+
+#include "G4Track.hh"
+#include "G4UserLimits.hh"
+#include "G4TrackVector.hh"
+#include "G4TrackingManager.hh"
+
+G4BeamTestUserTrackingAction::G4BeamTestUserTrackingAction(){}
+
+void G4BeamTestUserTrackingAction::PreUserTrackingAction(const G4Track*){}
+
+void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track)
+{
+ const G4LogicalVolume *volume = track->GetLogicalVolumeAtVertex();
+ G4UserLimits *limit = volume->GetUserLimits();
+ if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
+ G4double threshold = limit->GetUserMinEkine(*track);
+ G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
+ if(secondaries)
+ {
+ size_t nSeco = secondaries->size();
+ if(nSeco>0)
+ {
+ for(size_t i=0;i<nSeco;i++)
+ {
+ //check if secondary particle is a gamma
+ if((*secondaries)[i]->GetDefinition()->GetParticleName() == "gamma")
+ {
+ //check if particle energy is below threshold; if true, kill the particle
+ G4double energy = (*secondaries)[i]->GetTotalEnergy();
+ if(energy < threshold) (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+ }
+ }
+ }
+ }
+}
diff --git a/src/G4Interface.cxx b/src/G4Interface.cxx
new file mode 100644
index 0000000..36ba3d9
--- /dev/null
+++ b/src/G4Interface.cxx
@@ -0,0 +1,292 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4Interface.cxx 149388 2016-08-18 21:50:04Z jgonzalez $
+ *
+ * @file G4Interface.cxx
+ * @version $Rev: 149388 $
+ * @date $Date: 2016-08-18 22:50:04 +0100 (Thu, 18 Aug 2016) $
+ * @author Tilo Waldenmaier
+ */
+
+#include <dataclasses/physics/I3Particle.h>
+
+#include <g4-tankresponse/g4classes/G4Interface.h>
+#include <g4-tankresponse/g4classes/G4BeamTestDetectorConstruction.h>
+#include <g4-tankresponse/g4classes/G4BeamTestPhysicsList.h>
+#include <g4-tankresponse/g4classes/G4BeamTestUserTrackingAction.h>
+#include <g4-tankresponse/g4classes/G4BeamTestUserSteppingAction.h>
+
+#include <icetray/I3Logging.h>
+
+#ifdef G4VIS_USE
+#include <G4VisExecutive.hh>
+#endif
+
+#include <G4ParticleGun.hh>
+#include <G4ParticleTable.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4UImanager.hh>
+
+
+G4Interface* G4Interface::g4Interface_ = NULL;
+
+G4Interface::G4Interface(const std::string& visMacro):
+ detector_(NULL), initialized_(false),
+ eventInitialized_(false), visMacro_(visMacro)
+{
+ g4Interface_ = this;
+
+ // Visualization manager
+#ifdef G4VIS_USE
+ visManager_ = NULL;
+ if(!visMacro_.empty())
+ {
+ visManager_ = new G4VisExecutive();
+ visManager_->Initialize();
+ }
+#endif
+}
+
+
+G4Interface::~G4Interface()
+{
+ g4Interface_ = NULL;
+
+#ifdef G4VIS_USE
+ if(visManager_) delete visManager_;
+#endif
+}
+
+
+void G4Interface::InstallTank(G4BeamTestTank* tank)
+{
+ if(initialized_)
+ {
+ log_fatal("G4Interface aleady initialized. Cannot install tank!");
+ return;
+ }
+
+ if(!detector_) detector_ = new G4BeamTestDetectorConstruction();
+ detector_->InstallTank(tank);
+}
+
+
+void G4Interface::InitializeEvent()
+{
+ ///
+ /// An IceTray EVENT corresponds to a G4RUN
+ /// where each injected particle initiates an G4EVENT
+ ///
+
+ if(!initialized_)
+ {
+ Initialize();
+ }
+
+ if(!eventInitialized_)
+ {
+ runManager_.InitializeRun();
+ eventInitialized_ = true;
+ }
+}
+
+
+void G4Interface::InjectParticle(const I3Particle& particle)
+{
+ if(!eventInitialized_)
+ {
+ log_fatal("No event initialized. Cannot inject particle!");
+ return;
+ }
+
+ G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+ G4ParticleDefinition* particleDef = NULL;
+ switch(particle.GetType())
+ {
+ case I3Particle::Gamma:
+ particleDef = particleTable->FindParticle("gamma");
+ break;
+ case I3Particle::EMinus:
+ particleDef = particleTable->FindParticle("e-");
+ break;
+ case I3Particle::EPlus:
+ particleDef = particleTable->FindParticle("e+");
+ break;
+ case I3Particle::MuMinus:
+ particleDef = particleTable->FindParticle("mu-");
+ break;
+ case I3Particle::MuPlus:
+ particleDef = particleTable->FindParticle("mu+");
+ break;
+ case I3Particle::PPlus:
+ particleDef = particleTable->FindParticle("proton");
+ break;
+ case I3Particle::PMinus:
+ particleDef = particleTable->FindParticle("anti_proton");
+ break;
+ case I3Particle::Neutron:
+ particleDef = particleTable->FindParticle("neutron");
+ break;
+#ifdef I3PARTICLE_SUPPORTS_PDG_ENCODINGS
+ case I3Particle::NeutronBar:
+#else
+ case 25:
+#endif
+ particleDef = particleTable->FindParticle("anti_neutron");
+ break;
+ case I3Particle::PiPlus:
+ particleDef = particleTable->FindParticle("pi+");
+ break;
+ case I3Particle::PiMinus:
+ particleDef = particleTable->FindParticle("pi-");
+ break;
+ case I3Particle::Pi0:
+ particleDef = particleTable->FindParticle("pi0");
+ break;
+ case I3Particle::KPlus:
+ particleDef = particleTable->FindParticle("kaon+");
+ break;
+ case I3Particle::KMinus:
+ particleDef = particleTable->FindParticle("kaon-");
+ break;
+ case I3Particle::K0_Long:
+ particleDef = particleTable->FindParticle("kaon0L");
+ break;
+ case I3Particle::K0_Short:
+ particleDef = particleTable->FindParticle("kaon0S");
+ break;
+ case I3Particle::NuE:
+ particleDef = particleTable->FindParticle("nu_e");
+ break;
+ case I3Particle::NuEBar:
+ particleDef = particleTable->FindParticle("anti_nu_e");
+ break;
+ case I3Particle::NuMu:
+ particleDef = particleTable->FindParticle("nu_mu");
+ break;
+ case I3Particle::NuMuBar:
+ particleDef = particleTable->FindParticle("anti_nu_mu");
+ break;
+ case I3Particle::NuTau:
+ particleDef = particleTable->FindParticle("nu_tau");
+ break;
+ case I3Particle::NuTauBar:
+ particleDef = particleTable->FindParticle("anti_nu_tau");
+ break;
+ default:
+ log_warn("Man, check out that strange particle \"%s\" ?!", particle.GetTypeString().c_str());
+ return;
+ }
+
+ // Particle position in G4 units
+ G4ThreeVector position((particle.GetX() / I3Units::m) * CLHEP::m,
+ (particle.GetY() / I3Units::m) * CLHEP::m,
+ (particle.GetZ() / I3Units::m) * CLHEP::m);
+
+ // Transform I3 coorinates to world system
+ position -= detector_->GetWorldOrigin();
+
+ G4ThreeVector direction(particle.GetDir().GetX(),
+ particle.GetDir().GetY(),
+ particle.GetDir().GetZ());
+
+ G4ParticleGun gun(1);
+ gun.SetParticleDefinition(particleDef);
+ gun.SetParticleEnergy((particle.GetEnergy() / I3Units::GeV) * CLHEP::GeV);
+ gun.SetParticlePosition(position);
+ gun.SetParticleMomentumDirection(direction);
+
+ log_trace("Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV",
+ particle.GetTypeString().c_str(),
+ position.x() / CLHEP::m,
+ position.y() / CLHEP::m,
+ position.z() / CLHEP::m,
+ gun.GetParticleEnergy() / CLHEP::MeV);
+
+ runManager_.InjectParticle(&gun);
+}
+
+
+void G4Interface::TerminateEvent()
+{
+ ///
+ /// An IceTray EVENT corresponds to a G4RUN
+ /// where each injected particle initiates an G4EVENT
+ ///
+
+ if(eventInitialized_)
+ {
+ runManager_.TerminateRun();
+ eventInitialized_ = false;
+ }
+}
+
+
+void G4Interface::Initialize()
+{
+ if(initialized_)
+ {
+ log_error("G4Interface has already been initialized. Ignoring this call!");
+ return;
+ }
+
+ log_debug("Init geometry ...");
+ runManager_.SetUserInitialization(detector_);
+
+ log_debug("Init physics list ...");
+ runManager_.SetUserInitialization(new G4BeamTestPhysicsList());
+
+ log_debug("Init UserTrackingAction ...");
+ runManager_.SetUserAction(new G4BeamTestUserTrackingAction());
+
+ log_debug("Init UserSteppingAction ...");
+ runManager_.SetUserAction(new G4BeamTestUserSteppingAction());
+
+ // Initialize G4 kernel
+ log_debug("Init run manager ...");
+ runManager_.Initialize();
+
+ // Set verbosity
+ int verboseLevel = 0;
+ switch (GetIcetrayLogger()->LogLevelForUnit("G4Interface"))
+ {
+ case I3LOG_FATAL:
+ case I3LOG_ERROR:
+ case I3LOG_WARN:
+ case I3LOG_INFO:
+ case I3LOG_NOTICE:
+ default:
+ verboseLevel = 0;
+ break;
+ case I3LOG_DEBUG:
+ verboseLevel = 1;
+ break;
+ case I3LOG_TRACE:
+ verboseLevel = 2;
+ break;
+ }
+
+ runManager_.SetVerboseLevel(verboseLevel);
+ G4EventManager::GetEventManager()->SetVerboseLevel(verboseLevel);
+ G4EventManager::GetEventManager()->GetStackManager()->SetVerboseLevel(verboseLevel);
+ G4EventManager::GetEventManager()->GetTrackingManager()->SetVerboseLevel(verboseLevel);
+#ifdef G4VIS_USE
+ if(visManager_) visManager_->SetVerboseLevel(verboseLevel);
+#endif
+
+ // Execute visualization macro (if specified)
+ if(!visMacro_.empty())
+ {
+ G4UImanager* uim = G4UImanager::GetUIpointer();
+
+ // Checking geometry
+ uim->ApplyCommand("/geometry/test/grid_test");
+
+ // Execute visualization macro
+ std::string visCmd = "/control/execute " + visMacro_;
+ uim->ApplyCommand(visCmd.c_str());
+ }
+
+ initialized_ = true;
+}
diff --git a/src/G4TankIceSD.cxx b/src/G4TankIceSD.cxx
new file mode 100644
index 0000000..4a1be0b
--- /dev/null
+++ b/src/G4TankIceSD.cxx
@@ -0,0 +1,176 @@
+/**
+ * Copyright (C) 2009
+ * The IceCube collaboration
+ * ID: $Id: G4TankIceSD.cxx 152814 2017-01-19 21:34:52Z jgonzalez $
+ *
+ * @file G4TankIceSD.cxx
+ * @version $Rev: 152814 $
+ * @date $Date: 2017-01-19 21:34:52 +0000 (Thu, 19 Jan 2017) $
+ * @author Tilo Waldenmaier, Thomas Melzig
+ */
+
+
+#include <g4-tankresponse/g4classes/G4TankIceSD.h>
+#include <g4-tankresponse/g4classes/G4BeamTestTank.h>
+
+#include <G4Step.hh>
+#include <G4HCofThisEvent.hh>
+#include <G4TouchableHistory.hh>
+#include <G4ios.hh>
+#include <iostream>
+#include <ios>
+#include <fstream>
+
+#include <G4VProcess.hh>
+#include <G4OpticalPhoton.hh>
+#include "G4Poisson.hh"
+
+#include <icetray/I3Units.h>
+
+#include "G4Version.hh"
+
+#if G4VERSION_NUMBER >= 950
+// The G4MaterialPropertyVector is gone since 4.9.5.
+// It has been typedef'd to a G4UnorderedPhysicsVector
+// with a different interface. Try to support both
+// versions with an ifdef.
+#define MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+#endif
+
+
+G4TankIceSD::G4TankIceSD(G4String name, const std::map<OMKey, G4ThreeVector>& domPositions)
+ : G4VSensitiveDetector(name), domPositions_(domPositions)
+{}
+
+
+G4TankIceSD::~G4TankIceSD() {}
+
+
+void G4TankIceSD::Initialize(G4HCofThisEvent* HCE)
+{
+ std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+ for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+ {
+ sumEdep_[om_iter->first] = 0;
+ cogTime_[om_iter->first] = 0;
+ cherenkovCounter_[om_iter->first] = 0;
+ cherenkovCounterWeight_[om_iter->first] = 0;
+ }
+}
+
+
+G4bool G4TankIceSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // Get energy deposit/time --> CHECK WHAT THESE QUANTITIES ACTUALLY MEAN !!!!!!!!!!!!!!!!!!!!!!!!
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
+ G4double cherenkovNumber = GetCerenkovNumber(aStep);
+
+ if(edep<=0 && cherenkovNumber<=0) return false;
+
+ // Get Position relative to ice center
+ G4ThreeVector preStepPoint = aStep->GetPreStepPoint()->GetPosition();
+ G4ThreeVector postStepPoint = aStep->GetPostStepPoint()->GetPosition();
+ G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle();
+ G4ThreeVector worldPosition = (preStepPoint + postStepPoint) / 2.0;
+ G4ThreeVector localPosition = theTouchable->GetHistory()->
+ GetTopTransform().TransformPoint(worldPosition);
+
+ std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+ for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+ {
+ G4double radius = sqrt(pow(om_iter->second.x() - localPosition.x(), 2) +
+ pow(om_iter->second.y() - localPosition.y(), 2));
+ G4double height = (om_iter->second.z() - localPosition.z());
+
+ sumEdep_[om_iter->first] += edep;
+ cogTime_[om_iter->first] += edep*time;
+ cherenkovCounterWeight_[om_iter->first] += GetProbability(radius, height) * cherenkovNumber;
+ cherenkovCounter_[om_iter->first] += cherenkovNumber;
+ }
+ return true;
+}
+
+
+void G4TankIceSD::EndOfEvent(G4HCofThisEvent*)
+{
+ std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
+ for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
+ {
+ if(sumEdep_[om_iter->first]>0)
+ {
+ cogTime_[om_iter->first] /= sumEdep_[om_iter->first];
+ }
+ }
+}
+
+
+G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep)
+{
+ // same function as in "source/processes/electromagnetic/xrays/src/G4Cerenkov.cc" but without
+ // cerenkov angle and only for an energy independent refraction index
+ G4Track* aTrack = aStep->GetTrack();
+ const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
+ const G4Material* aMaterial = aTrack->GetMaterial();
+
+ G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
+ G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
+
+ G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
+ if (!aMaterialPropertiesTable) return 0.0;
+
+ G4MaterialPropertyVector* Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
+ if (!Rindex) return 0.0;
+
+ const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm);
+ const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
+ const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
+ if (beta <= 0.0) return 0.0;
+ G4double BetaInverse = 1. / beta;
+
+ // Min and Max photon energies
+#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+ G4double Pmin = Rindex->GetMinLowEdgeEnergy();
+ G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
+#else
+ G4double Pmin = Rindex->GetMinPhotonEnergy();
+ G4double Pmax = Rindex->GetMaxPhotonEnergy();
+#endif
+
+ // Min and Max Refraction Indices
+#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
+ G4double nMin = Rindex->GetMinValue();
+ G4double nMax = Rindex->GetMaxValue();
+#else
+ G4double nMin = Rindex->GetMinProperty();
+ G4double nMax = Rindex->GetMaxProperty();
+#endif
+ if (nMin!=nMax)
+ {
+ log_error("Support for energy dependent refraction index not yet implemented!");
+ return 0.0;
+ }
+
+ // If n(Pmax) < 1/Beta -- no photons generated
+ if (nMax < BetaInverse) return 0.0;
+
+ G4double MeanNumberOfPhotons = Rfact * charge/CLHEP::eplus * charge/CLHEP::eplus * (Pmax - Pmin) *
+ (1. - BetaInverse * BetaInverse / nMin / nMin);
+
+ G4double step_length = aStep->GetStepLength();
+
+ return MeanNumberOfPhotons * step_length;
+}
+
+
+G4double G4TankIceSD::GetProbability(G4double radius, G4double height)
+{
+ height = 0.90 - height / CLHEP::m;
+ radius = radius / CLHEP::m;
+ G4double p0 = 0.0340 + 0.0120 * height;
+ G4double p1 = 0.0400 + 0.000622 * exp(-1.13 + 8.03 * height);
+ G4double p2 = 0.917 + 2.34 * exp(-1.82 + 3.12 * height);
+ G4double p3 = -0.00325 - 0.00199 * height - 0.0121 * height*height;
+
+ return (p0 + p3 * radius) + (p1 - p0) * exp(-p2*p2 * radius*radius);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 15:43:10 +0000