00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00041
00042 #include <CLHEP/Random/Stat.h>
00043
00044 #include "G4DNABrownianTransportation.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "Randomize.hh"
00048 #include "G4Molecule.hh"
00049 #include "G4RandomDirection.hh"
00050 #include "G4ParticleTable.hh"
00051 #include "G4SafetyHelper.hh"
00052 #include "G4TransportationManager.hh"
00053 #include "G4UnitsTable.hh"
00054 #include "G4NistManager.hh"
00055 #include "G4DNAMolecularMaterial.hh"
00056
00057 using namespace std;
00058
00059 #ifndef State
00060 #define State(theXInfo) (fpBrownianState->theXInfo)
00061 #endif
00062
00063
00064
00065
00066
00067
00068
00069
00070 #define GREEN_ON_BLUE "\033[1;32;44m"
00071 #define RESET "\033[0m"
00072
00073 static double InvErf(double x)
00074 {
00075 return CLHEP::HepStat::inverseErf(x);
00076 }
00077
00078 static double InvErfc(double x)
00079 {
00080 return CLHEP::HepStat::inverseErf(1.-x);
00081 }
00082
00083 G4DNABrownianTransportation::G4DNABrownianTransportation(const G4String& aName, G4int verbosity) :
00084 G4ITTransportation(aName, verbosity),
00085 InitProcessState(fpBrownianState, fTransportationState)
00086 {
00087
00088 SetProcessSubType(61);
00089 verboseLevel = 1;
00090 fUseMaximumTimeBeforeReachingBoundary = true;
00091 fNistWater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00092 fpWaterDensity = 0;
00093 }
00094
00095 G4DNABrownianTransportation::~G4DNABrownianTransportation()
00096 {;}
00097
00098 G4DNABrownianTransportation::G4DNABrownianTransportation(const G4DNABrownianTransportation& right) :
00099 G4ITTransportation(right),
00100 InitProcessState(fpBrownianState, fTransportationState)
00101 {
00102
00103 SetProcessSubType(61);
00104 verboseLevel = right.verboseLevel;
00105 fUseMaximumTimeBeforeReachingBoundary = right.fUseMaximumTimeBeforeReachingBoundary;
00106 fNistWater = right.fNistWater;
00107 fpWaterDensity = right.fpWaterDensity ;
00108 }
00109
00110 G4DNABrownianTransportation& G4DNABrownianTransportation::operator=(const G4DNABrownianTransportation& rhs)
00111 {
00112 if (this == &rhs) return *this;
00113
00114 return *this;
00115 }
00116
00117 G4DNABrownianTransportation::G4ITBrownianState::G4ITBrownianState() : G4ITTransportationState()
00118 {
00119 fPathLengthWasCorrected = false;
00120 }
00121
00122 void G4DNABrownianTransportation::StartTracking(G4Track* track)
00123 {
00124 fpState = new G4ITBrownianState();
00125 SetInstantiateProcessState(false);
00126 G4ITTransportation::StartTracking(track);
00127 }
00128
00129 void G4DNABrownianTransportation::BuildPhysicsTable(const G4ParticleDefinition& particle)
00130 {
00131 if(verboseLevel > 0)
00132 {
00133 G4cout << G4endl << GetProcessName() << ": for "
00134 << setw(24) << particle.GetParticleName()
00135 << "\tSubType= " << GetProcessSubType() << G4endl;
00136 }
00137
00138
00139 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetDensityTableFor(G4Material::GetMaterial("G4_WATER"));
00140 }
00141
00142 void G4DNABrownianTransportation::ComputeStep(const G4Track& track,
00143 const G4Step& step,
00144 const double timeStep,
00145 double& spaceStep)
00146 {
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 if(GetIT(track)->GetTrackingInfo()->IsLeadingStep())
00160 {
00161 const G4VITProcess* ITProc = ((const G4VITProcess*) step.GetPostStepPoint()->GetProcessDefinedStep());
00162 bool makeException = true;
00163
00164 if(ITProc && ITProc->ProposesTimeStep()) makeException = false;
00165
00166 if(makeException)
00167 {
00168
00169 G4ExceptionDescription exceptionDescription ;
00170 exceptionDescription << "ComputeStep is called while the track has the minimum interaction time";
00171 exceptionDescription << " so it should not recompute a timeStep ";
00172 G4Exception("G4DNABrownianTransportation::ComputeStep","G4DNABrownianTransportation001",
00173 FatalErrorInArgument,exceptionDescription);
00174 }
00175 }
00176
00177 State(fGeometryLimitedStep) = false;
00178
00179
00180
00181 G4Molecule* molecule = GetMolecule(track) ;
00182 G4double diffCoeff = molecule->GetDiffusionCoefficient();
00183
00184 if(timeStep > 0)
00185 {
00186 spaceStep = DBL_MAX;
00187
00188 while(spaceStep > State(endpointDistance))
00189
00190
00191
00192 {
00193 G4double x = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
00194 G4double y = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
00195 G4double z = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
00196
00197 spaceStep = sqrt(x*x + y*y + z*z);
00198 }
00199
00200
00201
00202
00203 State(fTransportEndPosition)= spaceStep*step.GetPostStepPoint()->GetMomentumDirection() + track.GetPosition();
00204 }
00205 else
00206 {
00207 spaceStep = 0. ;
00208 State(fTransportEndPosition) = track.GetPosition() ;
00209 }
00210
00211 State(fCandidateEndGlobalTime) = step.GetPreStepPoint()->GetGlobalTime() + timeStep ;
00212 State(fEndGlobalTimeComputed) = true ;
00213
00214 #ifdef G4VERBOSE
00215
00216 if(fVerboseLevel>1)
00217 {
00218 G4cout<< GREEN_ON_BLUE
00219 << "G4ITBrownianTransportation::ComputeStep() : "
00220 << " trackID : " << track.GetTrackID()
00221 << " : Molecule name: " << molecule-> GetName()
00222 << G4endl
00223 << "Diffusion length : " << G4BestUnit(spaceStep, "Length")
00224 << " within time step : " << G4BestUnit(timeStep,"Time")
00225 << RESET
00226 << G4endl<< G4endl;
00227 }
00228 #endif
00229 }
00230
00231 G4VParticleChange* G4DNABrownianTransportation::PostStepDoIt( const G4Track& track, const G4Step& step)
00232 {
00233 G4ITTransportation::PostStepDoIt(track,step);
00234
00235 #ifdef G4VERBOSE
00236
00237 if(fVerboseLevel>1)
00238 {
00239 G4cout<< GREEN_ON_BLUE
00240 << "G4ITBrownianTransportation::PostStepDoIt() :"
00241 << " trackID : " << track.GetTrackID()
00242 << " Molecule name: " << GetMolecule(track)-> GetName() << G4endl;
00243 G4cout<< "Diffusion length : "<< G4BestUnit(step.GetStepLength(),"Length") <<" within time step : "
00244 << G4BestUnit(step.GetDeltaTime(),"Time") << "\t"
00245 << " Current global time : " << G4BestUnit(track.GetGlobalTime(),"Time")
00246 << RESET
00247 << G4endl<< G4endl;
00248 }
00249 #endif
00250 return &fParticleChange ;
00251 }
00252
00253 void G4DNABrownianTransportation::Diffusion(
00254 const G4Track& track)
00255 {
00256
00257 #ifdef G4VERBOSE
00258
00259 if (fVerboseLevel>1)
00260 {
00261 G4cout<< GREEN_ON_BLUE
00262 << setw(18)<< "G4DNABrownianTransportation::Diffusion :"
00263 << setw(8) << GetIT(track)->GetName()
00264 << "\t trackID:" << track.GetTrackID() <<"\t"
00265 << " Global Time = " << G4BestUnit(track.GetGlobalTime(),"Time")
00266 << RESET
00267 << G4endl<< G4endl;
00268 }
00269 #endif
00270
00271 G4Material* material = track.GetMaterial();
00272
00273
00274 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00275
00276 if(waterDensity == 0.0)
00277
00278 {
00279 G4cout << "A track is outside water material : trackID"<< track.GetTrackID() << " (" << GetMolecule(track)->GetName() <<")" << G4endl;
00280 G4cout << "Local Time : "<< (track.GetLocalTime()) /s<<G4endl;
00281 G4cout<< "Step Number :" << track.GetCurrentStepNumber() <<G4endl;
00282
00283 fParticleChange.ProposeEnergy(0.) ;
00284 fParticleChange.ProposeTrackStatus(fStopButAlive);
00285
00286
00287
00288
00289
00290 return ;
00291 }
00292
00293 G4double costheta = (2*G4UniformRand()-1);
00294 G4double theta = acos (costheta);
00295 G4double phi = 2*pi*G4UniformRand();
00296
00297 G4double xMomentum = cos(phi)* sin(theta);
00298 G4double yMomentum = sin(theta)*sin(phi);
00299 G4double zMomentum = costheta;
00300
00301 fParticleChange.ProposeMomentumDirection(xMomentum, yMomentum, zMomentum);
00302 State(fMomentumChanged) = true;
00303 fParticleChange.SetMomentumChanged(true) ;
00304
00305
00306
00307
00308
00309
00310 return;
00311 }
00312
00313
00314 G4double G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength(
00315 const G4Track& track,
00316 G4double previousStepSize,
00317 G4double currentMinimumStep,
00318 G4double& currentSafety,
00319 G4GPILSelection* selection)
00320 {
00321 G4double geometryStepLength = G4ITTransportation::AlongStepGetPhysicalInteractionLength(track,
00322 previousStepSize,
00323 currentMinimumStep,
00324 currentSafety,
00325 selection);
00326
00327 G4double diffusionCoefficient = GetMolecule(track)->GetDiffusionCoefficient();
00328
00329
00330 if(State(fGeometryLimitedStep))
00331 {
00332
00333
00334
00335
00336 if(fUseMaximumTimeBeforeReachingBoundary)
00337 {
00338 State(theInteractionTimeLeft) = (geometryStepLength*geometryStepLength)/(64 * diffusionCoefficient);
00339 }
00340 else
00341 {
00342 State(theInteractionTimeLeft) = 1/(4*diffusionCoefficient) * pow(geometryStepLength/InvErfc(G4UniformRand()),2);
00343 }
00344
00345 State(fCandidateEndGlobalTime) = track.GetGlobalTime() + State(theInteractionTimeLeft);
00346 State(fPathLengthWasCorrected) = false;
00347 }
00348 else
00349 {
00350 geometryStepLength = 2*sqrt(diffusionCoefficient*State(theInteractionTimeLeft) ) *InvErf(G4UniformRand());
00351 State(fPathLengthWasCorrected) = true;
00352 State(endpointDistance) = geometryStepLength;
00353 }
00354
00355 return geometryStepLength ;
00356 }
00357
00359
00360
00361
00362 G4VParticleChange* G4DNABrownianTransportation::AlongStepDoIt( const G4Track& track,
00363 const G4Step& step )
00364 {
00365 G4ITTransportation::AlongStepDoIt(track,step);
00366 Diffusion(track);
00367 return &fParticleChange;
00368 }