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