Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * 20 // * * 21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by QinetiQ Ltd * 22 // * under contract to the European Space Agen 22 // * under contract to the European Space Agency (ESA) are the * 23 // * intellectual property of ESA. Rights to u 23 // * intellectual property of ESA. Rights to use, copy, modify and * 24 // * redistribute this software for general pu 24 // * redistribute this software for general public use are granted * 25 // * in compliance with any licensing, distrib 25 // * in compliance with any licensing, distribution and development * 26 // * policy adopted by the Geant4 Collaboratio 26 // * policy adopted by the Geant4 Collaboration. This code has been * 27 // * written by QinetiQ Ltd for the European S 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA * 28 // * contract 17191/03/NL/LvH (Aurora Programm 28 // * contract 17191/03/NL/LvH (Aurora Programme). * 29 // * 29 // * * 30 // * By using, copying, modifying or distri 30 // * By using, copying, modifying or distributing the software (or * 31 // * any work based on the software) you ag 31 // * any work based on the software) you agree to acknowledge its * 32 // * use in resulting scientific publicati 32 // * use in resulting scientific publications, and indicate your * 33 // * acceptance of all terms of the Geant4 Sof 33 // * acceptance of all terms of the Geant4 Software license. * 34 // ******************************************* 34 // ******************************************************************** 35 // 35 // 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37 // 37 // 38 // MODULE: G4EMDissociation.cc 38 // MODULE: G4EMDissociation.cc 39 // 39 // 40 // Version: B.1 40 // Version: B.1 41 // Date: 15/04/04 41 // Date: 15/04/04 42 // Author: P R Truscott 42 // Author: P R Truscott 43 // Organisation: QinetiQ Ltd, UK 43 // Organisation: QinetiQ Ltd, UK 44 // Customer: ESA/ESTEC, NOORDWIJK 44 // Customer: ESA/ESTEC, NOORDWIJK 45 // Contract: 17191/03/NL/LvH 45 // Contract: 17191/03/NL/LvH 46 // 46 // 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 // 48 // 49 // CHANGE HISTORY 49 // CHANGE HISTORY 50 // -------------- 50 // -------------- 51 // 51 // 52 // 17 October 2003, P R Truscott, QinetiQ Ltd, 52 // 17 October 2003, P R Truscott, QinetiQ Ltd, UK 53 // Created. 53 // Created. 54 // 54 // 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK 56 // Beta release 56 // Beta release 57 // 57 // 58 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 58 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 59 ////////////////////////////////////////////// 59 //////////////////////////////////////////////////////////////////////////////// 60 // 60 // 61 #include "G4EMDissociation.hh" 61 #include "G4EMDissociation.hh" 62 #include "G4PhysicalConstants.hh" << 62 #include "G4Evaporation.hh" 63 #include "G4SystemOfUnits.hh" << 63 #include "G4FermiBreakUp.hh" >> 64 #include "G4StatMF.hh" 64 #include "G4ParticleDefinition.hh" 65 #include "G4ParticleDefinition.hh" 65 #include "G4LorentzVector.hh" 66 #include "G4LorentzVector.hh" 66 #include "G4PhysicsFreeVector.hh" 67 #include "G4PhysicsFreeVector.hh" 67 #include "G4EMDissociationCrossSection.hh" 68 #include "G4EMDissociationCrossSection.hh" 68 #include "G4Proton.hh" 69 #include "G4Proton.hh" 69 #include "G4Neutron.hh" 70 #include "G4Neutron.hh" >> 71 #include "G4ParticleTable.hh" 70 #include "G4IonTable.hh" 72 #include "G4IonTable.hh" >> 73 #include "G4GeneralPhaseSpaceDecay.hh" 71 #include "G4DecayProducts.hh" 74 #include "G4DecayProducts.hh" 72 #include "G4DynamicParticle.hh" 75 #include "G4DynamicParticle.hh" 73 #include "G4Fragment.hh" 76 #include "G4Fragment.hh" 74 #include "G4ReactionProductVector.hh" 77 #include "G4ReactionProductVector.hh" 75 #include "Randomize.hh" 78 #include "Randomize.hh" 76 #include "globals.hh" 79 #include "globals.hh" 77 #include "G4PhysicsModelCatalog.hh" << 80 //////////////////////////////////////////////////////////////////////////////// 78 << 81 // 79 G4EMDissociation::G4EMDissociation() : << 82 G4EMDissociation::G4EMDissociation():G4HadronicInteraction("EMDissociation") 80 G4HadronicInteraction("EMDissociation"), << 83 { 81 secID_projectileDissociation(-1), secID_targ << 84 // 82 { << 85 // 83 // Send message to stdout to advise that the << 86 // Send message to stdout to advise that the G4EMDissociation model is being 84 // used. << 87 // used. >> 88 // 85 PrintWelcomeMessage(); 89 PrintWelcomeMessage(); 86 << 90 // 87 // No de-excitation handler has been supplie << 91 // 88 theExcitationHandler = new G4Exci << 92 // No de-excitation handler has been supplied - define the default handler. >> 93 // >> 94 theExcitationHandler = new G4ExcitationHandler; >> 95 G4Evaporation * theEvaporation = new G4Evaporation; >> 96 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; >> 97 G4StatMF * theMF = new G4StatMF; >> 98 theExcitationHandler->SetEvaporation(theEvaporation); >> 99 theExcitationHandler->SetFermiModel(theFermiBreakUp); >> 100 theExcitationHandler->SetMultiFragmentation(theMF); >> 101 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); 89 theExcitationHandler->SetMinEForMultiFrag(5. 102 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); 90 handlerDefinedInternally = true; 103 handlerDefinedInternally = true; 91 << 104 // 92 // This EM dissociation model needs access t << 105 // 93 // G4EMDissociationCrossSection. << 106 // This EM dissociation model needs access to the cross-sections held in >> 107 // G4EMDissociationCrossSection. >> 108 // 94 dissociationCrossSection = new G4EMDissociat 109 dissociationCrossSection = new G4EMDissociationCrossSection; 95 thePhotonSpectrum = new G4EMDissociationSpec 110 thePhotonSpectrum = new G4EMDissociationSpectrum; 96 << 111 // 97 // Set the minimum and maximum range for the << 112 // 98 // is in energy per nucleon number). << 113 // Set the minimum and maximum range for the model (despite nomanclature, this >> 114 // is in energy per nucleon number). >> 115 // 99 SetMinEnergy(100.0*MeV); 116 SetMinEnergy(100.0*MeV); 100 SetMaxEnergy(500.0*GeV); 117 SetMaxEnergy(500.0*GeV); 101 << 118 // 102 // Set the default verbose level to 0 - no o << 119 // >> 120 // Set the default verbose level to 0 - no output. >> 121 // 103 verboseLevel = 0; 122 verboseLevel = 0; 104 << 105 // Creator model ID for the secondaries crea << 106 secID_projectileDissociation = G4PhysicsMode << 107 secID_targetDissociation = G4PhysicsMode << 108 } 123 } 109 << 124 //////////////////////////////////////////////////////////////////////////////// 110 G4EMDissociation::G4EMDissociation (G4Excitati << 125 // 111 G4HadronicInteraction("EMDissociation"), << 126 G4EMDissociation::G4EMDissociation (G4ExcitationHandler *aExcitationHandler) 112 secID_projectileDissociation(-1), secID_targ << 113 { 127 { 114 // Send message to stdout to advise that the << 128 // 115 // used. << 129 // >> 130 // Send message to stdout to advise that the G4EMDissociation model is being >> 131 // used. >> 132 // 116 PrintWelcomeMessage(); 133 PrintWelcomeMessage(); 117 134 118 theExcitationHandler = aExcitationHandle 135 theExcitationHandler = aExcitationHandler; 119 handlerDefinedInternally = false; 136 handlerDefinedInternally = false; 120 << 137 // 121 // This EM dissociation model needs access t << 138 // 122 // G4EMDissociationCrossSection. << 139 // This EM dissociation model needs access to the cross-sections held in >> 140 // G4EMDissociationCrossSection. >> 141 // 123 dissociationCrossSection = new G4EMDissociat 142 dissociationCrossSection = new G4EMDissociationCrossSection; 124 thePhotonSpectrum = new G4EMDissociationSpec 143 thePhotonSpectrum = new G4EMDissociationSpectrum; 125 << 144 // 126 // Set the minimum and maximum range for the << 145 // 127 // is in energy per nucleon number) << 146 // Set the minimum and maximum range for the model (despite nomanclature, this >> 147 // is in energy per nucleon number). >> 148 // 128 SetMinEnergy(100.0*MeV); 149 SetMinEnergy(100.0*MeV); 129 SetMaxEnergy(500.0*GeV); 150 SetMaxEnergy(500.0*GeV); >> 151 // >> 152 // >> 153 // Set the default verbose level to 0 - no output. >> 154 // 130 verboseLevel = 0; 155 verboseLevel = 0; 131 << 132 // Creator model ID for the secondaries crea << 133 secID_projectileDissociation = G4PhysicsMode << 134 secID_targetDissociation = G4PhysicsMode << 135 } 156 } 136 << 157 //////////////////////////////////////////////////////////////////////////////// 137 << 158 // 138 G4EMDissociation::~G4EMDissociation() { << 159 G4EMDissociation::~G4EMDissociation () >> 160 { 139 if (handlerDefinedInternally) delete theExci 161 if (handlerDefinedInternally) delete theExcitationHandler; 140 // delete dissociationCrossSection; << 162 delete dissociationCrossSection; 141 // Cross section deleted by G4CrossSectionRe << 142 // Bug reported by Gong Ding in Bug Report # << 143 delete thePhotonSpectrum; 163 delete thePhotonSpectrum; 144 } 164 } 145 << 165 //////////////////////////////////////////////////////////////////////////////// 146 << 166 // 147 G4HadFinalState *G4EMDissociation::ApplyYourse 167 G4HadFinalState *G4EMDissociation::ApplyYourself 148 (const G4HadProjectile &theTrack, G4Nucleus 168 (const G4HadProjectile &theTrack, G4Nucleus &theTarget) 149 { 169 { 150 // The secondaries will be returned in G4Had << 170 // 151 // initialise this. << 171 // 152 << 172 // The secondaries will be returned in G4HadFinalState &theParticleChange - >> 173 // initialise this. >> 174 // 153 theParticleChange.Clear(); 175 theParticleChange.Clear(); 154 theParticleChange.SetStatusChange(stopAndKil 176 theParticleChange.SetStatusChange(stopAndKill); 155 << 177 // 156 // Get relevant information about the projec << 178 // 157 // energy/nuc, momentum, velocity, Lorentz f << 179 // Get relevant information about the projectile and target (A, Z) and 158 // projectile. << 180 // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the 159 << 181 // projectile. >> 182 // 160 const G4ParticleDefinition *definitionP = th 183 const G4ParticleDefinition *definitionP = theTrack.GetDefinition(); 161 const G4double AP = definitionP->GetBaryonN 184 const G4double AP = definitionP->GetBaryonNumber(); 162 const G4double ZP = definitionP->GetPDGChar 185 const G4double ZP = definitionP->GetPDGCharge(); 163 G4LorentzVector pP = theTrack.Get4Momentum() 186 G4LorentzVector pP = theTrack.Get4Momentum(); 164 G4double E = theTrack.GetKineticEner 187 G4double E = theTrack.GetKineticEnergy()/AP; 165 G4double MP = theTrack.GetTotalEnergy 188 G4double MP = theTrack.GetTotalEnergy() - E*AP; 166 G4double b = pP.beta(); 189 G4double b = pP.beta(); 167 G4double AT = theTarget.GetA_asInt(); 190 G4double AT = theTarget.GetA_asInt(); 168 G4double ZT = theTarget.GetZ_asInt(); 191 G4double ZT = theTarget.GetZ_asInt(); 169 G4double MT = G4NucleiProperties::Get 192 G4double MT = G4NucleiProperties::GetNuclearMass(AT,ZT); 170 << 193 // 171 // Depending upon the verbosity level, outpu << 194 // 172 // projectile and target << 195 // Depending upon the verbosity level, output the initial information on the 173 if (verboseLevel >= 2) { << 196 // projectile and target. >> 197 // >> 198 if (verboseLevel >= 2) >> 199 { 174 G4cout.precision(6); 200 G4cout.precision(6); 175 G4cout <<"################################ 201 G4cout <<"########################################" 176 <<"################################ 202 <<"########################################" 177 <<G4endl; 203 <<G4endl; 178 G4cout <<"IN G4EMDissociation" <<G4endl; 204 G4cout <<"IN G4EMDissociation" <<G4endl; 179 G4cout <<"Initial projectile A=" <<AP 205 G4cout <<"Initial projectile A=" <<AP 180 <<", Z=" <<ZP 206 <<", Z=" <<ZP 181 <<G4endl; 207 <<G4endl; 182 G4cout <<"Initial target A=" <<AT 208 G4cout <<"Initial target A=" <<AT 183 <<", Z=" <<ZT 209 <<", Z=" <<ZT 184 <<G4endl; 210 <<G4endl; 185 G4cout <<"Projectile momentum and Energy/n 211 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl; 186 } 212 } 187 << 213 // 188 // Initialise the variables which will be us << 214 // 189 // to boost the secondaries from the interac << 215 // Initialise the variables which will be used with the phase-space decay and 190 << 216 // to boost the secondaries from the interaction. >> 217 // 191 G4ParticleDefinition *typeNucleon = NULL; 218 G4ParticleDefinition *typeNucleon = NULL; 192 G4ParticleDefinition *typeDaughter = NULL; 219 G4ParticleDefinition *typeDaughter = NULL; 193 G4double Eg = 0.0; 220 G4double Eg = 0.0; 194 G4double mass = 0.0; 221 G4double mass = 0.0; 195 G4ThreeVector boost = G4ThreeVector(0.0, 0.0 222 G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0); 196 << 223 // 197 // Determine the cross-sections at the giant << 224 // 198 // resonance energies for the projectile and << 225 // Determine the cross-sections at the giant dipole and giant quadrupole 199 // initially provided in the G4PhysicsFreeVe << 226 // resonance energies for the projectile and then target. The information is 200 // and E2 fields. These are then summed. << 227 // initially provided in the G4PhysicsFreeVector individually for the E1 201 << 228 // and E2 fields. These are then summed. >> 229 // 202 G4double bmin = thePhotonSpectrum->GetCloses 230 G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b); 203 G4PhysicsFreeVector *crossSectionP = dissoci 231 G4PhysicsFreeVector *crossSectionP = dissociationCrossSection-> 204 GetCrossSectionForProjectile(AP, ZP, AT, Z 232 GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin); 205 G4PhysicsFreeVector *crossSectionT = dissoci 233 G4PhysicsFreeVector *crossSectionT = dissociationCrossSection-> 206 GetCrossSectionForTarget(AP, ZP, AT, ZT, b 234 GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin); 207 235 208 G4double totCrossSectionP = (*crossSectionP) 236 G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1]; 209 G4double totCrossSectionT = (*crossSectionT) 237 G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1]; 210 << 238 // 211 // Now sample whether the interaction involv << 239 // 212 // or the target. << 240 // Now sample whether the interaction involved EM dissociation of the projectile 213 << 241 // or the target. 214 G4int secID = -1; // Creator model ID for t << 242 // 215 if (G4UniformRand() < 243 if (G4UniformRand() < 216 totCrossSectionP / (totCrossSectionP + tot << 244 totCrossSectionP / (totCrossSectionP + totCrossSectionT)) 217 << 245 { 218 // It was the projectile which underwent E << 246 // 219 // boost to be applied to the secondaries, << 247 // 220 // neutron was ejected. Then determine th << 248 // It was the projectile which underwent EM dissociation. Define the Lorentz 221 // which passed from the target nucleus .. << 249 // boost to be applied to the secondaries, and sample whether a proton or a 222 // excitation of the projectile. << 250 // neutron was ejected. Then determine the energy of the virtual gamma ray 223 << 251 // which passed from the target nucleus ... this will be used to define the 224 secID = secID_projectileDissociation; << 252 // excitation of the projectile. >> 253 // 225 mass = MP; 254 mass = MP; 226 if (G4UniformRand() < dissociationCrossSec 255 if (G4UniformRand() < dissociationCrossSection-> 227 GetWilsonProbabilityForProtonDissociatio 256 GetWilsonProbabilityForProtonDissociation (AP, ZP)) 228 { 257 { 229 if (verboseLevel >= 2) 258 if (verboseLevel >= 2) 230 G4cout <<"Projectile underwent EM diss 259 G4cout <<"Projectile underwent EM dissociation producing a proton" 231 <<G4endl; 260 <<G4endl; 232 typeNucleon = G4Proton::ProtonDefinition 261 typeNucleon = G4Proton::ProtonDefinition(); 233 typeDaughter = G4IonTable::GetIonTable() << 262 typeDaughter = G4ParticleTable::GetParticleTable()-> 234 GetIon((G4int) ZP-1, (G4int) AP-1, 0.0); 263 GetIon((G4int) ZP-1, (G4int) AP-1, 0.0); 235 } 264 } 236 else 265 else 237 { 266 { 238 if (verboseLevel >= 2) 267 if (verboseLevel >= 2) 239 G4cout <<"Projectile underwent EM diss 268 G4cout <<"Projectile underwent EM dissociation producing a neutron" 240 <<G4endl; 269 <<G4endl; 241 typeNucleon = G4Neutron::NeutronDefiniti 270 typeNucleon = G4Neutron::NeutronDefinition(); 242 typeDaughter = G4IonTable::GetIonTable() << 271 typeDaughter = G4ParticleTable::GetParticleTable()-> 243 GetIon((G4int) ZP, (G4int) AP-1, 0.0); 272 GetIon((G4int) ZP, (G4int) AP-1, 0.0); 244 } 273 } 245 if (G4UniformRand() < (*crossSectionP)[0]/ 274 if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP) 246 { 275 { 247 Eg = crossSectionP->GetLowEdgeEnergy(0); 276 Eg = crossSectionP->GetLowEdgeEnergy(0); 248 if (verboseLevel >= 2) 277 if (verboseLevel >= 2) 249 G4cout <<"Transition type was E1" <<G4 278 G4cout <<"Transition type was E1" <<G4endl; 250 } 279 } 251 else 280 else 252 { 281 { 253 Eg = crossSectionP->GetLowEdgeEnergy(1); 282 Eg = crossSectionP->GetLowEdgeEnergy(1); 254 if (verboseLevel >= 2) 283 if (verboseLevel >= 2) 255 G4cout <<"Transition type was E2" <<G4 284 G4cout <<"Transition type was E2" <<G4endl; 256 } 285 } 257 << 286 // 258 // We need to define a Lorentz vector with << 287 // 259 // energy includes the projectile and virt << 288 // We need to define a Lorentz vector with the original momentum, but total 260 // to calculate the boost required for the << 289 // energy includes the projectile and virtual gamma. This is then used 261 << 290 // to calculate the boost required for the secondaries. 262 pP.setE( std::sqrt( pP.vect().mag2() + (ma << 291 // >> 292 pP.setE(pP.e()+Eg); 263 boost = pP.findBoostToCM(); 293 boost = pP.findBoostToCM(); 264 } 294 } 265 else 295 else 266 { 296 { 267 // It was the target which underwent EM di << 297 // 268 // proton or a neutron was ejected. Then << 298 // 269 // gamma ray which passed from the project << 299 // It was the target which underwent EM dissociation. Sample whether a 270 // define the excitation of the target. << 300 // proton or a neutron was ejected. Then determine the energy of the virtual 271 << 301 // gamma ray which passed from the projectile nucleus ... this will be used to 272 secID = secID_targetDissociation; << 302 // define the excitation of the target. >> 303 // 273 mass = MT; 304 mass = MT; 274 if (G4UniformRand() < dissociationCrossSec 305 if (G4UniformRand() < dissociationCrossSection-> 275 GetWilsonProbabilityForProtonDissociatio 306 GetWilsonProbabilityForProtonDissociation (AT, ZT)) 276 { 307 { 277 if (verboseLevel >= 2) 308 if (verboseLevel >= 2) 278 G4cout <<"Target underwent EM dissocia 309 G4cout <<"Target underwent EM dissociation producing a proton" 279 <<G4endl; 310 <<G4endl; 280 typeNucleon = G4Proton::ProtonDefinition 311 typeNucleon = G4Proton::ProtonDefinition(); 281 typeDaughter = G4IonTable::GetIonTable() << 312 typeDaughter = G4ParticleTable::GetParticleTable()-> 282 GetIon((G4int) ZT-1, (G4int) AT-1, 0.0); 313 GetIon((G4int) ZT-1, (G4int) AT-1, 0.0); 283 } 314 } 284 else 315 else 285 { 316 { 286 if (verboseLevel >= 2) 317 if (verboseLevel >= 2) 287 G4cout <<"Target underwent EM dissocia 318 G4cout <<"Target underwent EM dissociation producing a neutron" 288 <<G4endl; 319 <<G4endl; 289 typeNucleon = G4Neutron::NeutronDefiniti 320 typeNucleon = G4Neutron::NeutronDefinition(); 290 typeDaughter = G4IonTable::GetIonTable() << 321 typeDaughter = G4ParticleTable::GetParticleTable()-> 291 GetIon((G4int) ZT, (G4int) AT-1, 0.0); 322 GetIon((G4int) ZT, (G4int) AT-1, 0.0); 292 } 323 } 293 if (G4UniformRand() < (*crossSectionT)[0]/ 324 if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT) 294 { 325 { 295 Eg = crossSectionT->GetLowEdgeEnergy(0); 326 Eg = crossSectionT->GetLowEdgeEnergy(0); 296 if (verboseLevel >= 2) 327 if (verboseLevel >= 2) 297 G4cout <<"Transition type was E1" <<G4 328 G4cout <<"Transition type was E1" <<G4endl; 298 } 329 } 299 else 330 else 300 { 331 { 301 Eg = crossSectionT->GetLowEdgeEnergy(1); 332 Eg = crossSectionT->GetLowEdgeEnergy(1); 302 if (verboseLevel >= 2) 333 if (verboseLevel >= 2) 303 G4cout <<"Transition type was E2" <<G4 334 G4cout <<"Transition type was E2" <<G4endl; 304 } 335 } 305 << 336 // 306 // Add the projectile to theParticleChange << 337 // 307 // not-so-virtual gamma-ray. Not that at << 338 // Add the projectile to theParticleChange, less the energy of the 308 // is transferred between the projectile a << 339 // not-so-virtual gamma-ray. Not that at the moment, no lateral momentum 309 << 340 // is transferred between the projectile and target nuclei. >> 341 // 310 G4ThreeVector v = pP.vect(); 342 G4ThreeVector v = pP.vect(); 311 v.setMag(1.0); 343 v.setMag(1.0); 312 G4DynamicParticle *changedP = new G4Dynami << 344 G4DynamicParticle *changedP = new G4DynamicParticle 313 theParticleChange.AddSecondary (changedP, << 345 (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg); >> 346 theParticleChange.AddSecondary (changedP); 314 if (verboseLevel >= 2) 347 if (verboseLevel >= 2) 315 { 348 { 316 G4cout <<"Projectile change:" <<G4endl; 349 G4cout <<"Projectile change:" <<G4endl; 317 changedP->DumpInfo(); 350 changedP->DumpInfo(); 318 } 351 } 319 } 352 } 320 << 353 // 321 // Perform a two-body decay based on the res << 354 // 322 // gamma-ray, and the masses of the daughter << 355 // Perform a two-body decay based on the restmass energy of the parent and 323 // the nucles, the angular distribution is s << 356 // gamma-ray, and the masses of the daughters. In the frame of reference of 324 // the nucleon and secondary nucleus are boo << 357 // the nucles, the angular distribution is sampled isotropically, but the 325 // projectile. << 358 // the nucleon and secondary nucleus are boosted if they've come from the 326 << 359 // projectile. >> 360 // 327 G4double e = mass + Eg; 361 G4double e = mass + Eg; 328 G4double mass1 = typeNucleon->GetPDGMass(); << 362 G4double m1 = typeNucleon->GetPDGMass(); 329 G4double mass2 = typeDaughter->GetPDGMass(); << 363 G4double m2 = typeDaughter->GetPDGMass(); 330 G4double pp = (e+mass1+mass2)*(e+mass1-mass2 << 364 G4double pp = (e+m1+m2)*(e+m1-m2)*(e-m1+m2)*(e-m1-m2)/(4.0*e*e); 331 (e-mass1+mass2)*(e-mass1-mass2 << 365 if (pp < 0.0) 332 if (pp < 0.0) { << 366 { 333 pp = 1.0*eV; 367 pp = 1.0*eV; 334 // if (verboseLevel >`= 1) 368 // if (verboseLevel >`= 1) 335 // { 369 // { 336 // G4cout <<"IN G4EMDissociation::ApplyYo 370 // G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl; 337 // G4cout <<"Error in mass of secondaries 371 // G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl; 338 // G4cout <<"Rest mass of primary = 372 // G4cout <<"Rest mass of primary = " <<mass <<" MeV" <<G4endl; 339 // G4cout <<"Virtual gamma energy = 373 // G4cout <<"Virtual gamma energy = " <<Eg <<" MeV" <<G4endl; 340 // G4cout <<"Rest mass of secondary #1 = << 374 // G4cout <<"Rest mass of secondary #1 = " <<m1 <<" MeV" <<G4endl; 341 // G4cout <<"Rest mass of secondary #2 = << 375 // G4cout <<"Rest mass of secondary #2 = " <<m2 <<" MeV" <<G4endl; 342 // } 376 // } 343 } 377 } 344 else 378 else 345 pp = std::sqrt(pp); 379 pp = std::sqrt(pp); 346 G4double costheta = 2.*G4UniformRand()-1.0; 380 G4double costheta = 2.*G4UniformRand()-1.0; 347 G4double sintheta = std::sqrt((1.0 - costhet 381 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); 348 G4double phi = 2.0*pi*G4UniformRand()*r 382 G4double phi = 2.0*pi*G4UniformRand()*rad; 349 G4ThreeVector direction(sintheta*std::cos(ph 383 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); 350 G4DynamicParticle *dynamicNucleon = 384 G4DynamicParticle *dynamicNucleon = 351 new G4DynamicParticle(typeNucleon, directi 385 new G4DynamicParticle(typeNucleon, direction*pp); 352 dynamicNucleon->Set4Momentum(dynamicNucleon- 386 dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost)); 353 G4DynamicParticle *dynamicDaughter = 387 G4DynamicParticle *dynamicDaughter = 354 new G4DynamicParticle(typeDaughter, -direc 388 new G4DynamicParticle(typeDaughter, -direction*pp); 355 dynamicDaughter->Set4Momentum(dynamicDaughte 389 dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost)); 356 << 390 // 357 // The "decay" products have to be transferr << 391 // 358 // Furthermore, the residual nucleus should << 392 // The "decay" products have to be transferred to the G4HadFinalState object. 359 << 393 // Furthermore, the residual nucleus should be de-excited. 360 theParticleChange.AddSecondary (dynamicNucle << 394 // 361 if (verboseLevel >= 2) { << 395 theParticleChange.AddSecondary (dynamicNucleon); >> 396 if (verboseLevel >= 2) >> 397 { 362 G4cout <<"Nucleon from the EMD process:" < 398 G4cout <<"Nucleon from the EMD process:" <<G4endl; 363 dynamicNucleon->DumpInfo(); 399 dynamicNucleon->DumpInfo(); 364 } 400 } 365 401 366 G4Fragment* theFragment = new << 402 G4Fragment *theFragment = new 367 G4Fragment(typeDaughter->GetBaryonNumber() << 403 G4Fragment((G4int) typeDaughter->GetBaryonNumber(), 368 G4lrint(typeDaughter->GetPDGCharge()/ << 404 (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum()); 369 dynamicDaughter->Get4Momentum()); << 405 if (verboseLevel >= 2) 370 << 406 { 371 if (verboseLevel >= 2) { << 372 G4cout <<"Dynamic properties of the prefra 407 G4cout <<"Dynamic properties of the prefragment:" <<G4endl; 373 G4cout.precision(6); 408 G4cout.precision(6); 374 dynamicDaughter->DumpInfo(); 409 dynamicDaughter->DumpInfo(); 375 G4cout <<"Nuclear properties of the prefra 410 G4cout <<"Nuclear properties of the prefragment:" <<G4endl; 376 G4cout <<theFragment <<G4endl; 411 G4cout <<theFragment <<G4endl; 377 } 412 } 378 << 413 G4ReactionProductVector *products = 379 G4ReactionProductVector* products = << 414 theExcitationHandler->BreakItUp(*theFragment); 380 theExcitationHandler->Br << 381 delete theFragment; 415 delete theFragment; 382 theFragment = NULL; 416 theFragment = NULL; 383 417 384 G4DynamicParticle* secondary = 0; << 385 G4ReactionProductVector::iterator iter; 418 G4ReactionProductVector::iterator iter; 386 for (iter = products->begin(); iter != produ << 419 for (iter = products->begin(); iter != products->end(); ++iter) 387 secondary = new G4DynamicParticle((*iter)- << 420 { 388 (*iter)->GetTotalEnergy(), (*iter)->GetMom << 421 G4DynamicParticle *secondary = 389 theParticleChange.AddSecondary (secondary, << 422 new G4DynamicParticle((*iter)->GetDefinition(), >> 423 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); >> 424 theParticleChange.AddSecondary (secondary); 390 } 425 } 391 delete products; << 392 << 393 delete crossSectionP; << 394 delete crossSectionT; << 395 426 396 if (verboseLevel >= 2) 427 if (verboseLevel >= 2) 397 G4cout <<"################################ 428 G4cout <<"########################################" 398 <<"################################ 429 <<"########################################" 399 <<G4endl; 430 <<G4endl; 400 431 401 return &theParticleChange; 432 return &theParticleChange; 402 } 433 } 403 << 434 //////////////////////////////////////////////////////////////////////////////// 404 << 435 // 405 void G4EMDissociation::PrintWelcomeMessage () 436 void G4EMDissociation::PrintWelcomeMessage () 406 { 437 { 407 G4cout <<G4endl; 438 G4cout <<G4endl; 408 G4cout <<" ********************************* 439 G4cout <<" ****************************************************************" 409 <<G4endl; 440 <<G4endl; 410 G4cout <<" EM dissociation model for nuclear 441 G4cout <<" EM dissociation model for nuclear-nuclear interactions activated" 411 <<G4endl; 442 <<G4endl; 412 G4cout <<" (Written by QinetiQ Ltd for the E 443 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)" 413 <<G4endl; 444 <<G4endl; 414 G4cout <<" ********************************* 445 G4cout <<" ****************************************************************" 415 <<G4endl; 446 <<G4endl; 416 G4cout << G4endl; 447 G4cout << G4endl; 417 448 418 return; 449 return; 419 } 450 } 420 << 451 //////////////////////////////////////////////////////////////////////////////// >> 452 // 421 453