From 773ec7f2d214e59b0ff6f8727d49ae8eb41bb9a2 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Wed, 4 Dec 2019 17:19:11 +0100
Subject: [PATCH] * Trying to fix DE transport
---
NPSimulation/Process/G4DETransport.cc | 63 ++++++++++-----------------
Projects/Minos/pp.reaction | 6 +--
2 files changed, 25 insertions(+), 44 deletions(-)
diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc
index afd4959df..12b51292a 100644
--- a/NPSimulation/Process/G4DETransport.cc
+++ b/NPSimulation/Process/G4DETransport.cc
@@ -147,18 +147,16 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep){
driftDir = driftDir.unit();
G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
- G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift);
- G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift);
-
+ // Speed sigma in transversal and longitudinal direction based on the
+ // step length, drift speed and long/trans spread parameters
+ G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift)/(aStep.GetDeltaTime());
+ G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift)/(aStep.GetDeltaTime());
+
+ // Building the modified drift vector
+ // i.e. the speed vector of the particle at the end of the step
G4double d_trans = G4RandGauss::shoot(0,sigmaTrans);
- G4double d_long=0;
- G4double d_drift=-1;
+ G4double d_long = G4RandGauss::shoot(v_drift,sigmaLong);
- while(d_drift<0){
- d_long = G4RandGauss::shoot(0,sigmaLong);
- d_drift = step_length+d_long;
- }
-
// The transverse component is perpendicular to the drift direction,
// to build an arbritary trans vector, we have to make a cross product with
// an arbritary, non parallel vector, and driftDir
@@ -174,37 +172,15 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep){
// Rotate randomly around driftDir for the Phi component
trans.rotate(driftDir,twopi*G4UniformRand());
- // new position is drift length*driftDir + the transversal movement
- G4ThreeVector d_Pos = (d_drift)*driftDir+trans;
-
- // Garanty that the electron does not jump outside the current volume
- G4double safety = m_SafetyHelper->ComputeSafety(x0,d_Pos.mag());
-
- // If the distance travelled if above safety, the step is not taken
- /* bool safety_mod = false; */
- /* if(d_Pos.mag()>safety){ */
- /* // multiply by 0.9 to take a 10% safety margin */
- /* safety_mod = true; */
- /* d_drift = d_drift*0.9*safety/d_Pos.mag(); */
- /* d_trans = d_trans*0.9*safety/d_Pos.mag(); */
- /* trans=trans.unit()*d_trans; */
- /* d_Pos = (d_drift)*driftDir+trans; */
- /* } */
-
- // Should be equal to delta length
- G4double step = (d_drift)/v_drift;
-
- // New particle Position with matching time
- G4ThreeVector pos = x0 + d_Pos;
- G4double time = t0 + step;
-
- G4double Energy= (0.5*((electron_mass_c2/joule)/(c_squared/((m*m)/(s*s))))*(v_drift/(m/s))*(v_drift/(m/s)))*joule;
+ // new position is drift vector
+ G4ThreeVector drift_vector=d_long*driftDir.unit()+trans.unit()*d_trans;
+ double speed = drift_vector.mag();
+// std::cout << " " << step_length << " " << speed << " " << d_trans << " " << d_long-v_drift << std::endl;
- aParticleChange.ProposeMomentumDirection(d_Pos.unit());
+ G4double Energy= (0.5*((electron_mass_c2/joule)/(c_squared/((m*m)/(s*s))))*(speed/(m/s))*(speed/(m/s)))*joule;
+ aParticleChange.ProposeMomentumDirection(drift_vector.unit());
aParticleChange.ProposeEnergy(Energy);
- aParticleChange.ProposePosition(pos);
-
- // m_SafetyHelper->ReLocateWithinVolume(pos);
+
return &aParticleChange;
}
@@ -219,7 +195,12 @@ G4double G4DETransport::GetMeanFreePath(const G4Track& track,
// Typical distance after which the electron position should be reevaluated
// to take into account the diffusivity of the charge cloud inside the
// medium
-
- return m_SafetyHelper->ComputeSafety(track.GetPosition());
+ double d = m_SafetyHelper->ComputeSafety(track.GetPosition());
+ if(d>1*mm)
+ return 1*mm;
+ else if (d>1*micrometer)
+ return d;
+ else
+ return 1*pc;
/* return 1.5*mm; */
}
diff --git a/Projects/Minos/pp.reaction b/Projects/Minos/pp.reaction
index 4fe2c02f9..0e5e9244c 100644
--- a/Projects/Minos/pp.reaction
+++ b/Projects/Minos/pp.reaction
@@ -3,7 +3,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
Particle= 1H
- Energy= 550 MeV
+ Energy= 180 MeV
SigmaEnergy= 0 MeV
SigmaThetaX= 0 deg
SigmaPhiY= 0 deg
@@ -24,8 +24,8 @@ TwoBodyReaction
Light= 1H
Heavy= 1H
ExcitationEnergy3= 0.0 MeV
- ExcitationEnergy4= 0 MeV
- CrossSectionPath= 11Li(d,3He)10He.txt CS10He
+ ExcitationEnergy4= 0.0 MeV
+ CrossSectionPath= distrib.txt CS10He
ShootLight= 1
ShootHeavy= 1
--
GitLab