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 #include "G4DNASecondOrderReaction.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4Molecule.hh"
00031 #include "G4DNAMolecularMaterial.hh"
00032 #include "G4MolecularConfiguration.hh"
00033 #include "G4DNADamages.hh"
00034 #include "G4UnitsTable.hh"
00035 #include "G4ITTrackHolder.hh"
00036
00037 void G4DNASecondOrderReaction::Create()
00038 {
00039 pParticleChange = &fParticleChange;
00040 enableAtRestDoIt = false;
00041 enableAlongStepDoIt = false;
00042 enablePostStepDoIt = true;
00043
00044 SetProcessSubType(60);
00045
00046 G4VITProcess::SetInstantiateProcessState(false);
00047
00048
00049 fIsInitialized = false;
00050 fpMolecularConfiguration = 0;
00051 fpMaterial = 0;
00052 fReactionRate = -1.;
00053 fConcentration = -1.;
00054 fMolarMassOfMaterial = -1.;
00055 fProposesTimeStep = true;
00056 fReturnedValue = DBL_MAX;
00057 fpMoleculeDensity = 0;
00058
00059 verboseLevel = 0;
00060 }
00061
00062 G4DNASecondOrderReaction::G4DNASecondOrderReaction(const G4String &aName, G4ProcessType type) :
00063 G4VITProcess(aName,type),
00064 InitProcessState(fpSecondOrderReactionState, fpState)
00065 {
00066 Create();
00067 }
00068
00069 G4DNASecondOrderReaction::G4DNASecondOrderReaction(const G4DNASecondOrderReaction& rhs):
00070 G4VITProcess(rhs),
00071 InitProcessState(fpSecondOrderReactionState, fpState)
00072 {
00073 Create();
00074 }
00075
00076 G4DNASecondOrderReaction::~G4DNASecondOrderReaction()
00077 {
00078 ;
00079 }
00080 G4DNASecondOrderReaction& G4DNASecondOrderReaction::operator=(const G4DNASecondOrderReaction& rhs)
00081 {
00082 if (this == &rhs) return *this;
00083
00084
00085 return *this;
00086 }
00087
00088 G4DNASecondOrderReaction::SecondOrderReactionState::SecondOrderReactionState() : G4ProcessState()
00089 {
00090 fPreviousTimeAtPreStepPoint = -1;
00091 fIsInGoodMaterial = false;
00092 }
00093
00094 void G4DNASecondOrderReaction::BuildPhysicsTable(const G4ParticleDefinition&)
00095 {
00096 fpMoleculeDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(fpMaterial);
00097 fMolarMassOfMaterial = fpMaterial->GetMassOfMolecule()*CLHEP::Avogadro*1e3;
00098 fIsInitialized = true;
00099 }
00100
00101 void
00102 G4DNASecondOrderReaction::StartTracking(G4Track* track)
00103 {
00104 G4VProcess::StartTracking(track);
00105 G4VITProcess::fpState = new SecondOrderReactionState();
00106 G4VITProcess::StartTracking(track);
00107 }
00108
00109 void
00110 G4DNASecondOrderReaction::SetReaction(const G4MolecularConfiguration* molConf,
00111 const G4Material* mat, double reactionRate)
00112 {
00113 if(fIsInitialized)
00114 {
00115 G4ExceptionDescription exceptionDescription ;
00116 exceptionDescription << "G4DNASecondOrderReaction was already initialised. ";
00117 exceptionDescription << "You cannot set a reaction after initialisation.";
00118 G4Exception("G4DNASecondOrderReaction::SetReaction","G4DNASecondOrderReaction001",
00119 FatalErrorInArgument,exceptionDescription);
00120 }
00121 fpMolecularConfiguration = molConf;
00122 fpMaterial = mat;
00123 fReactionRate = reactionRate;
00124 }
00125
00126 G4double G4DNASecondOrderReaction::PostStepGetPhysicalInteractionLength(const G4Track& track,
00127 G4double ,
00128 G4ForceCondition* pForceCond)
00129 {
00130
00131
00132
00133
00134
00135 const G4Material* material = track.GetMaterial();
00136
00137 G4Molecule* mol = GetMolecule(track);
00138 if(!mol) return DBL_MAX;
00139 if(mol->GetMolecularConfiguration() != fpMolecularConfiguration)
00140 {
00141
00142 return DBL_MAX;
00143 }
00144
00145 G4double molDensity = (*fpMoleculeDensity)[material->GetIndex()];
00146
00147 if(molDensity == 0.0)
00148 {
00149 if(fpSecondOrderReactionState->fIsInGoodMaterial)
00150 {
00151 ResetNumberOfInteractionLengthLeft();
00152 fpSecondOrderReactionState->fIsInGoodMaterial = false;
00153 }
00154
00155
00156
00157
00158
00159 return DBL_MAX;
00160 }
00161
00162
00163
00164 fpSecondOrderReactionState->fIsInGoodMaterial = true;
00165
00166
00167 fConcentration = molDensity/CLHEP::Avogadro;
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 G4double previousTimeStep(-1.);
00178
00179 if(track.GetCurrentStepNumber() > 0)
00180 previousTimeStep = track.GetGlobalTime() - fpSecondOrderReactionState->fPreviousTimeAtPreStepPoint ;
00181
00182 fpSecondOrderReactionState->fPreviousTimeAtPreStepPoint = track.GetGlobalTime();
00183
00184 if ( (previousTimeStep < 0.0) || (fpState->theNumberOfInteractionLengthLeft<=0.0)) {
00185
00186 ResetNumberOfInteractionLengthLeft();
00187 } else if ( previousTimeStep > 0.0) {
00188
00189 SubtractNumberOfInteractionLengthLeft(previousTimeStep );
00190 } else {
00191
00192
00193 }
00194
00195
00196 *pForceCond = NotForced;
00197
00198
00199 fpState->currentInteractionLength = 1/(fReactionRate*fConcentration);
00200
00201 G4double value;
00202 if (fpState->currentInteractionLength <DBL_MAX) {
00203 value = fpState->theNumberOfInteractionLengthLeft * (fpState->currentInteractionLength);
00204 } else {
00205 value = DBL_MAX;
00206 }
00207 #ifdef G4VERBOSE
00208 if (verboseLevel>2){
00209 G4cout << "G4VITRestDiscreteProcess::PostStepGetPhysicalInteractionLength ";
00210 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
00211 track.GetDynamicParticle()->DumpInfo();
00212 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
00213 G4cout << "InteractionLength= " << value/cm <<"[cm] " <<G4endl;
00214 }
00215 #endif
00216
00217
00218
00219
00220 if(value < fReturnedValue)
00221 fReturnedValue = value;
00222
00223 return value*-1;
00224
00225 }
00226
00227 G4VParticleChange* G4DNASecondOrderReaction::PostStepDoIt(const G4Track& track,const G4Step& )
00228 {
00229 G4Molecule* molecule = GetMolecule(track);
00230 #ifdef G4VERBOSE
00231 if(verboseLevel > 1)
00232 {
00233 G4cout << "___________" << G4endl;
00234 G4cout << ">>> Beginning of G4DNASecondOrderReaction verbose" << G4endl;
00235 G4cout << ">>> Returned value : " << G4BestUnit(fReturnedValue,"Time") << G4endl;
00236 G4cout << ">>> Time Step : " << G4BestUnit(G4ITTrackHolder::Instance()->GetTimeStep(),"Time") << G4endl;
00237 G4cout << ">>> Reaction : " << molecule->GetName() << " + " << fpMaterial->GetName() << G4endl;
00238 G4cout << ">>> End of G4DNASecondOrderReaction verbose <<<" << G4endl;
00239 }
00240 #endif
00241 fReturnedValue = DBL_MAX;
00242 fParticleChange.Initialize(track);
00243 fParticleChange.ProposeTrackStatus(fStopAndKill);
00244 G4DNADamages::Instance()->AddIndirectDamage(fpMaterial->GetName(),molecule,track.GetPosition(),track.GetGlobalTime());
00245 return &fParticleChange;
00246 }
00247