From 7cf84094dd06f8b56b05b38a8df33f0b29117dfd Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Wed, 4 Jul 2018 17:33:35 +0200
Subject: [PATCH] *Updating G4DE physics process for faster simulation
- Number of DE produce by ionisation is now fixed per mm - YIELD no
longer used - DE track carry weight of the DE - DE scorer
read out hte weight - Amplification also use weight and fixed amount
of generated electron This should insure smoother performance with
regards to different setup
---
NPSimulation/Detectors/Chio/Chio.cc | 3 +-
NPSimulation/Process/G4DEAmplification.cc | 65 ++++++++-------
NPSimulation/Process/G4DETransport.cc | 1 -
NPSimulation/Process/G4IonizationWithDE.cc | 88 ++++++++++++--------
NPSimulation/Scorers/DriftElectronScorers.cc | 2 +-
nptool.sh | 12 +--
6 files changed, 94 insertions(+), 77 deletions(-)
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index a02d0cdee..97d811c6e 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -153,7 +153,6 @@ G4LogicalVolume* Chio::BuildDetector(){
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- MPT->AddConstProperty("DE_YIELD",1e-2);
// MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
MPT->AddConstProperty("DE_DRIFTSPEED",11*cm/microsecond);
@@ -181,7 +180,7 @@ G4LogicalVolume* Chio::BuildDetector(){
logicGas,
"ChioGas",m_SquareDetector,false,0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(-2.5*cm,0,0)),
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(2.5*cm,0,0)),
logicGrid,
"ChioGrid",logicGas,false,0);
diff --git a/NPSimulation/Process/G4DEAmplification.cc b/NPSimulation/Process/G4DEAmplification.cc
index 8e8b3e272..22b273de0 100644
--- a/NPSimulation/Process/G4DEAmplification.cc
+++ b/NPSimulation/Process/G4DEAmplification.cc
@@ -98,19 +98,11 @@ G4DEAmplification::~G4DEAmplification(){}
G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
{
-// if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" )
-// return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
-
// Get the primary track
aParticleChange.Initialize(aTrack);
const G4Material* aMaterial = aTrack.GetMaterial();
- G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
-
- G4ThreeVector x0 = pPreStepPoint->GetPosition();
- G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
- G4double t0 = pPreStepPoint->GetGlobalTime();
// Get the material table
G4MaterialPropertiesTable* aMaterialPropertiesTable =
@@ -118,41 +110,50 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
if (!aMaterialPropertiesTable){ // Does the table exist
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
-
- if(!aMaterialPropertiesTable->ConstPropertyExists("DE_AMPLIFICATION") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD"))
+
+ if(!aMaterialPropertiesTable->ConstPropertyExists("DE_AMPLIFICATION"))
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
-
+
+
+ G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
+
+ G4ThreeVector x0 = pPreStepPoint->GetPosition();
+ G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
+ G4double t0 = pPreStepPoint->GetGlobalTime();
+
G4double pair_energy = pPreStepPoint->GetKineticEnergy();
G4double amplification = aMaterialPropertiesTable->GetConstProperty("DE_AMPLIFICATION");
- G4double Yield = aMaterialPropertiesTable->GetConstProperty("DE_YIELD");
- G4int number_electron = amplification*Yield;
- //if no electron leave
- if(number_electron<1){
+
+ double OriginalW = aTrack.GetWeight();
+ // Number of Physical electron to be created
+ G4int number_electron = amplification*OriginalW;
+ //Number of tracked electron effectively created instead
+ G4int tracked_electron = 10;
+
+ //if no electron leave
+ if(tracked_electron<1){
aParticleChange.SetNumberOfSecondaries(0);
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
- aParticleChange.SetNumberOfSecondaries(number_electron);
+ static G4ParticleDefinition* DEDef = G4DriftElectron::DriftElectron();
+ G4TouchableHandle handle = pPreStepPoint->GetTouchableHandle();
+ G4int ParentID = aTrack.GetTrackID();
+
+ aParticleChange.SetNumberOfSecondaries(tracked_electron);
+ aParticleChange.SetSecondaryWeightByProcess(true);
// Create the secondary tracks
- for(G4int i = 0 ; i < number_electron ; i++){
- // Random direction at creation
- G4double cost = 1-2*G4UniformRand();
- G4double theta = acos(cost);
- G4double phi = twopi*G4UniformRand();
- G4ThreeVector p;
- p.setRThetaPhi(1,theta,phi);
-
- G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy);
+ for(G4int i = 0 ; i < tracked_electron ; i++){
+ G4DynamicParticle* particle = new G4DynamicParticle(DEDef,p0, pair_energy);
+ // The secondary track follow the original one
G4Track* aSecondaryTrack = new G4Track(particle,t0,x0);
- aSecondaryTrack->SetTouchableHandle(
- aStep.GetPreStepPoint()->GetTouchableHandle());
aSecondaryTrack->SetCreatorProcess(this);
- aSecondaryTrack->SetParentID(aTrack.GetTrackID());
- aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
+ aSecondaryTrack->SetParentID(ParentID);
+ aSecondaryTrack->SetTouchableHandle(handle);
+ aSecondaryTrack->SetWeight(number_electron/tracked_electron);
aParticleChange.AddSecondary(aSecondaryTrack);
}
// to kill the primary track to avoid it being re-amplified
@@ -183,10 +184,10 @@ G4double G4DEAmplification::GetMeanLifeTime(const G4Track& aTrack,
G4ForceCondition* condition)
{
if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" ){
- *condition = NotForced;
+ *condition = NotForced;
return DBL_MAX;
}
-
+
*condition = StronglyForced;
return DBL_MAX;
}
diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc
index aa7a75fbc..f37b3995d 100644
--- a/NPSimulation/Process/G4DETransport.cc
+++ b/NPSimulation/Process/G4DETransport.cc
@@ -128,7 +128,6 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") )
diff --git a/NPSimulation/Process/G4IonizationWithDE.cc b/NPSimulation/Process/G4IonizationWithDE.cc
index b4a2abe6c..21b70fda9 100644
--- a/NPSimulation/Process/G4IonizationWithDE.cc
+++ b/NPSimulation/Process/G4IonizationWithDE.cc
@@ -117,7 +117,6 @@ void G4IonizationWithDE::BuildPhysicsTable(const G4ParticleDefinition&)
//
G4VParticleChange*
G4IonizationWithDE::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
-
// This routine simply calls the equivalent PostStepDoIt since all the
// necessary information resides in aStep.GetTotalEnergyDeposit()
@@ -140,15 +139,6 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
aParticleChange.Initialize(aTrack);
const G4Material* aMaterial = aTrack.GetMaterial();
-
- G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
-
- G4ThreeVector x0 = pPreStepPoint->GetPosition();
- G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
- G4double t0 = pPreStepPoint->GetGlobalTime();
-
- G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
-
// Get the material table
G4MaterialPropertiesTable* aMaterialPropertiesTable =
aMaterial->GetMaterialPropertiesTable();
@@ -157,59 +147,87 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD") ||
!aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") )
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
+ G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
+ const G4TouchableHandle& handle = pPreStepPoint->GetTouchableHandle();
+
+ G4ThreeVector x0 = pPreStepPoint->GetPosition();
+ G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
+ G4double t0 = pPreStepPoint->GetGlobalTime();
+
+ G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
+
+
G4double v_drift=
aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
G4double pair_energy=
aMaterialPropertiesTable->GetConstProperty("DE_PAIRENERGY");
- G4double IonizationWithDEYield = 0;
- IonizationWithDEYield=
- aMaterialPropertiesTable->GetConstProperty("DE_YIELD");
- G4int number_electron = IonizationWithDEYield*TotalEnergyDeposit/pair_energy;
+ // Physical number of electron produced
+ G4int number_electron = TotalEnergyDeposit/pair_energy;
number_electron = G4Poisson(number_electron);
- //if no electron leave
- if(number_electron<1){
+
+ // Tracked electron produced
+ // 100 per mm per step to have a good statistical accuracy
+
+ G4int tracked_electron = (aStep.GetStepLength()/mm)*5;
+
+ //if no electron leave
+ if(tracked_electron<1){
aParticleChange.SetNumberOfSecondaries(0);
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
- aParticleChange.SetNumberOfSecondaries(number_electron);
+ aParticleChange.SetNumberOfSecondaries(tracked_electron);
+ aParticleChange.SetSecondaryWeightByProcess(true);
+
+
- // Create the secondary tracks
- for(G4int i = 0 ; i < number_electron ; i++){
// Electron follow the field direction
// The field direction is taken from the field manager
- G4double* fieldArr = new G4double[6];
- G4double Point[4]={x0.x(),x0.y(),x0.z(),t0};
- G4FieldManager* fMng = pPreStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->
- GetFieldManager();
+ //Everything common to all created DE:
+ static G4double fieldArr[6];
+ static G4double Point[4];
+ Point[0] = x0.x();
+ Point[1] = x0.y();
+ Point[2] = x0.z();
+ Point[3] = t0;
+
+ G4FieldManager* fMng = handle->GetVolume()->GetLogicalVolume()->GetFieldManager();
G4ElectroMagneticField* field = (G4ElectroMagneticField*)fMng->GetDetectorField();
- field->GetFieldValue(Point,fieldArr) ;
+ field->GetFieldValue(Point,&fieldArr[0]) ;
// Electron move opposite to the field direction, hance the minus sign
G4ThreeVector p(-fieldArr[3],-fieldArr[4],-fieldArr[5]);
// Normalised the drift direction
p = p.unit();
- // Random Position along the step with matching time
+ G4ThreeVector delta = aStep.GetDeltaPosition();
+ G4double deltaT = aStep.GetDeltaTime();
+ static G4ParticleDefinition* DEDefinition = G4DriftElectron::DriftElectron();
+ G4double velocity = v_drift/c_light;
+ G4int parentID = aTrack.GetTrackID();
+
+ // Create the secondary tracks
+ for(G4int i = 0 ; i < tracked_electron ; i++){ // always create 100 electron, but with weight
+ // Random Position along the step with matching time
G4double rand = G4UniformRand();
- G4ThreeVector pos = x0 + rand * aStep.GetDeltaPosition();
- G4double time = t0+ rand* aStep.GetDeltaTime();
-
- G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy);
- G4Track* aSecondaryTrack = new G4Track(particle,time,pos);
- aSecondaryTrack->SetVelocity(v_drift/c_light);
-
- aSecondaryTrack->SetParentID(aTrack.GetTrackID());
- aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
+ G4ThreeVector pos = x0 + rand * delta;
+ G4double time = t0 + rand* deltaT;
+
+ G4DynamicParticle* particle = new G4DynamicParticle(DEDefinition,p, pair_energy);
+ G4Track* aSecondaryTrack = new G4Track(particle,time,pos);
+ aSecondaryTrack->SetVelocity(velocity);
+ aSecondaryTrack->SetParentID(parentID);
+ aSecondaryTrack->SetTouchableHandle(handle);
+ // Set the weight, ie, how many electron the track represents
+ aSecondaryTrack->SetWeight(number_electron/tracked_electron);
aParticleChange.AddSecondary(aSecondaryTrack);
}
diff --git a/NPSimulation/Scorers/DriftElectronScorers.cc b/NPSimulation/Scorers/DriftElectronScorers.cc
index 2f66eaeec..75489de5c 100644
--- a/NPSimulation/Scorers/DriftElectronScorers.cc
+++ b/NPSimulation/Scorers/DriftElectronScorers.cc
@@ -57,7 +57,7 @@ G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){
G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName();
if(PID=="driftelectron"){
- Infos[0] = 1;
+ Infos[0] = aStep->GetTrack()->GetWeight();
}
// Check if the particle has interact before, if yes, add up the number of electron.
diff --git a/nptool.sh b/nptool.sh
index 3b15a1257..bb3fdb8ce 100755
--- a/nptool.sh
+++ b/nptool.sh
@@ -37,15 +37,15 @@ NPARCH=$(uname)
# mac os x case
if [ "${NPARCH}" = "Darwin" ] ;
then
- ${CMD} DYLD_LIBRARY_PATH${SEP}$DYLD_LIBRARY_PATH:$NPTOOL/NPLib/lib
- ${CMD} DYLD_LIBRARY_PATH${SEP}$DYLD_LIBRARY_PATH:$NPTOOL/NPSimulation/lib
+ ${CMD} DYLD_LIBRARY_PATH${SEP}$NPTOOL/NPLib/lib:$DYLD_LIBRARY_PATH
+ ${CMD} DYLD_LIBRARY_PATH${SEP}$NPTOOL/NPSimulation/lib:$DYLD_LIBRARY_PATH
else
- ${CMD} LD_LIBRARY_PATH${SEP}$LD_LIBRARY_PATH:$NPTOOL/NPLib/lib
- ${CMD} LD_LIBRARY_PATH${SEP}$LD_LIBRARY_PATH:$NPTOOL/NPSimulation/lib
+ ${CMD} LD_LIBRARY_PATH${SEP}$NPTOOL/NPLib/lib:$LD_LIBRARY_PATH
+ ${CMD} LD_LIBRARY_PATH${SEP}$NPTOOL/NPSimulation/lib:$LD_LIBRARY_PATH
fi
-${CMD} PATH=$PATH:$NPTOOL/NPLib/bin
-${CMD} PATH=$PATH:$NPTOOL/NPSimulation/bin
+${CMD} PATH=$NPTOOL/NPLib/bin:$PATH
+${CMD} PATH=$NPTOOL/NPSimulation/bin:$PATH
alias npt='cd $NPTOOL'
alias npl='cd $NPTOOL/NPLib'
--
GitLab