Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4IntraNucleiCascader.cc @ 1007

Last change on this file since 1007 was 1007, checked in by garnier, 15 years ago
update to geant4.9.2
File size: 11.9 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26
27	#define RUN
28
29	#include "G4IntraNucleiCascader.hh"
30	#include "G4InuclElementaryParticle.hh"
31	#include "G4InuclNuclei.hh"
32	#include "G4LorentzConvertor.hh"
33	#include "G4ParticleLargerEkin.hh"
34	#include "G4NucleiModel.hh"
35	#include "G4CascadParticle.hh"
36	#include "Randomize.hh"
37	#include <algorithm>
38
39	typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
40
41	G4IntraNucleiCascader::G4IntraNucleiCascader()
42	: verboseLevel(1) {
43
44	if (verboseLevel > 3) {
45	G4cout << " >>> G4IntraNucleiCascader::G4IntraNucleiCascader" << G4endl;
46	}
47
48	}
49
50	G4CollisionOutput G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
51	G4InuclParticle* target) {
52
53	if (verboseLevel > 3) {
54	G4cout << " >>> G4IntraNucleiCascader::collide" << G4endl;
55	}
56
57	const G4int itry_max = 1000;
58	const G4int reflection_cut = 500;
59	// const G4double eexs_cut = 0.0001;
60
61	if (verboseLevel > 3) {
62	bullet->printParticle();
63	target->printParticle();
64	}
65
66	G4CollisionOutput output;
67
68	#ifdef RUN
69	G4InuclNuclei* tnuclei = dynamic_cast<G4InuclNuclei*>(target);
70	G4InuclNuclei* bnuclei = dynamic_cast<G4InuclNuclei*>(bullet);
71	G4InuclElementaryParticle* bparticle =
72	dynamic_cast<G4InuclElementaryParticle*>(bullet);
73	G4NucleiModel model(tnuclei);
74	G4double coulombBarrier = 0.00126*tnuclei->getZ()/
75	(1.+std::pow(tnuclei->getA(),0.333));
76
77	G4CascadeMomentum momentum_in = bullet->getMomentum();
78
79	momentum_in[0] += tnuclei->getMass();
80
81	G4double ekin_in;
82
83	if (verboseLevel > 3) {
84	model.printModel();
85	G4cout << " intitial momentum E " << momentum_in[0] << " Px " << momentum_in[1]
86	<< " Py " << momentum_in[2] << " Pz " << momentum_in[3] << G4endl;
87	}
88
89	G4int itry = 0;
90
91	while (itry < itry_max) {
92	itry++;
93	model.reset();
94
95	std::vector<G4CascadParticle> cascad_particles;
96	G4ExitonConfiguration theExitonConfiguration;
97	std::vector<G4InuclElementaryParticle> output_particles;
98	G4double afin = tnuclei->getA();
99	G4double zfin = tnuclei->getZ();
100
101	if (inter_case == 1) { // particle with nuclei
102	ekin_in = bparticle->getKineticEnergy();
103	zfin += bparticle->getCharge();
104
105	if (bparticle->baryon()) afin += 1.0;
106
107	cascad_particles.push_back(model.initializeCascad(bparticle));
108
109	} else { // nuclei with nuclei
110
111	ekin_in = bnuclei->getKineticEnergy();
112
113	G4double ab = bnuclei->getA();
114	G4double zb = bnuclei->getZ();
115
116	afin += ab;
117	zfin += zb;
118
119	std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> >
120	all_particles = model.initializeCascad(bnuclei, tnuclei);
121
122	cascad_particles = all_particles.first;
123
124	for (G4int ip = 0; ip < G4int(all_particles.second.size()); ip++)
125	output_particles.push_back(all_particles.second[ip]);
126
127	if (cascad_particles.size() == 0) { // compound nuclei
128	G4int ia = G4int(ab + 0.5);
129	G4int iz = G4int(zb + 0.5);
130	G4int i;
131
132	for (i = 0; i < ia; i++) {
133	G4int knd = i < iz ? 1 : 2;
134	theExitonConfiguration.incrementQP(knd);
135	};
136
137	G4int ihn = G4int(2.0 * (ab - zb) * inuclRndm() + 0.5);
138	G4int ihz = G4int(2.0 * zb * inuclRndm() + 0.5);
139
140	for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
141
142	for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
143	};
144	};
145
146	std::vector<G4CascadParticle> new_cascad_particles;
147	G4int iloop = 0;
148
149	while (!cascad_particles.empty() && !model.empty()) {
150	iloop++;
151
152	if (verboseLevel > 3) {
153	G4cout << " Number of cparticles " << cascad_particles.size() << G4endl;
154	cascad_particles.back().print();
155	}
156
157	new_cascad_particles = model.generateParticleFate(cascad_particles.back(),
158	theElementaryParticleCollider);
159	if (verboseLevel > 2) {
160	G4cout << " ew particles " << new_cascad_particles.size() << G4endl;
161	}
162
163	// handle the result of a new step
164
165	if (new_cascad_particles.size() == 1) { // last particle goes without interaction
166	cascad_particles.pop_back();
167
168	if (model.stillInside(new_cascad_particles[0])) { // particle survives
169
170	if (verboseLevel > 3) {
171	G4cout << " still inside " << G4endl;
172	}
173
174	if (new_cascad_particles[0].getNumberOfReflections() < reflection_cut &&
175	model.worthToPropagate(new_cascad_particles[0])) { // it's ok
176
177	if (verboseLevel > 3) {
178	G4cout << " survives " << G4endl;
179	}
180
181	cascad_particles.push_back(new_cascad_particles[0]);
182
183	} else { // it becomes an exiton
184
185	if (verboseLevel > 3) {
186	G4cout << " becomes an exiton " << G4endl;
187	}
188
189	theExitonConfiguration.incrementQP(new_cascad_particles[0].getParticle().type());
190	};
191
192	} else { // particle about to leave nucleus - check for Coulomb barrier
193
194	if (verboseLevel > 3) {
195	G4cout << " Goes out " << G4endl;
196	new_cascad_particles[0].print();
197	}
198	G4InuclElementaryParticle currentParticle = new_cascad_particles[0].getParticle();
199	G4double KE = currentParticle.getKineticEnergy();
200	G4double mass = currentParticle.getMass();
201	G4double Q = currentParticle.getCharge();
202	if (KE < Q*coulombBarrier) {
203	// Calculate barrier penetration
204	G4double CBP = 0.0;
205	// if (KE > 0.0001) CBP = std::exp(-0.00126tnuclei->getZ()0.25*
206	// (1./KE - 1./coulombBarrier));
207	if (KE > 0.0001) CBP = std::exp(-0.01810.5tnuclei->getZ()*
208	(1./KE - 1./coulombBarrier)*
209	std::sqrt(mass*(coulombBarrier-KE)) );
210	if (G4UniformRand() < CBP) {
211	output_particles.push_back(currentParticle);
212	} else {
213	theExitonConfiguration.incrementQP(currentParticle.type());
214	}
215	} else {
216	output_particles.push_back(currentParticle);
217	}
218	}
219
220	} else { // interaction
221
222	cascad_particles.pop_back();
223
224	for (G4int i = 0; i < G4int(new_cascad_particles.size()); i++)
225	cascad_particles.push_back(new_cascad_particles[i]);
226
227	std::pair<G4int, G4int> holes = model.getTypesOfNucleonsInvolved();
228
229	theExitonConfiguration.incrementHoles(holes.first);
230
231	if (holes.second > 0) theExitonConfiguration.incrementHoles(holes.second);
232
233	};
234	};
235
236	// Cascade is finished. Check if it's OK.
237
238	if (verboseLevel > 3) {
239	G4cout << " Cascade finished " << G4endl
240	<< " output_particles " << output_particles.size() << G4endl;
241	}
242
243	G4CascadeMomentum momentum_out;
244	particleIterator ipart;
245
246	for (ipart = output_particles.begin(); ipart != output_particles.end(); ipart++) {
247	const G4CascadeMomentum& mom = ipart->getMomentum();
248
249	for (G4int j = 0; j < 4; j++) momentum_out[j] += mom[j];
250
251	zfin -= ipart->getCharge();
252	if (ipart->baryon()) afin -= 1.0;
253	};
254
255	if (verboseLevel > 3) {
256	G4cout << " afin " << afin << " zfin " << zfin << G4endl;
257	}
258
259	if (afin > 1.0) {
260	G4InuclNuclei outgoing_nuclei(afin, zfin);
261	G4double mass = outgoing_nuclei.getMass();
262	momentum_out[0] += mass;
263
264	for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
265
266	if (verboseLevel > 3) {
267	G4cout << " Eex + Ekin " << momentum_out[0] << G4endl;
268	}
269
270	if (momentum_out[0] > 0.0) { // Eex + Ekin > 0.0
271	G4double pnuc = momentum_out[1] * momentum_out[1] +
272	momentum_out[2] * momentum_out[2] +
273	momentum_out[3] * momentum_out[3];
274	G4double ekin = std::sqrt(mass * mass + pnuc) - mass;
275	G4double Eex = 1000.0 * (momentum_out[0] - ekin);
276
277	if (verboseLevel > 3) {
278	G4cout << " Eex " << Eex << G4endl;
279	}
280
281	if (goodCase(afin, zfin, Eex, ekin_in)) { // ok, exitation energy > cut
282	std::sort(output_particles.begin(), output_particles.end(), G4ParticleLargerEkin());
283	output.addOutgoingParticles(output_particles);
284	outgoing_nuclei.setMomentum(momentum_out);
285	outgoing_nuclei.setEnergy();
286	outgoing_nuclei.setExitationEnergy(Eex);
287	outgoing_nuclei.setExitonConfiguration(theExitonConfiguration);
288	output.addTargetFragment(outgoing_nuclei);
289
290	return output;
291	};
292	};
293
294	} else { // special case, when one has no nuclei after the cascad
295
296	if (afin == 1.0) { // recoiling nucleon
297
298	for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
299
300	G4InuclElementaryParticle last_particle;
301
302	if (zfin == 1.0) { // recoiling proton
303	last_particle.setType(1);
304
305	} else { // neutron
306	last_particle.setType(2);
307	};
308
309	last_particle.setMomentum(momentum_out);
310	output_particles.push_back(last_particle);
311	};
312
313	std::sort(output_particles.begin(), output_particles.end(), G4ParticleLargerEkin());
314	output.addOutgoingParticles(output_particles);
315
316	return output;
317	};
318	};
319
320	#else
321
322	// special branch to avoid the cascad generation but to get the input for evaporation etc
323
324	G4CascadeMomentum momentum_out;
325	G4InuclNuclei outgoing_nuclei(169, 69);
326
327	outgoing_nuclei.setMomentum(momentum_out);
328	outgoing_nuclei.setEnergy();
329	outgoing_nuclei.setExitationEnergy(150.0);
330
331	G4ExitonConfiguration theExitonConfiguration(3.0, 3.0, 5.0, 6.0);
332
333	outgoing_nuclei.setExitonConfiguration(theExitonConfiguration);
334	output.addTargetFragment(outgoing_nuclei);
335
336	return output;
337
338	/*
339	G4InuclElementaryParticle* bparticle = dynamic_cast<G4InuclElementaryParticle*>
340	(bullet);
341	G4InuclNuclei* tnuclei = dynamic_cast<G4InuclNuclei*>(target);
342	output.addOutgoingParticle(*bparticle);
343	output.addTargetFragment(*tnuclei);
344	return output;
345	*/
346
347	#endif
348
349	if (verboseLevel > 3) {
350	G4cout << " IntraNucleiCascader-> no inelastic interaction after " << itry_max << " attempts "
351	<< G4endl;
352	}
353
354	output.trivialise(bullet, target);
355
356	return output;
357	}
358
359	G4bool G4IntraNucleiCascader::goodCase(G4double a,
360	G4double z,
361	G4double eexs,
362	G4double ein) const {
363
364	if (verboseLevel > 3) {
365	G4cout << " >>> G4IntraNucleiCascader::goodCase" << G4endl;
366	}
367
368	const G4double eexs_cut = 0.0001;
369	const G4double reason_cut = 7.0;
370	const G4double ediv_cut = 5.0;
371
372	G4bool good = false;
373
374	if (eexs > eexs_cut) {
375	G4double eexs_max0z = 1000.0 * ein / ediv_cut;
376	G4double dm = bindingEnergy(a, z);
377	G4double eexs_max = eexs_max0z > reason_cutdm ? eexs_max0z : reason_cut dm;
378
379	if(eexs < eexs_max) good = true;
380
381	if (verboseLevel > 3) {
382	G4cout << " eexs " << eexs << " max " << eexs_max << " dm " << dm << G4endl;
383	}
384
385	};
386
387	return good;
388	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: