diff --git a/Inputs/DetectorConfiguration/csi.detector b/Inputs/DetectorConfiguration/csi.detector
index 13e38411fe6a94d6f74da4f20c14ad3476284a70..0d24f8f85966ba43c3f4f9f7b067674f1fa96cd8 100644
--- a/Inputs/DetectorConfiguration/csi.detector
+++ b/Inputs/DetectorConfiguration/csi.detector
@@ -8,7 +8,7 @@ Target
THICKNESS= 0.001
ANGLE= 0
RADIUS= 2
- MATERIAL= CH2
+ MATERIAL= CH2
X= 0
Y= 0
Z= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 1fed2da9732e26e71a224917753c257bb81ffbde..b482984f1198ad3b5705566b25d7262f4a9fbcf5 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -4,14 +4,14 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 5.2
- EnergyHigh= 5.2
+ EnergyLow= 10
+ EnergyHigh= 10
HalfOpenAngleMin= 0
HalfOpenAngleMax= 45
x0= 0
y0= 0
z0= 0
- Particle= alpha
+ Particle= alpha
ExcitationEnergy= 0
% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.cc b/NPSimulation/Core/EventGeneratorIsotropic.cc
index 7c018a95d3d386ac39b8c6deb5aadf6173f6bdbd..7437ab295e5063cda414ed99fef14d1bc7728305 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.cc
+++ b/NPSimulation/Core/EventGeneratorIsotropic.cc
@@ -32,7 +32,6 @@
// NPS headers
#include "EventGeneratorIsotropic.hh"
-
// NPl headers
#include "RootOutput.h"
#include "NPNucleus.h"
@@ -202,7 +201,7 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
void EventGeneratorIsotropic::GenerateEvent(G4Event*){
if(m_particle==NULL){
- if(m_particleName=="gamma" || m_particleName=="neutron"){
+ if(m_particleName=="gamma" || m_particleName=="neutron" || m_particleName=="opticalphoton"){
m_particle = G4ParticleTable::GetParticleTable()->FindParticle(m_particleName.c_str());
}
else{
@@ -212,8 +211,7 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
}
}
-
- // Clear TInitialConditions
+
G4double cos_theta_min = cos(m_HalfOpenAngleMin);
G4double cos_theta_max = cos(m_HalfOpenAngleMax);
G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index 666f6df87d89fbeac085e4ea670cc2c9b61cf70b..4ac63cfad317ef7dcab04823257ff402b4b70f5c 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -256,14 +256,28 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
material->AddElement(GetElementFromLibrary("Cs"),1);
material->AddElement(GetElementFromLibrary("I"),1);
// Adding Scintillation property:
- int NumberOfPoints = 100;
+ int NumberOfPoints = 10;
double wlmin = 0.25*um;
double wlmax = 67*um;
double step = (wlmax-wlmin)/NumberOfPoints;
- double* energy = new double[NumberOfPoints];
+ double* energy_r = new double[NumberOfPoints];
double* rindex = new double[NumberOfPoints];
double* absorption= new double[NumberOfPoints];
+
+ double* energy_e = new double[2];
+ double* fast = new double[2];
+ double* slow = new double[2];
+ double* scint = new double[2];
+
+ energy_e[0] = h_Planck*c_light / (550*nm);
+ energy_e[1] = h_Planck*c_light / (550*nm);
+
+ fast[0] = 1 ; fast[1]=1;
+ slow[0] = 1 ; slow[1]=1;
+ scint[0] = 1; scint[1] = 1;
+
double wl;
+
for(int i = 0 ; i < NumberOfPoints ;i++){
wl= wlmin+i*step;
// Formula from www.refractiveindex.info
@@ -276,18 +290,24 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
+0.01918/(1-pow(0.218/wl,2))
+3.38229/(1-pow(161.29/wl,2))) ;
- energy[i] = h_Planck*c_light / wl;
+ energy_r[i] = h_Planck*c_light / wl;
// To be defined properly
- // absorption[i] = 344.8*cm;
+ absorption[i] = 344.8*cm;
}
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
// From St Gobain
MPT -> AddConstProperty("SCINTILLATIONYIELD",54/keV);
- MPT -> AddProperty("RINDEX",energy,rindex,NumberOfPoints)->SetSpline(true);
- // MPT -> AddProperty("ABSLENGTH",energy,absorption,NumberOfPoints)->SetSpline(true);
- delete energy; delete rindex ; delete absorption;
+ MPT -> AddProperty("SCINTILLATION",energy_e,scint,NumberOfPoints) ;
+ MPT -> AddProperty("RINDEX",energy_r,rindex,NumberOfPoints) ;
+ MPT -> AddProperty("ABSLENGTH",energy_r,absorption,NumberOfPoints);
+ MPT->AddProperty("FASTCOMPONENT", energy_e, fast, NumberOfPoints);
+ MPT->AddProperty("SLOWCOMPONENT", energy_e, slow, NumberOfPoints);
+ MPT->AddConstProperty("RESOLUTIONSCALE",1.0);
+ MPT->AddConstProperty("FASTTIMECONSTANT",1000*ns);
+ MPT->AddConstProperty("SLOWTIMECONSTANT",1000*ns);
+ MPT->AddConstProperty("YIELDRATIO",1.0);
material -> SetMaterialPropertiesTable(MPT);
m_Material[Name]=material;
return material;
diff --git a/NPSimulation/Core/PhysicsList.cc b/NPSimulation/Core/PhysicsList.cc
index 37d45597d2adf399b79717ac3e6be1bc4b78bbc7..7678b309218cff73bacf93613c381093fbc82ec0 100644
--- a/NPSimulation/Core/PhysicsList.cc
+++ b/NPSimulation/Core/PhysicsList.cc
@@ -23,33 +23,20 @@
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
-#include "G4Region.hh"
-#include "G4RegionStore.hh"
#include "PhysicsList.hh"
#include "G4PhysListFactory.hh"
#include "G4VPhysicsConstructor.hh"
// Physics List
-#include "G4NeutronTrackingCut.hh"
#include "G4LossTableManager.hh"
#include "G4UnitsTable.hh"
#include "G4ProcessManager.hh"
-#include "G4IonFluctuations.hh"
-#include "G4IonParametrisedLossModel.hh"
#include "G4EmProcessOptions.hh"
-#include "G4ParallelWorldPhysics.hh"
-#include "G4EmLivermorePhysics.hh"
-#include "G4AutoDelete.hh"
-
/////////////////////////////////////////////////////////////////////////////
PhysicsList::PhysicsList() : G4VModularPhysicsList(){
ReadConfiguration("PhysicsListOption.txt");
G4LossTableManager::Instance();
- defaultCutValue = 1*pc;
- cutForGamma = defaultCutValue;
- cutForElectron = defaultCutValue;
- cutForPositron = defaultCutValue;
-
+ defaultCutValue = 1*mm;
SetVerboseLevel(0);
// ****** Definition of defaults for the physics processes *****
@@ -66,6 +53,7 @@ PhysicsList::PhysicsList() : G4VModularPhysicsList(){
// Elecromagnetic physics
// Using the more accurate option4
emPhysicsList = new G4EmStandardPhysics_option4();
+ opticalPhysicsList = NULL;
// Hadronic physics
if(m_IonBinaryCascadePhysics){
@@ -84,16 +72,20 @@ PhysicsList::PhysicsList() : G4VModularPhysicsList(){
m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP();
// Optical Photon for scintillator simulation
- if(m_OpticalPhysics)
- m_PhysList["OpticalPhysics"] = new G4OpticalPhysics();
-
+ if(m_OpticalPhysics){
+ opticalPhysicsList = new G4OpticalPhysics(0);
+ opticalPhysicsList->SetMaxNumPhotonsPerStep(100);
+ opticalPhysicsList->SetScintillationYieldFactor(0.01);
+ opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true);
+ opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true);
+ //RegisterPhysics(opticalPhysicsList);
+ }
// Decay physics
// Add Radioactive decay
// Gamma decay of known states
if(m_Decay){
- std::cout << "Decay is activated: m_Decay=" << m_Decay << std::endl;
decay_List = new G4DecayPhysics();
radioactiveDecay_List = new G4RadioactiveDecayPhysics() ;
m_PhysList["decay_list"]= decay_List;
@@ -148,7 +140,7 @@ void PhysicsList::ReadConfiguration(std::string filename){
// Most special process need decay to be activated
if( (m_IonBinaryCascadePhysics || m_EmExtraPhysics || m_HadronElasticPhysics
- || m_StoppingPhysics || m_OpticalPhysics || m_HadronPhysicsQGSP_BIC_HP) && !m_Decay){
+ || m_StoppingPhysics || m_HadronPhysicsQGSP_BIC_HP) && !m_Decay){
m_Decay = true;
std::cout << "Information: Selected process require Decay to be activated." << std::endl;
}
@@ -160,11 +152,18 @@ PhysicsList::~PhysicsList(){
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructParticle(){
+ if(m_OpticalPhysics){
+ //G4OpticalPhoton::OpticalPhotonDefinition();
+ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle();
+
+ }
+
if(decay_List){
decay_List -> ConstructParticle();
radioactiveDecay_List->ConstructParticle();
}
- else{
+
+ else{
// If decay is not activated we have to declare the particle ourself
G4He3::He3Definition();
G4Deuteron::DeuteronDefinition();
@@ -193,6 +192,7 @@ void PhysicsList::ConstructParticle(){
G4AntiKaonZero ::AntiKaonZeroDefinition();
G4KaonZeroLong ::KaonZeroLongDefinition();
G4KaonZeroShort::KaonZeroShortDefinition();
+ // Ion
G4IonConstructor ionConstructor ;
ionConstructor.ConstructParticle() ;
}
@@ -205,15 +205,16 @@ void PhysicsList::ConstructProcess(){
// Electromagnetic physics
emPhysicsList -> ConstructProcess();
- em_config.AddModels();
-
+
+ if(opticalPhysicsList){
+ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess();
+ }
// Hadronic physics
std::map<std::string,G4VPhysicsConstructor*>::iterator it;
for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){
it->second -> ConstructProcess();
}
BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
- SetCuts();
return;
}
/////////////////////////////////////////////////////////////////////////////
@@ -230,7 +231,7 @@ void PhysicsList::SetCuts(){
// Special Cut for optical photon to be emmitted
SetCutsWithDefault();
- SetCutValue(1*um,"opticalphoton");
+// SetCutValue(10*um,"opticalphoton");
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/NPSimulation/Core/PhysicsList.hh b/NPSimulation/Core/PhysicsList.hh
index 4dff0558c54e185e7ef8f08a31b58ede4c4e2162..bee34a6cfd4e8c8c15751ed6c0ac2275b01fc960 100644
--- a/NPSimulation/Core/PhysicsList.hh
+++ b/NPSimulation/Core/PhysicsList.hh
@@ -38,7 +38,6 @@
#include "G4EmStandardPhysics_option4.hh"
#include "G4EmExtraPhysics.hh"
#include "G4StoppingPhysics.hh"
-#include "G4DecayPhysics.hh"
#include "G4OpticalPhysics.hh"
#include "G4HadronElasticPhysics.hh"
#include "G4HadronElasticPhysicsHP.hh"
@@ -57,10 +56,6 @@ class PhysicsList: public G4VModularPhysicsList{
void ReadConfiguration(std::string filename);
void ConstructParticle();
void SetCuts();
- void SetCutForGamma(G4double);
- void SetCutForElectron(G4double);
- void SetCutForPositron(G4double);
- void SetDetectorCut(G4double cut);
void ConstructProcess();
void AddStepMax();
void AddPackage(const G4String& name);
@@ -72,9 +67,7 @@ class PhysicsList: public G4VModularPhysicsList{
G4EmConfigurator em_config;
private: // Cuts
- G4double cutForGamma;
- G4double cutForElectron;
- G4double cutForPositron;
+ G4OpticalPhysics* opticalPhysicsList;
G4VPhysicsConstructor* emPhysicsList;
G4VPhysicsConstructor* decay_List;
G4VPhysicsConstructor* radioactiveDecay_List;
diff --git a/NPSimulation/CsI/CsI.cc b/NPSimulation/CsI/CsI.cc
index cbd8dee268fb977e7beb7846a8a6851db79415c8..45537564cc18580d1087953c8cf8dc76fd13e73f 100644
--- a/NPSimulation/CsI/CsI.cc
+++ b/NPSimulation/CsI/CsI.cc
@@ -433,7 +433,7 @@ void CsI::VolumeMaker(G4ThreeVector Det_pos, int DetNumber, G4LogicalVolume* wor
logicPD->SetVisAttributes(PDVisAtt);
new G4PVPlacement(0 ,
- Det_pos+(m_CsIThickness[i]/2+PhotoDiodeThickness/2)*Det_pos.unit() ,
+ Det_pos+(m_CsIThickness[i]*0.5+PhotoDiodeThickness*0.5)*Det_pos.unit() ,
logicPD ,
NamePD ,
world ,
diff --git a/NPSimulation/PhysicsListOption.txt b/NPSimulation/PhysicsListOption.txt
index 768664c8715dce1bac4d080ea8d8e886d8452db8..847ee8d2e01ee6fe3adf4434843e86a2b58c7e06 100644
--- a/NPSimulation/PhysicsListOption.txt
+++ b/NPSimulation/PhysicsListOption.txt
@@ -1,7 +1,7 @@
IonBinaryCascadePhysics 0
-EmExtraPhysics 0
+EmExtraPhysics 0
HadronElasticPhysics 0
StoppingPhysics 0
OpticalPhysics 1
HadronPhysicsQGSP_BIC_HP 0
-Decay 1
+Decay 0
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 992787b018f6997f2923fdc32941616e0c579455..98034832234df5e3c65dbcb6d651e627ec47396b 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -72,19 +72,11 @@ int main(int argc, char** argv){
runManager->SetUserInitialization(detector);
PhysicsList* physicsList = new PhysicsList();
- runManager->SetUserInitialization(physicsList);
-
- // Test for Built in physics list
- // G4PhysListFactory *physListFactory = new G4PhysListFactory();
- //G4VUserPhysicsList *physicsList =
- //physListFactory->GetReferencePhysList("QGSP_BERT");
-
runManager->SetUserInitialization(physicsList);
PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
// Initialize Geant4 kernel
runManager->Initialize();
- //physicsList->MyOwnConstruction();
///////////////////////////////////////////////////////////////
/////////// Define UI terminal for interactive mode ///////////