42 #include <CLHEP/Random/Stat.h>
64 #define State(theXInfo) (GetState<G4ITBrownianState>()->theXInfo)
74 #define RED "\033[0;31m"
75 #define LIGHT_RED "\33[1;31m"
76 #define GREEN "\033[32;40m"
77 #define GREEN_ON_BLUE "\033[1;32;44m"
78 #define RESET_COLOR "\033[0m"
83 #define GREEN_ON_BLUE ""
84 #define RESET_COLOR ""
97 return CLHEP::HepStat::inverseErf(x);
102 return CLHEP::HepStat::inverseErf(1. - x);
112 return 1-CLHEP::HepStat::erf(1. - x);
116 #define State(theXInfo) (GetState<G4ITTransportationState>()->theXInfo)
164 if (
this == &rhs)
return *
this;
205 const double timeStep,
221 if (
GetIT(track)->GetTrackingInfo()->IsLeadingStep())
225 bool makeException =
true;
233 exceptionDescription <<
"ComputeStep is called while the track has"
234 "the minimum interaction time";
235 exceptionDescription <<
" so it should not recompute a timeStep ";
236 G4Exception(
"G4DNABrownianTransportation::ComputeStep",
238 exceptionDescription);
242 State(fGeometryLimitedStep) =
false;
254 static double sqrt_2 = sqrt(2.);
255 double sqrt_Dt = sqrt(diffCoeff*timeStep);
256 double sqrt_2Dt = sqrt_2*sqrt_Dt;
258 if(
State(fTimeStepReachedLimit)==
true)
261 State(fGeometryLimitedStep) =
true;
263 spaceStep =
State(fEndPointDistance);
272 spaceStep = sqrt(x*x + y*y + z*z);
274 if(spaceStep >=
State(fEndPointDistance))
278 State(fGeometryLimitedStep) =
true;
293 <<
"G4ITBrownianTransportation::ComputeStep() : "
294 <<
"Step was limited to boundary"
300 if(
State(fRandomNumber)>=0)
320 double invErfc =
InvErfc(value);
321 spaceStep = invErfc*2*sqrt_Dt;
323 if(
State(fTimeStepReachedLimit)==
false)
326 State(fGeometryLimitedStep) =
false;
347 double min_randomNumber =
Erfc(
State(fEndPointDistance)/2*sqrt_Dt);
348 double value = min_randomNumber+(1-min_randomNumber)*
G4UniformRand();
349 double invErfc =
InvErfc(value);
350 spaceStep = invErfc*2*sqrt_Dt;
351 if(spaceStep >=
State(fEndPointDistance))
354 State(fGeometryLimitedStep) =
true;
357 else if(
State(fTimeStepReachedLimit)==
false)
360 State(fGeometryLimitedStep) =
false;
367 State(fGeometryLimitedStep) =
true;
369 spaceStep =
State(fEndPointDistance);
386 State(fTransportEndPosition)= spaceStep*
394 State(fGeometryLimitedStep) =
false;
405 State(fGeometryLimitedStep) =
false;
410 State(fEndGlobalTimeComputed) =
true;
417 <<
"G4ITBrownianTransportation::ComputeStep() : "
418 <<
" trackID : " << track.
GetTrackID() <<
" : Molecule name: "
423 <<
"Final position:" <<
G4BestUnit(
State(fTransportEndPosition),
"Length")
427 <<
"Final magnitude:" <<
G4BestUnit(
State(fTransportEndPosition).mag(),
"Length")
429 <<
"Diffusion length : "
431 <<
" within time step : " <<
G4BestUnit(timeStep,
"Time")
433 <<
"State(fTimeStepReachedLimit)= " <<
State(fTimeStepReachedLimit) <<
G4endl
434 <<
"State(fGeometryLimitedStep)=" <<
State(fGeometryLimitedStep) <<
G4endl
435 <<
"End point distance was: " <<
G4BestUnit(
State(fEndPointDistance),
"Length")
460 <<
" trackID : " << track.
GetTrackID() <<
" Molecule name: "
462 G4cout <<
"Diffusion length : "
465 <<
"\t Current global time : "
486 <<
"G4DNABrownianTransportation::Diffusion :" << setw(8)
488 <<
"\t" <<
" Global Time = "
509 if (waterDensity == 0.0)
523 G4cout <<
"A track is outside water material : trackID = "
540 mem_diff = mem_intermediaire-mem_first;
541 G4cout <<
"\t\t\t >> || MEM || In G4DNABrownianTransportation::Diffusion "
542 "after dealing with waterDensity for "<< track.
GetTrackID()
543 <<
", diff is : " << mem_diff <<
G4endl;
547 State(fMomentumChanged) =
true;
552 mem_diff = mem_intermediaire-mem_first;
553 G4cout <<
"\t\t\t >> || MEM || In G4DNABrownianTransportation::"
554 "After proposing new direction to fParticleChange for "
590 G4cout <<
" G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength - track ID: "
599 track, previousStepSize, currentMinimumStep, currentSafety,
602 if(geometryStepLength==0)
605 if(
State(fGeometryLimitedStep))
611 State(fCurrentTouchableHandle)->GetHistory());
623 State(fCurrentTouchableHandle) = newTouchable;
638 track, previousStepSize, currentMinimumStep, currentSafety,
680 State(fComputeLastPosition) =
false;
681 State(fTimeStepReachedLimit) =
false;
683 if (
State(fGeometryLimitedStep))
694 / (diffusionCoefficient);
699 / (diffusionCoefficient);
706 State(fComputeLastPosition) =
true;
714 * pow(geometryStepLength /
InvErfc(
State(fRandomNumber)),2);
716 State(fTransportEndPosition) = geometryStepLength*
731 State(fTimeStepReachedLimit) =
true;
732 State(fComputeLastPosition) =
true;
738 State(fTimeStepReachedLimit) =
true;
742 State(fComputeLastPosition) =
true;
746 State(fCandidateEndGlobalTime) =
749 State(fEndGlobalTimeComputed) =
true;
751 State(fPathLengthWasCorrected) =
false;
756 geometryStepLength = 2
759 State(fPathLengthWasCorrected) =
true;
761 State(fTransportEndPosition) = geometryStepLength*
770 <<
"G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength = "
782 return geometryStepLength;
794 MemStat mem_first, mem_second, mem_diff;
801 if (
GetIT(track)->GetTrackingInfo()->IsLeadingStep()
802 &&
State(fComputeLastPosition))
815 spaceStep =
State(fEndPointDistance);
816 State(fGeometryLimitedStep) =
true;
821 GetDiffusionCoefficient();
828 spaceStep = sqrt(x * x + y * y + z * z);
830 if(spaceStep >=
State(fEndPointDistance))
832 State(fGeometryLimitedStep) =
true;
836 && spaceStep >=
State(fEndPointDistance))
838 spaceStep =
State(fEndPointDistance);
843 State(fGeometryLimitedStep) =
false;
859 <<
"G4DNABrownianTransportation::AlongStepDoIt: "
860 "GeometryLimitedStep = "
861 <<
State(fGeometryLimitedStep)
875 mem_diff = mem_intermediaire-mem_first;
876 G4cout <<
"\t\t\t >> || MEM || After calling G4ITTransportation::"
877 "AlongStepDoIt for "<< track.
GetTrackID() <<
", diff is : "
926 mem_diff = mem_intermediaire-mem_first;
927 G4cout <<
"\t\t\t >> || MEM || After calling G4DNABrownianTransportation::"
928 "Diffusion for "<< track.
GetTrackID() <<
", diff is : "
virtual void UpdateYourself(G4VPhysicalVolume *pPhysVol, const G4NavigationHistory *history=0)
virtual void Transport(const G4Track &, G4ParticleChangeForTransport &)=0
ThreeVector shoot(const G4int Ap, const G4int Af)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4ParticleChangeForTransport fParticleChange
G4DNABrownianTransportation & operator=(const G4DNABrownianTransportation &other)
std::ostringstream G4ExceptionDescription
static G4VScheduler * Instance()
G4double GetStepLength() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &)
G4double CheckNextStep(const G4ThreeVector &position, const G4ThreeVector &direction, const G4double currentMaxStep, G4double &newSafety)
const G4ThreeVector & GetPosition() const
virtual void LoadTrackState(G4TrackStateManager &manager)
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
G4double GetDiffusionCoefficient() const
Returns the diffusion coefficient D.
virtual void StartTracking(G4Track *aTrack)
G4ThreeVector G4RandomDirection()
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
virtual const G4String & GetName() const =0
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double, G4double currentMinimumStep, G4double ¤tSafety, G4GPILSelection *selection)
const G4Step * GetStep() const
static double InvErfc(double x)
static G4NistManager * Instance()
const G4String & GetParticleName() const
G4double fInternalMinTimeStep
G4ITNavigator * fLinearNavigator
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
G4StepPoint * GetPreStepPoint() const
const G4String & GetName() const
Returns the name of the molecule.
G4bool fUseMaximumTimeBeforeReachingBoundary
static double Erfc(double x)
virtual void StartTracking(G4Track *aTrack)
G4double ComputeGeomLimit(const G4Track &track, G4double &presafety, G4double limit)
G4IT * GetIT(const G4Track *track)
G4GLOB_DLL std::ostream G4cout
G4bool fPathLengthWasCorrected
G4int GetCurrentStepNumber() const
virtual void ResetTrackState()
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
G4bool fComputeLastPosition
G4bool fTimeStepReachedLimit
void SetInstantiateProcessState(G4bool flag)
void SetProcessSubType(G4int)
G4shared_ptr< G4ProcessState > fpState
G4double GetDeltaTime() const
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
const G4String & GetProcessName() const
G4ITSafetyHelper * fpSafetyHelper
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4TrackingInformation * GetTrackingInfo()
static G4DNAMolecularMaterial * Instance()
const G4VProcess * GetProcessDefinedStep() const
const G4ThreeVector & GetMomentumDirection() const
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
virtual double GetLimitingTimeStep() const
G4StepPoint * GetPostStepPoint() const
virtual void ComputeStep(const G4Track &, const G4Step &, const double, double &)
G4bool fUseSchedulerMinTimeSteps
void ProposeEnergy(G4double finalEnergy)
G4VPhysicalVolume * GetVolume() const
void Diffusion(const G4Track &track)
G4double GetGlobalTime() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
Class Description The dynamic molecule holds all the data that change for a molecule It has a pointer...
static const double picosecond
static double InvErf(double x)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeTrackStatus(G4TrackStatus status)
virtual ~G4DNABrownianTransportation()
G4bool ProposesTimeStep() const
G4VITProcess inherits from G4VProcess.
void SetMomentumChanged(G4bool b)
G4double * theInteractionTimeLeft
G4int GetProcessSubType() const
G4BrownianAction * fpBrownianAction
const std::vector< G4double > * fpWaterDensity
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
G4double GetStepLength() const
G4VPhysicalVolume * GetWorldVolume()
G4DNABrownianTransportation(const G4String &aName="DNABrownianTransportation", G4int verbosityLevel=0)