From 2925dbae7f488c3113b5f574eca1ceaba0ffdaac Mon Sep 17 00:00:00 2001 From: shivesh Date: Tue, 21 Aug 2018 15:36:34 +0100 Subject: Tue 21 Aug 15:36:34 BST 2018 --- src/G4BeamTestDetectorConstruction.cxx | 93 +++---- src/G4BeamTestEMPhysics.cxx | 13 - src/G4BeamTestGeneralPhysics.cxx | 11 - src/G4BeamTestHadronPhysics.cxx | 425 ----------------------------- src/G4BeamTestHadronPhysics.cxx.backup | 412 +++++++++++++++++++++++++++++ src/G4BeamTestIonPhysics.cxx | 112 -------- src/G4BeamTestIonPhysics.cxx.backup | 100 +++++++ src/G4BeamTestMuonPhysics.cxx | 15 +- src/G4BeamTestPhysicsList.cxx | 12 - src/G4BeamTestRunManager.cxx | 14 +- src/G4BeamTestTank.cxx | 470 +++++++++++++-------------------- src/G4BeamTestUserSteppingAction.cxx | 16 +- src/G4BeamTestUserTrackingAction.cxx | 16 +- src/G4Interface.cxx | 244 ++++++----------- src/G4TankIceSD.cxx | 80 +++--- 15 files changed, 840 insertions(+), 1193 deletions(-) delete mode 100644 src/G4BeamTestHadronPhysics.cxx create mode 100644 src/G4BeamTestHadronPhysics.cxx.backup delete mode 100644 src/G4BeamTestIonPhysics.cxx create mode 100644 src/G4BeamTestIonPhysics.cxx.backup (limited to 'src') diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx index 6eeaa54..cf75551 100644 --- a/src/G4BeamTestDetectorConstruction.cxx +++ b/src/G4BeamTestDetectorConstruction.cxx @@ -30,71 +30,31 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct() 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); + /* // World origin in IceCube coordinates */ + /* origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength); */ // 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); + G4double xWorld = 20.0 * CLHEP::m; + G4double yWorld = 20.0 * CLHEP::m; + G4double zWorld = 20.0 * CLHEP::m; // Create world volume - G4Box* world_box = new G4Box("solid_world", xHalfLength, yHalfLength, zHalfLength); + G4Box* world_box = new G4Box("solid_world", xWorld, yWorld, zWorld); G4LogicalVolume* worldLog = - new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0); + new G4LogicalVolume(world_box, G4Material::GetMaterial("G4_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); + /* // 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_); - } + // Install tank + /* BOOST_FOREACH(G4BeamTestTank* tank, tankList_) */ + /* { */ + tank_->InstallTank(worldPhys, origin_); + // } // User limits (energy cutoffs) // Do not create photons or electrons below cherenkov threshold @@ -102,7 +62,7 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct() G4UserLimits* energyLimit = new G4UserLimits(); energyLimit->SetUserMinEkine(280.0 * CLHEP::keV); // Cherenkov threshold of electrons in ice TODO(shivesh) worldLog->SetUserLimits(energyLimit); - snowLog->SetUserLimits(energyLimit); + /* snowLog->SetUserLimits(energyLimit); */ return worldPhys; } @@ -128,7 +88,7 @@ void G4BeamTestDetectorConstruction::CreateMaterials() void G4BeamTestDetectorConstruction::CreateAir() { G4NistManager* nistManager = G4NistManager::Instance(); - nistManager->ConstructNewGasMaterial("Air","G4_AIR"/* , (273.15 - 40.0)*CLHEP::kelvin, 670.0E-3*CLHEP::bar */); + nistManager->ConstructNewGasMaterial("Air","G4_AIR" , (273.15 + 20.0)*CLHEP::kelvin, 1.0*CLHEP::atmosphere); } /*****************************************************************/ @@ -136,10 +96,19 @@ void G4BeamTestDetectorConstruction::CreateAir() 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"); + G4Material* water = new G4Material("Water", 1.0 * CLHEP::g / CLHEP::cm3, 2, kStateLiquid); + water->AddElement(nistManager->FindOrBuildElement("H"), 2); + water->AddElement(nistManager->FindOrBuildElement("O"), 1); + + // pmt spectral response 300-650nm + // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf + // TODO(shivesh): add more properties? + const G4int water_bins = 2; + G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV}; + G4double water_refr[water_bins] = {1.33, 1.33}; + G4MaterialPropertiesTable *mpt_water = new G4MaterialPropertiesTable (); + mpt_water->AddProperty ("RINDEX", water_ephot, water_refr, water_bins); + water->SetMaterialPropertiesTable(mpt_water); } /*****************************************************************/ diff --git a/src/G4BeamTestEMPhysics.cxx b/src/G4BeamTestEMPhysics.cxx index d36207a..55ecf99 100644 --- a/src/G4BeamTestEMPhysics.cxx +++ b/src/G4BeamTestEMPhysics.cxx @@ -1,16 +1,3 @@ -/* - * 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 , Peter Nießen (tanktop) - * Last changed by: $LastChangedBy: kislat $ - */ - - #include #include #include diff --git a/src/G4BeamTestGeneralPhysics.cxx b/src/G4BeamTestGeneralPhysics.cxx index 2def9b4..93304a8 100644 --- a/src/G4BeamTestGeneralPhysics.cxx +++ b/src/G4BeamTestGeneralPhysics.cxx @@ -1,14 +1,3 @@ -/* - * 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 Last changed by: $LastChangedBy: jgonzalez $ - */ - #include "G4BeamTestGeneralPhysics.h" #include diff --git a/src/G4BeamTestHadronPhysics.cxx b/src/G4BeamTestHadronPhysics.cxx deleted file mode 100644 index 5150b67..0000000 --- a/src/G4BeamTestHadronPhysics.cxx +++ /dev/null @@ -1,425 +0,0 @@ -/* - * 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 Last changed by: $LastChangedBy: jgonzalez $ - * - * Copied from tanktop - */ - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#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/G4BeamTestHadronPhysics.cxx.backup b/src/G4BeamTestHadronPhysics.cxx.backup new file mode 100644 index 0000000..1e37ed5 --- /dev/null +++ b/src/G4BeamTestHadronPhysics.cxx.backup @@ -0,0 +1,412 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#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 deleted file mode 100644 index 3c7e553..0000000 --- a/src/G4BeamTestIonPhysics.cxx +++ /dev/null @@ -1,112 +0,0 @@ -/* - * 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 , Peter Nießen (tanktop) - * Last changed by: $LastChangedBy: kislat $ - */ - -#include - -#include -#include -#include -#include -#include -#include - -#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/G4BeamTestIonPhysics.cxx.backup b/src/G4BeamTestIonPhysics.cxx.backup new file mode 100644 index 0000000..a8cc573 --- /dev/null +++ b/src/G4BeamTestIonPhysics.cxx.backup @@ -0,0 +1,100 @@ +#include + +#include +#include +#include +#include +#include +#include + +#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 index 92264a2..2c5404e 100644 --- a/src/G4BeamTestMuonPhysics.cxx +++ b/src/G4BeamTestMuonPhysics.cxx @@ -1,16 +1,3 @@ -/* - * 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 , Peter Nießen (tanktop) - * Last changed by: $LastChangedBy: kislat $ - */ - - #include #include #include @@ -78,7 +65,7 @@ void G4BeamTestMuonPhysics::ConstructProcess() pManager->AddProcess(&muMinusBremsstrahlung_, -1, 3, 3); pManager->AddProcess(&muMinusPairProduction_, -1, 4, 4); - pManager->AddRestProcess(&muMinusCaptureAtRest_); + pManager->AddRestProcess(&muMinusCapture_); // Tau Plus Physics pManager = G4TauPlus::TauPlus()->GetProcessManager(); diff --git a/src/G4BeamTestPhysicsList.cxx b/src/G4BeamTestPhysicsList.cxx index 5a5befd..0c27b0b 100644 --- a/src/G4BeamTestPhysicsList.cxx +++ b/src/G4BeamTestPhysicsList.cxx @@ -1,15 +1,3 @@ -/** - * 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 #include #include "G4BeamTestPhysicsList.h" diff --git a/src/G4BeamTestRunManager.cxx b/src/G4BeamTestRunManager.cxx index 753deb8..5ed2e05 100644 --- a/src/G4BeamTestRunManager.cxx +++ b/src/G4BeamTestRunManager.cxx @@ -1,20 +1,8 @@ -/** - * 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 +#include "G4BeamTestRunManager.h" #include #include diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx index 4f12c7a..2fc5bed 100644 --- a/src/G4BeamTestTank.cxx +++ b/src/G4BeamTestTank.cxx @@ -1,25 +1,3 @@ -/** - * 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 -#include -#include - -#include -#include -#include -#include -#include - #include #include #include @@ -33,7 +11,10 @@ #include -#include +/* #include */ + +#include "G4BeamTestTank.h" +#include "G4TankIceSD.h" //prevent gcc to make something stupid with pretended unused variables #ifdef __GNUC__ @@ -43,121 +24,26 @@ #endif -G4BeamTestTank::G4BeamTestTank(const TankKey& tankKey, const I3Geometry& geometry): -tankKey_(tankKey), geometry_(geometry) +G4BeamTestTank::G4BeamTestTank() { - 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()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; + tankThickness_ = 0.5*CLHEP::cm; // TODO(shivesh) : check thickness + tankHeight_ = 76.83 * 2.54 * CLHEP::cm; + innerRadius_ = 32 * 2.54 * CLHEP::cm; 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_; - + + // Get fill height + fillHeight_ = 55.1 * 2.54 * CLHEP::cm; + // 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() { - + } @@ -166,217 +52,219 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const // 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 ?????? + // TODO(shivesh): 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(tankKey_); - - // Define plastic frame + energyLimit->SetUserMinEkine(264.1 * CLHEP::keV); // Cherenkov threshold of electrons in water + + // std::string tankName=boost::lexical_cast(tankKey_); + std::string tankName = "BTT"; + + // Define plastic frame TODO(shivesh): plastic or polyethelene? 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(), + + // Define water volume + G4Material* water = G4Material::GetMaterial("Water"); + G4Tubs* solidWater = new G4Tubs(("solid_water_" + 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_, + G4LogicalVolume* logWater = + new G4LogicalVolume(solidWater, water, ("log_water_" + tankName).c_str(), 0, 0, 0); + G4ThreeVector physWaterPosition(0, 0, -0.5*tankHeight_ + tankThickness_ + 0.5*fillHeight_); + G4VPhysicalVolume* physWater ATTRIBUTE_UNUSED = + new G4PVPlacement(0, physWaterPosition, logWater, + ("water_" + tankName).c_str(), tankLog_, false, 0); + + // Define air volume + G4Material* air = G4Material::GetMaterial("Air"); + G4Tubs* solidAir = new G4Tubs(("solid_air_" + tankName).c_str(), + 0.0 * CLHEP::m, innerRadius_, 0.5 * airHeight_, 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); - + G4LogicalVolume* logAir = + new G4LogicalVolume(solidAir, air, ("log_air_" + tankName).c_str(), 0, 0, 0); + G4ThreeVector physAirPosition(0, 0, -0.5 * tankHeight_ + 0.5 * CLHEP::cm + fillHeight_ + + 0.5 * airHeight_); + G4VPhysicalVolume* physAir_UNUSED = + new G4PVPlacement(0, physAirPosition, logAir, + ("air_" + 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 domPosIce; - std::map::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(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, + + // std::map domPosIce; + // std::map::const_iterator om_iter; + // for(om_iter=relDomPositions_.begin(); om_iter!=relDomPositions_.end(); ++om_iter) + // { + /* const OMKey& omKey = om_iter->first; */ + + G4ThreeVector upperDOMpos(0, 0, -0.5 * airHeight_); + G4ThreeVector lowerDOMpos(0, 0, 0.5 * fillHeight_); + + // domPosIce[omKey] = lowerDOMpos; + + // std::string omName=boost::lexical_cast(omKey); + std::string omName="BTOM"; + + 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 *lowerdomsphere = new G4Sphere (("solid_inside_dom_lo_" + omName).c_str(), - 0.0 * CLHEP::m, domInnerRadius, + 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); - - 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 + + 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(), logAir, false, 0); + G4VPhysicalVolume* physLowerGlass ATTRIBUTE_UNUSED = + new G4PVPlacement(0, lowerDOMpos, logLowerGlass, + ("dom_lo_" + omName).c_str(), logWater, 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 TODO(shivesh): make the PMT the SD G4SDManager* sdManager = G4SDManager::GetSDMpointer(); - iceSD_ = new G4TankIceSD(("ice_SD_" + tankName).c_str(), domPosIce); + iceSD_ = new G4TankIceSD(("ice_SD_" + tankName).c_str()); sdManager->AddNewDetector(iceSD_); - logIce->SetSensitiveDetector(iceSD_); - + logWater->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 + // 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.); -} - + logAir->SetUserLimits(energyLimit); + logWater->SetUserLimits(energyLimit); -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; + return tankPhys; } -double G4BeamTestTank::GetY_I3() +double G4BeamTestTank::GetNumCherenkov(/* const OMKey& omKey */) { - return (position_.y() / CLHEP::m) * I3Units::m; + return std::max(iceSD_->GetNumCherenkov(/* omKey */), 0.); } -double G4BeamTestTank::GetZ_I3() +double G4BeamTestTank::GetNumCherenkovWeight(/* const OMKey& omKey */) { - return (position_.z() / CLHEP::m) * I3Units::m; + return std::max(iceSD_->GetNumCherenkovWeight(/* omKey */), 0.); } -double G4BeamTestTank::GetTankHeight_I3() +double G4BeamTestTank::GetEDep_G4(/* const OMKey& omKey */) { - return (tankHeight_ / CLHEP::m) * I3Units::m; + return std::max(iceSD_->GetEDep(/* omKey */), 0.); } -double G4BeamTestTank::GetTankRadius_I3() +double G4BeamTestTank::GetTime_G4(/* const OMKey& omKey */) { - return (outerRadius_ / CLHEP::m) * I3Units::m; + return iceSD_->GetTime(/* omKey */); } -double G4BeamTestTank::GetSnowHeight_I3() -{ - return (snowHeight_ / CLHEP::m) * I3Units::m; -} +/* 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) */ +/* } */ +/* */ +/* */ +/* double G4BeamTestTank::GetTankRadius_I3() */ +/* { */ +/* return (outerRadius_ / CLHEP::m) */ +/* } */ +/* */ +/* */ +/* double G4BeamTestTank::GetFillHeight_I3() */ +/* { */ +/* return (fillHeight_ / CLHEP::m) */ +/* } */ diff --git a/src/G4BeamTestUserSteppingAction.cxx b/src/G4BeamTestUserSteppingAction.cxx index b16098c..e064628 100644 --- a/src/G4BeamTestUserSteppingAction.cxx +++ b/src/G4BeamTestUserSteppingAction.cxx @@ -1,18 +1,4 @@ -/** - * Copyright (C) 2011 - * The IceCube collaboration - * ID: $Id$ - * - * @file G4BeamTestUserSteppingAction.cxx - * @version $Revision$ - * @date $Date$ - * @author Thomas Melzig - * - * $LastChangedBy$ - */ - - -#include +#include "G4BeamTestUserSteppingAction.h" #include "G4Step.hh" #include "G4Track.hh" diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx index d554152..c6317df 100644 --- a/src/G4BeamTestUserTrackingAction.cxx +++ b/src/G4BeamTestUserTrackingAction.cxx @@ -1,18 +1,4 @@ -/** - * Copyright (C) 2011 - * The IceCube collaboration - * ID: $Id$ - * - * @file G4BeamTestUserTrackingAction.cxx - * @version $Revision$ - * @date $Date$ - * @author Thomas Melzig - * - * $LastChangedBy$ - */ - - -#include +#include "G4BeamTestUserTrackingAction.h" #include "G4Track.hh" #include "G4UserLimits.hh" diff --git a/src/G4Interface.cxx b/src/G4Interface.cxx index 36ba3d9..9f44d3e 100644 --- a/src/G4Interface.cxx +++ b/src/G4Interface.cxx @@ -1,23 +1,12 @@ -/** - * 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 - -#include -#include -#include -#include -#include - -#include +/* #include */ + +#include "G4Interface.h" +#include "G4BeamTestDetectorConstruction.h" +#include "G4BeamTestPhysicsList.h" +#include "G4BeamTestUserTrackingAction.h" +#include "G4BeamTestUserSteppingAction.h" + +/* #include */ #ifdef G4VIS_USE #include @@ -26,7 +15,7 @@ #include #include #include -#include +/* #include */ G4Interface* G4Interface::g4Interface_ = NULL; @@ -40,11 +29,8 @@ G4Interface::G4Interface(const std::string& visMacro): // Visualization manager #ifdef G4VIS_USE visManager_ = NULL; - if(!visMacro_.empty()) - { - visManager_ = new G4VisExecutive(); - visManager_->Initialize(); - } + visManager_ = new G4VisExecutive(); + visManager_->Initialize(); #endif } @@ -63,7 +49,8 @@ void G4Interface::InstallTank(G4BeamTestTank* tank) { if(initialized_) { - log_fatal("G4Interface aleady initialized. Cannot install tank!"); + /* log_fatal("G4Interface aleady initialized. Cannot install tank!"); */ + G4cout << "G4Interface aleady initialized. Cannot install tank!" << G4endl; return; } @@ -92,117 +79,37 @@ void G4Interface::InitializeEvent() } -void G4Interface::InjectParticle(const I3Particle& particle) +void G4Interface::InjectParticle( + const std::string& particleName, const G4ThreeVector& particlePosition, + const G4ThreeVector& particleDirection, const G4double particleEnergy + ) { if(!eventInitialized_) { - log_fatal("No event initialized. Cannot inject particle!"); + /* log_fatal("No event initialized. Cannot inject particle!"); */ + G4cout << "No event initialized. Cannot inject particle!" << G4endl; 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); + particleDef = particleTable->FindParticle(particleName); - // Transform I3 coorinates to world system - position -= detector_->GetWorldOrigin(); - - G4ThreeVector direction(particle.GetDir().GetX(), - particle.GetDir().GetY(), - particle.GetDir().GetZ()); + // Transform coordinates to world system + G4ThreeVector position = particlePosition - detector_->GetWorldOrigin(); G4ParticleGun gun(1); gun.SetParticleDefinition(particleDef); - gun.SetParticleEnergy((particle.GetEnergy() / I3Units::GeV) * CLHEP::GeV); + gun.SetParticleEnergy(particleEnergy); gun.SetParticlePosition(position); - gun.SetParticleMomentumDirection(direction); + gun.SetParticleMomentumDirection(particleDirection); - log_trace("Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV", - particle.GetTypeString().c_str(), + G4cout << "Injecting %s: x=%.2f m, y=%.2f m, z=%.2f m, E=%.3f MeV", + particleName.c_str(), position.x() / CLHEP::m, position.y() / CLHEP::m, position.z() / CLHEP::m, - gun.GetParticleEnergy() / CLHEP::MeV); + gun.GetParticleEnergy() / CLHEP::MeV; runManager_.InjectParticle(&gun); } @@ -227,66 +134,73 @@ void G4Interface::Initialize() { if(initialized_) { - log_error("G4Interface has already been initialized. Ignoring this call!"); + /* log_error("G4Interface has already been initialized. Ignoring this call!"); */ + G4cout << "G4Interface has already been initialized. Ignoring this call!" << G4endl; return; } - log_debug("Init geometry ..."); + /* log_debug("Init geometry ..."); */ + G4cout << "Init geometry ..." << G4endl; runManager_.SetUserInitialization(detector_); - log_debug("Init physics list ..."); + /* log_debug("Init physics list ..."); */ + G4cout << "Init physics list ..." << G4endl; runManager_.SetUserInitialization(new G4BeamTestPhysicsList()); - log_debug("Init UserTrackingAction ..."); + /* log_debug("Init UserTrackingAction ..."); */ + G4cout << "Init UserTrackingAction ..." << G4endl; runManager_.SetUserAction(new G4BeamTestUserTrackingAction()); - log_debug("Init UserSteppingAction ..."); + /* log_debug("Init UserSteppingAction ..."); */ + G4cout << "Init UserSteppingAction ..." << G4endl; runManager_.SetUserAction(new G4BeamTestUserSteppingAction()); // Initialize G4 kernel - log_debug("Init run manager ..."); + /* log_debug("Init run manager ..."); */ + G4cout << "Init run manager ..." << G4endl; 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 + /* // 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; */ + /* } */ + + // TODO(shivesh): verbosity +/* 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()); - } + /* 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 index 4a1be0b..2974837 100644 --- a/src/G4TankIceSD.cxx +++ b/src/G4TankIceSD.cxx @@ -1,17 +1,5 @@ -/** - * 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 -#include +#include "G4TankIceSD.h" +#include "G4BeamTestTank.h" #include #include @@ -25,7 +13,7 @@ #include #include "G4Poisson.hh" -#include +/* #include */ #include "G4Version.hh" @@ -38,8 +26,8 @@ #endif -G4TankIceSD::G4TankIceSD(G4String name, const std::map& domPositions) - : G4VSensitiveDetector(name), domPositions_(domPositions) +G4TankIceSD::G4TankIceSD(G4String name/* , const std::map& domPositions */) + : G4VSensitiveDetector(name)/* , domPositions_(domPositions) */ {} @@ -48,14 +36,14 @@ G4TankIceSD::~G4TankIceSD() {} void G4TankIceSD::Initialize(G4HCofThisEvent* HCE) { - std::map::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; - } + // std::map::const_iterator om_iter; + // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) + // { + sumEdep_ = 0; + cogTime_ = 0; + cherenkovCounter_ = 0; + cherenkovCounterWeight_ = 0; + // } } @@ -76,32 +64,32 @@ G4bool G4TankIceSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) G4ThreeVector localPosition = theTouchable->GetHistory()-> GetTopTransform().TransformPoint(worldPosition); - std::map::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; - } + // std::map::const_iterator om_iter; + // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) + // { + G4double radius = sqrt(pow(domPositions_.x() - localPosition.x(), 2) + + pow(domPositions_.y() - localPosition.y(), 2)); + G4double height = (domPositions_.z() - localPosition.z()); + + sumEdep_ += edep; + cogTime_ += edep*time; + cherenkovCounterWeight_ += GetProbability(radius, height) * cherenkovNumber; + cherenkovCounter_ += cherenkovNumber; + // } return true; } void G4TankIceSD::EndOfEvent(G4HCofThisEvent*) { - std::map::const_iterator om_iter; - for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) + // std::map::const_iterator om_iter; + // for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter) + // { + if(sumEdep_>0) { - if(sumEdep_[om_iter->first]>0) - { - cogTime_[om_iter->first] /= sumEdep_[om_iter->first]; - } + cogTime_ /= sumEdep_; } + // } } @@ -122,7 +110,7 @@ G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep) G4MaterialPropertyVector* Rindex = aMaterialPropertiesTable->GetProperty("RINDEX"); if (!Rindex) return 0.0; - const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm); + const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm); // TODO(shivesh): check this const G4double charge = aParticle->GetDefinition()->GetPDGCharge(); const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.; if (beta <= 0.0) return 0.0; @@ -147,7 +135,8 @@ G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep) #endif if (nMin!=nMax) { - log_error("Support for energy dependent refraction index not yet implemented!"); + /* log_error("Support for energy dependent refraction index not yet implemented!"); */ + G4cout << "Support for energy dependent refraction index not yet implemented!"; return 0.0; } @@ -163,6 +152,7 @@ G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep) } +// TODO(shivesh): what is this G4double G4TankIceSD::GetProbability(G4double radius, G4double height) { height = 0.90 - height / CLHEP::m; -- cgit v1.2.3