From 117da1b2b3197f41621242cdb7f92e52989ed81f Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Thu, 5 Jan 2017 10:31:33 +0100
Subject: [PATCH] * Drift Electron transport no correctly working -
Adjustable step length - good speed at 1mm step - take into
account the local field direction
---
NPSimulation/Detectors/Chio/Chio.cc | 47 +++++++++++++------
NPSimulation/Process/G4DEAmplification.cc | 10 +---
NPSimulation/Process/G4DETransport.cc | 53 ++++++++++++++++------
NPSimulation/Process/G4DETransport.hh | 4 +-
NPSimulation/Process/G4DriftElectron.cc | 6 ++-
NPSimulation/Process/G4IonizationWithDE.cc | 10 ++--
Projects/BeLise/PhysicsListOption.txt | 2 +-
7 files changed, 86 insertions(+), 46 deletions(-)
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index 4ed7fa3ee..9fd9f085c 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -68,11 +68,11 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Chio_NS{
// Energy and time Resolution
-// const double EnergyThreshold = 0.1*MeV;
+ // const double EnergyThreshold = 0.1*MeV;
//const double ResoTime = 4.5*ns ;
- // const double ResoEnergy = 1.0*MeV ;
+ // const double ResoEnergy = 1.0*MeV ;
//const double Radius = 50*mm ;
- // const double Width = 100*mm ;
+ // const double Width = 100*mm ;
const double Thickness = 300*mm ;
const string Material = "BC400";
}
@@ -133,12 +133,12 @@ G4LogicalVolume* Chio::BuildDetector(){
G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe");
G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
-
+
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- MPT->AddConstProperty("DE_YIELD",1e-3);
-// MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
- MPT->AddConstProperty("DE_ABSLENGTH",1*km);
+ MPT->AddConstProperty("DE_YIELD",1e-4);
+ // MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
+ MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
MPT->AddConstProperty("DE_DRIFTSPEED",8e-3*mm/ns);
MPT->AddConstProperty("DE_TRANSVERSALSPREAD",6e-5*mm2/ns);
MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-5*mm2/ns);
@@ -162,14 +162,31 @@ G4LogicalVolume* Chio::BuildDetector(){
logicWindows,
"ChioEntranceWindows",m_SquareDetector,false,0);
-
-/* G4Region* DriftRegion = new G4Region("DriftRegion");
- DriftRegion->AddRootLogicalVolume(logicGas);
- G4ProductionCuts* cuts = new G4ProductionCuts;
- cuts->SetProductionCut(1e-9*micrometer); // same cuts for gamma, e- and e+
- DriftRegion->SetProductionCuts(cuts);
-
- new DriftElectron("DriftElectron",DriftRegion);*/
+ G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-1000*volt/cm,0.0));
+ // Create an equation of motion for this field
+ G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
+ G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 );
+
+ // Get the global field manager
+ G4FieldManager* FieldManager= new G4FieldManager();
+ // Set this field to the global field manager
+ FieldManager->SetDetectorField(field );
+ logicGas->SetFieldManager(FieldManager,true);
+
+ G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm,
+ Stepper,
+ Stepper->GetNumberOfVariables() );
+
+ G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
+ FieldManager->SetChordFinder( ChordFinder );
+
+ /* G4Region* DriftRegion = new G4Region("DriftRegion");
+ DriftRegion->AddRootLogicalVolume(logicGas);
+ G4ProductionCuts* cuts = new G4ProductionCuts;
+ cuts->SetProductionCut(1e-9*micrometer); // same cuts for gamma, e- and e+
+ DriftRegion->SetProductionCuts(cuts);
+
+ new DriftElectron("DriftElectron",DriftRegion);*/
m_SquareDetector->SetVisAttributes(m_VisChamber);
logicGas->SetVisAttributes(m_VisGas);
logicWindows->SetVisAttributes(m_VisWindows);
diff --git a/NPSimulation/Process/G4DEAmplification.cc b/NPSimulation/Process/G4DEAmplification.cc
index a171c3187..2db666141 100644
--- a/NPSimulation/Process/G4DEAmplification.cc
+++ b/NPSimulation/Process/G4DEAmplification.cc
@@ -108,8 +108,6 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
-
-
G4ThreeVector x0 = pPreStepPoint->GetPosition();
G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
G4double t0 = pPreStepPoint->GetGlobalTime();
@@ -158,13 +156,8 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4ThreeVector p;
p.setRThetaPhi(1,theta,phi);
- // 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);
+ G4Track* aSecondaryTrack = new G4Track(particle,t0,x0);
aSecondaryTrack->SetTouchableHandle(
aStep.GetPreStepPoint()->GetTouchableHandle());
@@ -174,6 +167,7 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
aParticleChange.AddSecondary(aSecondaryTrack);
}
+ // The original track is left intact
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc
index 0e721759d..1247dabac 100644
--- a/NPSimulation/Process/G4DETransport.cc
+++ b/NPSimulation/Process/G4DETransport.cc
@@ -49,6 +49,8 @@
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4DETransport.hh"
+#include "G4FieldManager.hh"
+#include "G4ElectroMagneticField.hh"
/////////////////////////
// Class Implementation
@@ -90,7 +92,6 @@ G4DETransport::~G4DETransport(){}
////////////
// Methods
-////////////
// PostStepDoIt
// -------------
@@ -99,6 +100,7 @@ G4VParticleChange*
G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
{
aParticleChange.Initialize(aTrack);
+
G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
G4ThreeVector x0 = pPreStepPoint->GetPosition();
@@ -107,6 +109,13 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4double Energy = pPreStepPoint->GetKineticEnergy();
+ // The time scale is imposed by the drift speed
+ G4double step = aStep.GetDeltaTime();
+
+ if(!step){ // allow internal relocation of the track by the kernel taking a 0 length intermediate step
+ return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
+ }
+
// Get the material table
const G4Material* aMaterial = aTrack.GetMaterial();
@@ -125,15 +134,24 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
// Electron follow the field direction
- G4ThreeVector driftDir(0,1,0);
+ // 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();
+
+ G4ElectroMagneticField* field = (G4ElectroMagneticField*)fMng->GetDetectorField();
+ field->GetFieldValue(Point,fieldArr) ;
+
+ // Electron move opposite to the field direction, hance the minus sign
+ G4ThreeVector driftDir(-fieldArr[3],-fieldArr[4],-fieldArr[5]);
+ // Normalised the drift direction
+ driftDir = driftDir.unit();
G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
- // The time scale is imposed by the drift speed
- G4double step = aStep.GetDeltaTime();
- if(!step)
- return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
-
- // Should be equal to delta pos y
- G4double step_length = step*v_drift;
+ // Should be equal to delta length
+ G4double step_length = aStep.GetDeltaPosition().mag();
+ step = step_length/v_drift;
+
G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift);
G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift);
@@ -143,18 +161,22 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4ThreeVector trans(G4RandGauss::shoot(0,d_trans),0,0);
trans.rotateY(twopi*G4UniformRand());
G4ThreeVector d_Pos = (step_length+d_long)*driftDir+trans;
- // Random Position along the step with matching time
+
+ // New particle Position with matching time
G4ThreeVector pos = x0 + d_Pos;
G4double time = t0 + step;
-
// Garanty that the electron does not jump outside the current volume
- m_SafetyHelper->ReLocateWithinVolume(pos);
+ G4double safety = m_SafetyHelper->ComputeSafety(x0,d_Pos.mag());
+
+ if(d_Pos.mag()>safety){
+ pos = x0+safety*d_Pos.unit();
+ }
aParticleChange.ProposeMomentumDirection(d_Pos.unit());
aParticleChange.ProposeEnergy(Energy);
aParticleChange.ProposePosition(pos);
aParticleChange.ProposeGlobalTime(time);
-
+ aParticleChange.ProposeVelocity(v_drift/c_light);
return &aParticleChange;
}
@@ -166,5 +188,8 @@ G4double G4DETransport::GetMeanFreePath(const G4Track& ,
G4double ,
G4ForceCondition* )
{
- return 0.5*mm;
+ // Typical distance after which the electron position should be reevaluated
+ // to take into account the diffusivity of the charge cloud inside the
+ // medium
+ return 1*mm;
}
diff --git a/NPSimulation/Process/G4DETransport.hh b/NPSimulation/Process/G4DETransport.hh
index c85e2f588..d885efe8b 100644
--- a/NPSimulation/Process/G4DETransport.hh
+++ b/NPSimulation/Process/G4DETransport.hh
@@ -105,8 +105,8 @@ public:
G4double GetMeanFreePath(const G4Track& aTrack,
G4double ,
G4ForceCondition* );
- // Returns the transport length for bulk transport of drift
- // electron in media with a specified attenuation length.
+ // Returns the typical travel length for transport of drift
+ // electron in media
G4VParticleChange* PostStepDoIt(const G4Track& aTrack,
const G4Step& aStep);
diff --git a/NPSimulation/Process/G4DriftElectron.cc b/NPSimulation/Process/G4DriftElectron.cc
index 04a09993a..d64f55a76 100644
--- a/NPSimulation/Process/G4DriftElectron.cc
+++ b/NPSimulation/Process/G4DriftElectron.cc
@@ -53,6 +53,9 @@ G4DriftElectron* G4DriftElectron::Definition()
// search in particle table]
G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* anInstance = pTable->FindParticle(name);
+ // Note: The charge of drift electron is null so they are not afected by the
+ // electro magnetic field directly. This is taken care of "manually" by the
+ // G4DETransport process that get the local field direction
if (anInstance ==0)
{
// create particle
@@ -65,7 +68,7 @@ G4DriftElectron* G4DriftElectron::Definition()
// stable lifetime decay table
// shortlived subType anti_encoding
anInstance = new G4ParticleDefinition(
- name, electron_mass_c2, 0.0*MeV, -1.*eplus,
+ name, electron_mass_c2, 0.0*MeV, -0.*eplus,
1, 0, 0,
0, 0, 0,
"lepton", 1, 0, 11,
@@ -73,6 +76,7 @@ G4DriftElectron* G4DriftElectron::Definition()
false, "e"
);
}
+
theInstance = reinterpret_cast<G4DriftElectron*>(anInstance);
return theInstance;
}
diff --git a/NPSimulation/Process/G4IonizationWithDE.cc b/NPSimulation/Process/G4IonizationWithDE.cc
index 2856283b4..b2aa522a4 100644
--- a/NPSimulation/Process/G4IonizationWithDE.cc
+++ b/NPSimulation/Process/G4IonizationWithDE.cc
@@ -163,14 +163,16 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
- G4double pair_energy=0;
- pair_energy=
+ 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;
+ number_electron = G4Poisson(number_electron);
//if no electron leave
if(number_electron<1){
aParticleChange.SetNumberOfSecondaries(0);
@@ -195,9 +197,7 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy);
G4Track* aSecondaryTrack = new G4Track(particle,time,pos);
-
- aSecondaryTrack->SetTouchableHandle(
- aStep.GetPreStepPoint()->GetTouchableHandle());
+ aSecondaryTrack->SetVelocity(v_drift/c_light);
aSecondaryTrack->SetParentID(aTrack.GetTrackID());
aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
diff --git a/Projects/BeLise/PhysicsListOption.txt b/Projects/BeLise/PhysicsListOption.txt
index d85f8d5db..8c522e9f8 100644
--- a/Projects/BeLise/PhysicsListOption.txt
+++ b/Projects/BeLise/PhysicsListOption.txt
@@ -1,5 +1,5 @@
EmPhysicsList Option4
-DefaultCutOff 100
+DefaultCutOff 1000
IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0
EmExtraPhysics 0
--
GitLab