ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNABrownianTransportation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNABrownianTransportation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
28 //
29 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear in the next releases.
31 //
32 // History:
33 // -----------
34 // 10 Oct 2011 M.Karamitros created
35 //
36 // -------------------------------------------------------------------
37 
40 
41 #include <CLHEP/Random/Stat.h>
42 
43 #include "G4DNABrownianTransportation.hh"
44 
45 #include <G4Scheduler.hh>
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh"
49 #include "G4Molecule.hh"
50 #include "G4RandomDirection.hh"
51 #include "G4ParticleTable.hh"
52 #include "G4SafetyHelper.hh"
53 #include "G4TransportationManager.hh"
54 #include "G4UnitsTable.hh"
55 #include "G4NistManager.hh"
56 #include "G4DNAMolecularMaterial.hh"
57 #include "G4ITNavigator.hh"
58 #include "G4ITSafetyHelper.hh" // Not used yet
59 #include "G4TrackingInformation.hh"
60 
61 using namespace std;
62 
63 #ifndef State
64 #define State(theXInfo) (GetState<G4ITBrownianState>()->theXInfo)
65 #endif
66 
67 //#ifndef State
68 //#define State(theXInfo) (fTransportationState->theXInfo)
69 //#endif
70 
71 //#define USE_COLOR 1
72 
73 #ifdef USE_COLOR
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"
79 #else
80 #define RED ""
81 #define LIGHT_RED ""
82 #define GREEN ""
83 #define GREEN_ON_BLUE ""
84 #define RESET_COLOR ""
85 #endif
86 
87 //#define DEBUG_MEM 1
88 
89 #ifdef DEBUG_MEM
90 #include "G4MemStat.hh"
91 using namespace G4MemStat;
92 using G4MemStat::MemStat;
93 #endif
94 
95 static double InvErf(double x)
96 {
97  return CLHEP::HepStat::inverseErf(x);
98 }
99 
100 static double InvErfc(double x)
101 {
102  return CLHEP::HepStat::inverseErf(1. - x);
103 }
104 
105 //static double Erf(double x)
106 //{
107 // return CLHEP::HepStat::erf(x);
108 //}
109 
110 static double Erfc(double x)
111 {
112  return 1 - CLHEP::HepStat::erf(1. - x);
113 }
114 
115 #ifndef State
116 #define State(theXInfo) (GetState<G4ITTransportationState>()->theXInfo)
117 #endif
118 
119 G4DNABrownianTransportation::G4DNABrownianTransportation(const G4String& aName,
120  G4int verbosity) :
121  G4ITTransportation(aName, verbosity)
122 {
123 
124  fVerboseLevel = 0;
125 
126  fpState.reset(new G4ITBrownianState());
127 
128  //ctor
129  SetProcessSubType(61);
130 
131  fNistWater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
132  fpWaterDensity = 0;
133 
134  fUseMaximumTimeBeforeReachingBoundary = true;
135  fUseSchedulerMinTimeSteps = false;
136  fSpeedMeUp = true;
137 
138  fInternalMinTimeStep = 1*picosecond;
139  fpBrownianAction = 0;
140 }
141 
142 G4DNABrownianTransportation::~G4DNABrownianTransportation()
143 {
144  if(fpBrownianAction) delete fpBrownianAction;
145 }
146 
147 G4DNABrownianTransportation::G4DNABrownianTransportation(const G4DNABrownianTransportation& right) :
148  G4ITTransportation(right)
149 {
150  //copy ctor
151  SetProcessSubType(61);
152  fUseMaximumTimeBeforeReachingBoundary = right
153  .fUseMaximumTimeBeforeReachingBoundary;
154  fUseSchedulerMinTimeSteps = right.fUseSchedulerMinTimeSteps;
155  fNistWater = right.fNistWater;
156  fpWaterDensity = right.fpWaterDensity;
157  fInternalMinTimeStep = right.fInternalMinTimeStep;
158  fSpeedMeUp = right.fSpeedMeUp;
159  fpBrownianAction = right.fpBrownianAction;
160 }
161 
162 G4DNABrownianTransportation& G4DNABrownianTransportation::operator=(const G4DNABrownianTransportation& rhs)
163 {
164  if(this == &rhs) return *this; // handle self assignment
165  //assignment operator
166  return *this;
167 }
168 
169 G4DNABrownianTransportation::G4ITBrownianState::G4ITBrownianState() :
170  G4ITTransportationState()
171 {
172  fPathLengthWasCorrected = false;
173  fTimeStepReachedLimit = false;
174  fComputeLastPosition = false;
175  fRandomNumber = -1;
176 }
177 
178 void G4DNABrownianTransportation::StartTracking(G4Track* track)
179 {
180  fpState.reset(new G4ITBrownianState());
181 // G4cout << "G4DNABrownianTransportation::StartTracking : "
182 // "Initialised track State" << G4endl;
183  SetInstantiateProcessState(false);
184  G4ITTransportation::StartTracking(track);
185 }
186 
187 void G4DNABrownianTransportation::BuildPhysicsTable(const G4ParticleDefinition& particle)
188 {
189  if(verboseLevel > 0)
190  {
191  G4cout << G4endl<< GetProcessName() << ": for "
192  << setw(24) << particle.GetParticleName()
193  << "\tSubType= " << GetProcessSubType() << G4endl;
194  }
195  // Initialize water density pointer
196  fpWaterDensity = G4DNAMolecularMaterial::Instance()->
197  GetDensityTableFor(G4Material::GetMaterial("G4_WATER"));
198 
199  fpSafetyHelper->InitialiseHelper();
200  G4ITTransportation::BuildPhysicsTable(particle);
201 }
202 
203 void G4DNABrownianTransportation::ComputeStep(const G4Track& track,
204  const G4Step& step,
205  const double timeStep,
206  double& spaceStep)
207 {
208  // G4cout << "G4ITBrownianTransportation::ComputeStep" << G4endl;
209 
210  /* If this method is called, this step
211  * cannot be geometry limited.
212  * In case the step is limited by the geometry,
213  * this method should not be called.
214  * The fTransportEndPosition calculated in
215  * the method AlongStepIL should be taken
216  * into account.
217  * In order to do so, the flag IsLeadingStep
218  * is on. Meaning : this track has the minimum
219  * interaction length over all others.
220  */
221  if(GetIT(track)->GetTrackingInfo()->IsLeadingStep())
222  {
223  const G4VITProcess* ITProc = ((const G4VITProcess*) step.GetPostStepPoint()
224  ->GetProcessDefinedStep());
225  bool makeException = true;
226 
227  if(ITProc && ITProc->ProposesTimeStep()) makeException = false;
228 
229  if(makeException)
230  {
231 
232  G4ExceptionDescription exceptionDescription;
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",
237  "G4DNABrownianTransportation001", FatalErrorInArgument,
238  exceptionDescription);
239  }
240  }
241 
242  State(fGeometryLimitedStep) = false;
243 
244  G4Molecule* molecule = GetMolecule(track);
245 
246  if(timeStep > 0)
247  {
248  spaceStep = DBL_MAX;
249 
250  // TODO : generalize this process to all kind of Brownian objects
251  G4double diffCoeff = molecule->GetDiffusionCoefficient(track.GetMaterial(),
252  track.GetMaterial()->GetTemperature());
253 
254  static double sqrt_2 = sqrt(2.);
255  double sqrt_Dt = sqrt(diffCoeff*timeStep);
256  double sqrt_2Dt = sqrt_2*sqrt_Dt;
257 
258  if(State(fTimeStepReachedLimit)== true)
259  {
260  //========================================================================
261  State(fGeometryLimitedStep) = true;// important
262  //========================================================================
263  spaceStep = State(fEndPointDistance);
264  // G4cout << "State(fTimeStepReachedLimit)== true" << G4endl;
265  }
266  else
267  {
268  G4double x = G4RandGauss::shoot(0,sqrt_2Dt);
269  G4double y = G4RandGauss::shoot(0,sqrt_2Dt);
270  G4double z = G4RandGauss::shoot(0,sqrt_2Dt);
271 
272  spaceStep = sqrt(x*x + y*y + z*z);
273 
274  if(spaceStep >= State(fEndPointDistance))
275  {
276  //G4cout << "spaceStep >= State(fEndPointDistance)" << G4endl;
277  //======================================================================
278  State(fGeometryLimitedStep) = true;// important
279  //======================================================================
280 /*
281  if(fSpeedMeUp)
282  {
283  G4cout << "SpeedMeUp" << G4endl;
284  }
285  else
286 */
287  if(fUseSchedulerMinTimeSteps == false)// jump over barrier NOT used
288  {
289 #ifdef G4VERBOSE
290  if (fVerboseLevel > 1)
291  {
292  G4cout << GREEN_ON_BLUE
293  << "G4ITBrownianTransportation::ComputeStep() : "
294  << "Step was limited to boundary"
295  << RESET_COLOR
296  << G4endl;
297  }
298 #endif
299  //TODO
300  if(State(fRandomNumber)>=0) // CDF is used
301  {
302  /*
303  //=================================================================
304  State(fGeometryLimitedStep) = true;// important
305  //=================================================================
306  spaceStep = State(fEndPointDistance);
307  */
308 
309  //==================================================================
310  // BE AWARE THAT THE TECHNIQUE USED BELOW IS A 1D APPROXIMATION
311  // Cumulative density function for the 3D case is not yet
312  // implemented
313  //==================================================================
314 // G4cout << GREEN_ON_BLUE
315 // << "G4ITBrownianTransportation::ComputeStep() : "
316 // << "A random number was selected"
317 // << RESET_COLOR
318 // << G4endl;
319  double value = State(fRandomNumber)+(1-State(fRandomNumber))*G4UniformRand();
320  double invErfc = InvErfc(value);
321  spaceStep = invErfc*2*sqrt_Dt;
322 
323  if(State(fTimeStepReachedLimit)== false)
324  {
325  //================================================================
326  State(fGeometryLimitedStep) = false;// important
327  //================================================================
328  }
329  //==================================================================
330  // DEBUG
331 // if(spaceStep > State(fEndPointDistance))
332 // {
333 // G4cout << "value = " << value << G4endl;
334 // G4cout << "invErfc = " << invErfc << G4endl;
335 // G4cout << "spaceStep = " << G4BestUnit(spaceStep, "Length")
336 // << G4endl;
337 // G4cout << "end point distance= " << G4BestUnit(State(fEndPointDistance), "Length")
338 // << G4endl;
339 // }
340 //
341 // assert(spaceStep <= State(fEndPointDistance));
342  //==================================================================
343 
344  }
345  else if(fUseMaximumTimeBeforeReachingBoundary == false) // CDF is used
346  {
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))
352  {
353  //================================================================
354  State(fGeometryLimitedStep) = true;// important
355  //================================================================
356  }
357  else if(State(fTimeStepReachedLimit)== false)
358  {
359  //================================================================
360  State(fGeometryLimitedStep) = false;// important
361  //================================================================
362  }
363  }
364  else // CDF is NOT used
365  {
366  //==================================================================
367  State(fGeometryLimitedStep) = true;// important
368  //==================================================================
369  spaceStep = State(fEndPointDistance);
370  //TODO
371 
372  /*
373  //==================================================================
374  // 1D approximation to place the brownian between its starting point
375  // and the geometry boundary
376  //==================================================================
377  double min_randomNumber = Erfc(State(fEndPointDistance)/2*sqrt_Dt);
378  double value = State(fRandomNumber)+(1-State(fRandomNumber))*G4UniformRand();
379  double invErfc = InvErfc(value*G4UniformRand());
380  spaceStep = invErfc*2*sqrt_Dt;
381  State(fGeometryLimitedStep) = false;
382  */
383  }
384  }
385 
386  State(fTransportEndPosition)= spaceStep*
387 // step.GetPostStepPoint()->GetMomentumDirection()
388  track.GetMomentumDirection()
389  + track.GetPosition();
390  }
391  else
392  {
393  //======================================================================
394  State(fGeometryLimitedStep) = false;// important
395  //======================================================================
396  State(fTransportEndPosition)= spaceStep*step.GetPostStepPoint()->
397  GetMomentumDirection() + track.GetPosition();
398  }
399  }
400  }
401  else
402  {
403  spaceStep = 0.;
404  State(fTransportEndPosition) = track.GetPosition();
405  State(fGeometryLimitedStep) = false;
406  }
407 
408  State(fCandidateEndGlobalTime) = step.GetPreStepPoint()->GetGlobalTime()
409  + timeStep;
410  State(fEndGlobalTimeComputed) = true;
411 
412 #ifdef G4VERBOSE
413  // DEBUG
414  if (fVerboseLevel > 1)
415  {
416  G4cout << GREEN_ON_BLUE
417  << "G4ITBrownianTransportation::ComputeStep() : "
418  << " trackID : " << track.GetTrackID() << " : Molecule name: "
419  << molecule->GetName() << G4endl
420  << "Initial position:" << G4BestUnit(track.GetPosition(), "Length")
421  << G4endl
422  << "Initial direction:" << track.GetMomentumDirection() << G4endl
423  << "Final position:" << G4BestUnit(State(fTransportEndPosition), "Length")
424  << G4endl
425  << "Initial magnitude:" << G4BestUnit(track.GetPosition().mag(), "Length")
426  << G4endl
427  << "Final magnitude:" << G4BestUnit(State(fTransportEndPosition).mag(), "Length")
428  << G4endl
429  << "Diffusion length : "
430  << G4BestUnit(spaceStep, "Length")
431  << " within time step : " << G4BestUnit(timeStep,"Time")
432  << G4endl
433  << "State(fTimeStepReachedLimit)= " << State(fTimeStepReachedLimit) << G4endl
434  << "State(fGeometryLimitedStep)=" << State(fGeometryLimitedStep) << G4endl
435  << "End point distance was: " << G4BestUnit(State(fEndPointDistance), "Length")
436  << G4endl
437  << RESET_COLOR
438  << G4endl<< G4endl;
439  }
440 #endif
441 
442 //==============================================================================
443 // DEBUG
444 //assert(spaceStep < State(fEndPointDistance)
445 // || (spaceStep >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
446 //assert(track.GetMomentumDirection() == State(fTransportEndMomentumDir));
447 //==============================================================================
448 }
449 
450 G4VParticleChange* G4DNABrownianTransportation::PostStepDoIt(const G4Track& track,
451  const G4Step& step)
452 {
453  G4ITTransportation::PostStepDoIt(track, step);
454 
455 #ifdef G4VERBOSE
456  // DEBUG
457  if (fVerboseLevel > 1)
458  {
459  G4cout << GREEN_ON_BLUE << "G4ITBrownianTransportation::PostStepDoIt() :"
460  << " trackID : " << track.GetTrackID() << " Molecule name: "
461  << GetMolecule(track)->GetName() << G4endl;
462  G4cout << "Diffusion length : "
463  << G4BestUnit(step.GetStepLength(), "Length")
464  <<" within time step : " << G4BestUnit(step.GetDeltaTime(),"Time")
465  << "\t Current global time : "
466  << G4BestUnit(track.GetGlobalTime(),"Time")
467  << RESET_COLOR
468  << G4endl<< G4endl;
469  }
470 #endif
471  return &fParticleChange;
472 }
473 
474 void G4DNABrownianTransportation::Diffusion(const G4Track& track)
475 {
476 
477 #ifdef DEBUG_MEM
478  MemStat mem_first = MemoryUsage();
479 #endif
480 
481 #ifdef G4VERBOSE
482  // DEBUG
483  if (fVerboseLevel > 1)
484  {
485  G4cout << GREEN_ON_BLUE << setw(18)
486  << "G4DNABrownianTransportation::Diffusion :" << setw(8)
487  << GetIT(track)->GetName() << "\t trackID:" << track.GetTrackID()
488  << "\t" << " Global Time = "
489  << G4BestUnit(track.GetGlobalTime(), "Time")
490  << RESET_COLOR
491  << G4endl
492  << G4endl;
493  }
494 #endif
495 
496 /*
497  fParticleChange.ProposePosition(State(fTransportEndPosition));
498  //fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy));
499  fParticleChange.SetMomentumChanged(State(fMomentumChanged));
500 
501  fParticleChange.ProposeGlobalTime(State(fCandidateEndGlobalTime));
502  fParticleChange.ProposeLocalTime(State(fCandidateEndGlobalTime));
503  fParticleChange.ProposeTrueStepLength(track.GetStepLength());
504 */
505  G4Material* material = track.GetMaterial();
506 
507  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
508 
509  if (waterDensity == 0.0)
510  {
511  if(fpBrownianAction)
512  {
513  // Let the user Brownian action class decide what to do
514  fpBrownianAction->Transport(track,
515  fParticleChange);
516  return;
517  }
518  else
519  {
520 #ifdef G4VERBOSE
521  if(fVerboseLevel)
522  {
523  G4cout << "A track is outside water material : trackID = "
524  << track.GetTrackID() << " (" << GetMolecule(track)->GetName() <<")"
525  << G4endl;
526  G4cout << "Local Time : " << G4BestUnit(track.GetGlobalTime(), "Time")
527  << G4endl;
528  G4cout << "Step Number :" << track.GetCurrentStepNumber() << G4endl;
529  }
530 #endif
531  fParticleChange.ProposeEnergy(0.);
532  fParticleChange.ProposeTrackStatus(fStopAndKill);
533  return;// &fParticleChange is the final returned object
534  }
535  }
536 
537 
538  #ifdef DEBUG_MEM
539  MemStat mem_intermediaire = MemoryUsage();
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;
544  #endif
545 
546  fParticleChange.ProposeMomentumDirection(G4RandomDirection());
547  State(fMomentumChanged) = true;
548  fParticleChange.SetMomentumChanged(true);
549  //
550  #ifdef DEBUG_MEM
551  mem_intermediaire = MemoryUsage();
552  mem_diff = mem_intermediaire-mem_first;
553  G4cout << "\t\t\t >> || MEM || In G4DNABrownianTransportation::"
554  "After proposing new direction to fParticleChange for "
555  << track.GetTrackID() << ", diff is : " << mem_diff << G4endl;
556  #endif
557 
558  return;// &fParticleChange is the final returned object
559 }
560 
561 // NOT USED
562 G4double G4DNABrownianTransportation::ComputeGeomLimit(const G4Track& track,
563  G4double& presafety,
564  G4double limit)
565 {
566  G4double res = DBL_MAX;
567  if(track.GetVolume() != fpSafetyHelper->GetWorldVolume())
568  {
569  G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
570  ->GetTrackStateManager();
571  fpSafetyHelper->LoadTrackState(trackStateMan);
572  res = fpSafetyHelper->CheckNextStep(
573  track.GetStep()->GetPreStepPoint()->GetPosition(),
574  track.GetMomentumDirection(),
575  limit, presafety);
576  fpSafetyHelper->ResetTrackState();
577  }
578  return res;
579 }
580 
581 G4double G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength(const G4Track& track,
582  G4double previousStepSize,
583  G4double currentMinimumStep,
584  G4double& currentSafety,
585  G4GPILSelection* selection)
586 {
587 #ifdef G4VERBOSE
588  if(fVerboseLevel)
589  {
590  G4cout << " G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength - track ID: "
591  << track.GetTrackID() << G4endl;
592  G4cout << "In volume : " << track.GetVolume()->GetName()
593  << " position : " << G4BestUnit(track.GetPosition() , "Length") << G4endl;
594  }
595 #endif
596 
597  G4double geometryStepLength =
598  G4ITTransportation::AlongStepGetPhysicalInteractionLength(
599  track, previousStepSize, currentMinimumStep, currentSafety,
600  selection);
601 
602  if(geometryStepLength==0)
603  {
604 // G4cout << "geometryStepLength==0" << G4endl;
605  if(State(fGeometryLimitedStep))
606  {
607 // G4cout << "if(State(fGeometryLimitedStep))" << G4endl;
608  G4TouchableHandle newTouchable = new G4TouchableHistory;
609 
610  newTouchable->UpdateYourself(State(fCurrentTouchableHandle)->GetVolume(),
611  State(fCurrentTouchableHandle)->GetHistory());
612 
613  fLinearNavigator->SetGeometricallyLimitedStep();
614  fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
615  track.GetPosition(), track.GetMomentumDirection(),
616  newTouchable, true);
617 
618  if(newTouchable->GetVolume() == 0)
619  {
620  return 0;
621  }
622 
623  State(fCurrentTouchableHandle) = newTouchable;
624 
625  //=======================================================================
626  // TODO: speed up navigation update
627 // geometryStepLength = fLinearNavigator->ComputeStep(track.GetPosition(),
628 // track.GetMomentumDirection(),
629 // currentMinimumStep,
630 // currentSafety);
631  //=======================================================================
632 
633 
634  //=======================================
635  // Longer but safer ...
636  geometryStepLength =
637  G4ITTransportation::AlongStepGetPhysicalInteractionLength(
638  track, previousStepSize, currentMinimumStep, currentSafety,
639  selection);
640 
641  }
642  }
643 
644  //============================================================================
645  // DEBUG
646  // G4cout << "geometryStepLength: " << G4BestUnit(geometryStepLength, "Length")
647  // << " | trackID: " << track.GetTrackID()
648  // << G4endl;
649  //============================================================================
650 
651  G4double diffusionCoefficient = 0;
652 
653  /*
654  G4Material* material = track.GetMaterial();
655  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
656 
657  if (waterDensity == 0.0)
658  {
659  if(fpBrownianAction)
660  {
661  diffusionCoefficient = fpBrownianAction->GetDiffusionCoefficient(material,
662  GetMolecule(track));
663  }
664 
665  if(diffusionCoefficient <= 0)
666  {
667  State(fGeometryLimitedStep) = false;
668  State(theInteractionTimeLeft) = 0;
669  State(fTransportEndPosition) = track.GetPosition();
670  return 0;
671  }
672 
673  }
674  else
675  {
676  diffusionCoefficient = GetMolecule(track)->GetDiffusionCoefficient();
677  }
678  */
679 
680  diffusionCoefficient = GetMolecule(track)->GetDiffusionCoefficient();
681 
682  // To avoid divide by zero of diffusionCoefficient
683  if(diffusionCoefficient <= 0)
684  {
685  State(fGeometryLimitedStep) = false;
686  State(theInteractionTimeLeft) = DBL_MAX;
687  State(fTransportEndPosition) = track.GetPosition();
688  return 0;
689  }
690 
691 
692  State(fComputeLastPosition) = false;
693  State(fTimeStepReachedLimit) = false;
694 
695  if (State(fGeometryLimitedStep))
696  {
697  // 95 % of the space step distribution is lower than
698  // d_95 = 2 * sqrt(2*D*t)
699  // where t is the corresponding time step
700  // so by inversion :
701  if (fUseMaximumTimeBeforeReachingBoundary)
702  {
703  if(fSpeedMeUp)
704  {
705  State(theInteractionTimeLeft) = (geometryStepLength * geometryStepLength)
706  / (diffusionCoefficient); // d_50 - use straight line
707  }
708  else
709  {
710  State(theInteractionTimeLeft) = (currentSafety * currentSafety)
711  / (diffusionCoefficient); // d_50 - use safety
712 
713  //=====================================================================
714  // State(theInteractionTimeLeft) = (currentSafety * currentSafety)
715  // / (8 * diffusionCoefficient); // d_95
716  //=====================================================================
717  }
718  State(fComputeLastPosition) = true;
719  }
720  else
721  // Will use a random time - this is precise but long to compute in certain
722  // circumstances (many particles - small volumes)
723  {
724  State(fRandomNumber) = G4UniformRand();
725  State(theInteractionTimeLeft) = 1 / (4 * diffusionCoefficient)
726  * pow(geometryStepLength / InvErfc(State(fRandomNumber)),2);
727 
728  State(fTransportEndPosition) = geometryStepLength*
729  track.GetMomentumDirection() + track.GetPosition();
730  }
731 
732  if (fUseSchedulerMinTimeSteps)
733  {
734  double minTimeStepAllowed = G4VScheduler::Instance()->GetLimitingTimeStep();
735  //========================================================================
736  // TODO
737 // double currentMinTimeStep = G4VScheduler::Instance()->GetTimeStep();
738  //========================================================================
739 
740  if (State(theInteractionTimeLeft) < minTimeStepAllowed)
741  {
742  State(theInteractionTimeLeft) = minTimeStepAllowed;
743  State(fTimeStepReachedLimit) = true;
744  State(fComputeLastPosition) = true;
745  }
746  }
747  else if(State(theInteractionTimeLeft) < fInternalMinTimeStep)
748  // TODO: find a better way when fForceLimitOnMinTimeSteps is not used
749  {
750  State(fTimeStepReachedLimit) = true;
751  State(theInteractionTimeLeft) = fInternalMinTimeStep;
752  if (fUseMaximumTimeBeforeReachingBoundary)
753  {
754  State(fComputeLastPosition) = true;
755  }
756  }
757 
758  State(fCandidateEndGlobalTime) =
759  track.GetGlobalTime() + State(theInteractionTimeLeft);
760 
761  State(fEndGlobalTimeComputed) = true; // MK: ADDED ON 20/11/2014
762 
763  State(fPathLengthWasCorrected) = false;
764  }
765  else
766  {
767  // Transform geometrical step
768  geometryStepLength = 2
769  * sqrt(diffusionCoefficient * State(theInteractionTimeLeft))
770  * InvErf(G4UniformRand());
771  State(fPathLengthWasCorrected) = true;
772  //State(fEndPointDistance) = geometryStepLength;
773  State(fTransportEndPosition) = geometryStepLength*
774  track.GetMomentumDirection() + track.GetPosition();
775  }
776 
777 #ifdef G4VERBOSE
778  // DEBUG
779  if (fVerboseLevel > 1)
780  {
781  G4cout << GREEN_ON_BLUE
782  << "G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength = "
783  << G4BestUnit(geometryStepLength, "Length")
784  << " | trackID = "
785  << track.GetTrackID()
786  << RESET_COLOR
787  << G4endl;
788  }
789 #endif
790 
791 // assert(geometryStepLength < State(fEndPointDistance)
792 // || (geometryStepLength >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
793 
794  return geometryStepLength;
795 }
796 
798 //
799 // Initialize ParticleChange (by setting all its members equal
800 // to corresponding members in G4Track)
801 G4VParticleChange*
802 G4DNABrownianTransportation::AlongStepDoIt(const G4Track& track,
803  const G4Step& step)
804 {
805 #ifdef DEBUG_MEM
806  MemStat mem_first, mem_second, mem_diff;
807 #endif
808 
809 #ifdef DEBUG_MEM
810  mem_first = MemoryUsage();
811 #endif
812 
813  if (GetIT(track)->GetTrackingInfo()->IsLeadingStep()
814  && State(fComputeLastPosition))
815  {
816  //==========================================================================
817  // DEBUG
818  //
819 // assert(fabs(State(theInteractionTimeLeft)-
820 // G4VScheduler::Instance()->GetTimeStep()) < DBL_EPSILON);
821  //==========================================================================
822 
823  double spaceStep = DBL_MAX;
824 
825  if(State(theInteractionTimeLeft) <= fInternalMinTimeStep)
826  {
827  spaceStep = State(fEndPointDistance);
828  State(fGeometryLimitedStep) = true;
829  }
830  else
831  {
832  G4double diffusionCoefficient = GetMolecule(track)->
833  GetDiffusionCoefficient();
834 
835  double sqrt_2Dt= sqrt(2 * diffusionCoefficient * State(theInteractionTimeLeft));
836  G4double x = G4RandGauss::shoot(0, sqrt_2Dt);
837  G4double y = G4RandGauss::shoot(0, sqrt_2Dt);
838  G4double z = G4RandGauss::shoot(0, sqrt_2Dt);
839 
840  spaceStep = sqrt(x * x + y * y + z * z);
841 
842  if(spaceStep >= State(fEndPointDistance))
843  {
844  State(fGeometryLimitedStep) = true;
845  if (
846  //fSpeedMeUp == false&&
847  fUseSchedulerMinTimeSteps == false
848  && spaceStep >= State(fEndPointDistance))
849  {
850  spaceStep = State(fEndPointDistance);
851  }
852  }
853  else
854  {
855  State(fGeometryLimitedStep) = false;
856  }
857  }
858 
859 // assert( (spaceStep < State(fEndPointDistance) && State(fGeometryLimitedStep) == false)
860 //|| (spaceStep >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
861 
862  // Calculate final position
863  //
864  State(fTransportEndPosition) = track.GetPosition()
865  + spaceStep * track.GetMomentumDirection();
866  }
867 
868  if(fVerboseLevel)
869  {
870  G4cout << GREEN_ON_BLUE
871  << "G4DNABrownianTransportation::AlongStepDoIt: "
872  "GeometryLimitedStep = "
873  << State(fGeometryLimitedStep)
874  << RESET_COLOR
875  << G4endl;
876  }
877 
878 // static G4ThreadLocal G4int noCalls = 0;
879 // noCalls++;
880 
881 // fParticleChange.Initialize(track);
882 
883  G4ITTransportation::AlongStepDoIt(track, step);
884 
885 #ifdef DEBUG_MEM
886  MemStat mem_intermediaire = MemoryUsage();
887  mem_diff = mem_intermediaire-mem_first;
888  G4cout << "\t\t\t >> || MEM || After calling G4ITTransportation::"
889  "AlongStepDoIt for "<< track.GetTrackID() << ", diff is : "
890  << mem_diff << G4endl;
891 #endif
892 
893  if(track.GetStepLength() != 0 // && State(fGeometryLimitedStep)
894  //========================================================================
895  // TODO: avoid changing direction after too small time steps
896 // && (G4VScheduler::Instance()->GetTimeStep() > fInternalMinTimeStep
897 // || fSpeedMeUp == false)
898  //========================================================================
899  )
900  {
901  Diffusion(track);
902  }
903  //else
904  //{
905  // fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir));
906  //}
907 /*
908  if (State(fParticleIsLooping))
909  {
910  if ((State(fNoLooperTrials)>= fThresholdTrials))
911  {
912  fParticleChange.ProposeTrackStatus(fStopAndKill);
913  State(fNoLooperTrials) = 0;
914 #ifdef G4VERBOSE
915  if ((fVerboseLevel > 1))
916  {
917  G4cout
918  << " G4DNABrownianTransportation is killing track that is looping or stuck "
919  << G4endl;
920  G4cout << " Number of trials = " << State(fNoLooperTrials)
921  << " No of calls to AlongStepDoIt = " << noCalls
922  << G4endl;
923  }
924 #endif
925  }
926  else
927  {
928  State(fNoLooperTrials)++;
929  }
930  }
931  else
932  {
933  State(fNoLooperTrials)=0;
934  }
935 */
936 #ifdef DEBUG_MEM
937  mem_intermediaire = MemoryUsage();
938  mem_diff = mem_intermediaire-mem_first;
939  G4cout << "\t\t\t >> || MEM || After calling G4DNABrownianTransportation::"
940  "Diffusion for "<< track.GetTrackID() << ", diff is : "
941  << mem_diff << G4endl;
942 #endif
943 
944  return &fParticleChange;
945 }
946