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