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
00039 #include "G4DNAMolecularDecay.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4Track.hh"
00042 #include "G4Molecule.hh"
00043 #include "G4ITManager.hh"
00044 #include "G4ParticleChange.hh"
00045
00046 using namespace std;
00047
00048 G4DNAMolecularDecay::G4DNAMolecularDecay(const G4String& processName,
00049 G4ProcessType type) : G4VITRestProcess(processName, type)
00050 {
00051
00052 SetProcessSubType(59);
00053 enableAlongStepDoIt = false;
00054 enablePostStepDoIt = false;
00055 enableAtRestDoIt=true;
00056
00057 fVerbose = 0 ;
00058
00059 #ifdef G4VERBOSE
00060 if (verboseLevel>1)
00061 {
00062 G4cout << "G4MolecularDecayProcess constructor " << " Name:" << processName << G4endl;
00063 }
00064 #endif
00065
00066 pParticleChange = &aParticleChange;
00067
00068 fDecayAtFixedTime = true ;
00069 }
00070
00071 G4DNAMolecularDecay::~G4DNAMolecularDecay()
00072 {
00073 DecayDisplacementMap::iterator it = fDecayDisplacementMap.begin();
00074
00075 for( ; it != fDecayDisplacementMap.end() ; it++)
00076 {
00077 if(it->second)
00078 {
00079 delete it->second;
00080 it->second = 0;
00081 }
00082 }
00083 fDecayDisplacementMap.clear();
00084 }
00085
00086 G4DNAMolecularDecay::G4DNAMolecularDecay(const G4DNAMolecularDecay &right) :
00087 G4VITRestProcess(right)
00088 {
00089 fDecayAtFixedTime = right . fDecayAtFixedTime;
00090 fDecayDisplacementMap = right.fDecayDisplacementMap;
00091 fVerbose = right.fVerbose ;
00092 }
00093
00094 G4bool G4DNAMolecularDecay::IsApplicable(const G4ParticleDefinition& aParticleType)
00095 {
00096 if(aParticleType.GetParticleType()=="Molecule")
00097 {
00098 #ifdef G4VERBOSE
00099 if(fVerbose>1)
00100 {
00101 G4cout<<"G4MolecularDecay::IsApplicable(";
00102 G4cout<<aParticleType.GetParticleName()<<",";
00103 G4cout<<aParticleType.GetParticleType()<<")"<<G4endl;
00104 }
00105 #endif
00106 return(true);
00107 }
00108 else
00109 {
00110 return false;
00111 }
00112 }
00113
00114 G4double G4DNAMolecularDecay::GetMeanLifeTime(const G4Track& track ,
00115 G4ForceCondition*)
00116 {
00117 G4double output = GetMolecule(track)-> GetDecayTime() - track.GetProperTime() ;
00118 return (output > 0 ? output : 0 );
00119 }
00120
00121 G4VParticleChange* G4DNAMolecularDecay::DecayIt(
00122 const G4Track& track,
00123 const G4Step& )
00124 {
00125
00126
00127
00128 aParticleChange.Initialize(track);
00129 const G4Molecule * theMotherMolecule = GetMolecule(track);
00130 const G4MoleculeDefinition* moleculeDefinition = theMotherMolecule->GetDefinition();
00131
00132
00133
00134
00135
00136
00137 if(moleculeDefinition-> GetDecayTable())
00138 {
00139 const vector<const G4MolecularDecayChannel*>* DecayVector =
00140 (theMotherMolecule -> GetDecayChannel());
00141
00142 if(DecayVector == 0)
00143 {
00144 G4ExceptionDescription exceptionDescription;
00145 theMotherMolecule->GetElectronOccupancy()->DumpInfo();
00146 exceptionDescription << "No decay channel was found for the molecule : " << theMotherMolecule-> GetName() << G4endl;
00147 G4Exception("G4DNAMolecularDecay::DecayIt", "G4DNAMolecularDecay::NoDecayChannel",FatalException,exceptionDescription);
00148 return &aParticleChange;
00149 }
00150
00151 G4int DecayVectorSize = DecayVector-> size();
00152
00153
00154 G4double RdmValue = G4UniformRand();
00155
00156 const G4MolecularDecayChannel* decayChannel(0);
00157 G4int i=0;
00158 do
00159 {
00160 decayChannel = (*DecayVector)[i];
00161 if(RdmValue < decayChannel->GetProbability()) break;
00162 RdmValue -= decayChannel->GetProbability();
00163 i++;
00164 }
00165 while(i< DecayVectorSize);
00166
00167
00168
00169
00170 G4double decayEnergy = decayChannel->GetEnergy();
00171 G4int nbProducts = decayChannel->GetNbProducts();
00172
00173 if(decayEnergy)
00174 {
00175
00176
00177
00178 aParticleChange.ProposeLocalEnergyDeposit(decayChannel->GetEnergy());
00179 }
00180
00181 if(nbProducts)
00182 {
00183
00184
00185
00186
00187 vector<G4ThreeVector> ProductsDisplacement(nbProducts);
00188 G4ThreeVector theMotherMoleculeDisplacement;
00189
00190 DecayDisplacementMap::iterator it = fDecayDisplacementMap.find(moleculeDefinition);
00191
00192 if(it!=fDecayDisplacementMap.end())
00193 {
00194 G4VMolecularDecayDisplacer* displacer = it->second;
00195 ProductsDisplacement = displacer->GetProductsDisplacement(decayChannel);
00196 theMotherMoleculeDisplacement = displacer-> GetMotherMoleculeDisplacement(decayChannel);
00197 }
00198 else
00199 {
00200 G4ExceptionDescription errMsg;
00201 errMsg << "No G4MolecularDecayProcess::theDecayDisplacementMap["
00202 << theMotherMolecule->GetName() +"]" ;
00203 G4Exception("G4MolecularDecayProcess::DecayIt","DNAMolecularDecay001",FatalErrorInArgument, errMsg);
00204 }
00205
00206 aParticleChange.SetNumberOfSecondaries(nbProducts);
00207
00208 #ifdef G4VERBOSE
00209 if(fVerbose)
00210 {
00211 G4cout<<"Decay Process : "
00212 << theMotherMolecule->GetName()
00213 << " (trackID :" << track.GetTrackID() << ") "
00214 << decayChannel->GetName()
00215 << G4endl;
00216 }
00217 #endif
00218
00219 for (G4int j=0; j<nbProducts ; j++)
00220 {
00221 G4Molecule* product = new G4Molecule(*decayChannel->GetProduct(j));
00222
00223
00224
00225
00226
00227 G4Track* secondary = product->BuildTrack(picosecond,track.GetPosition()
00228 + theMotherMoleculeDisplacement + ProductsDisplacement[j]);
00229
00230 secondary-> SetTrackStatus(fAlive);
00231 #ifdef G4VERBOSE
00232 if(fVerbose)
00233 {
00234 G4cout<<"Product : "<< product->GetName()<<G4endl;
00235 }
00236 #endif
00237
00238 aParticleChange.G4VParticleChange::AddSecondary(secondary);
00239 G4ITManager<G4Molecule>::Instance()->Push(secondary);
00240 }
00241 #ifdef G4VERBOSE
00242 if(fVerbose)
00243 G4cout<<"-------------"<<G4endl;
00244 #endif
00245 }
00246
00247
00248
00249
00250
00251
00252 else if(!decayEnergy && !nbProducts)
00253 {
00254 G4ExceptionDescription errMsg;
00255 errMsg << "There is no products and no energy specified in the molecular decay channel";
00256 G4Exception("G4MolecularDecayProcess::DecayIt","DNAMolecularDecay002",FatalErrorInArgument, errMsg);
00257 }
00258 }
00259
00260 aParticleChange.ProposeTrackStatus(fStopAndKill);
00261
00262 return &aParticleChange;
00263 }
00264
00265 void G4DNAMolecularDecay::SetDecayDisplacer (const G4ParticleDefinition* molDef, G4VMolecularDecayDisplacer* aDisplacer)
00266 {
00267 fDecayDisplacementMap[molDef] = aDisplacer;
00268 }
00269
00270 G4VMolecularDecayDisplacer* G4DNAMolecularDecay::GetDecayDisplacer(const G4ParticleDefinition* molDef)
00271 {
00272 return fDecayDisplacementMap[molDef] ;
00273 }