Context Navigation

G4InclAblaCascadeInterface.cc @ 1358

Last change on this file since 1358 was 1347, checked in by garnier, 14 years ago
geant4 tag 9.4
File size: 22.1 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	// $Id: G4InclAblaCascadeInterface.cc,v 1.20 2010/11/17 20:19:09 kaitanie Exp $
27	// Translation of INCL4.2/ABLA V3
28	// Pekka Kaitaniemi, HIP (translation)
29	// Christelle Schmidt, IPNL (fission code)
30	// Alain Boudard, CEA (contact person INCL/ABLA)
31	// Aatos Heikkinen, HIP (project coordination)
32
33	//#define DEBUGINCL 1
34
35	#include "G4InclAblaCascadeInterface.hh"
36	#include "G4FermiBreakUp.hh"
37	#include "math.h"
38	#include "G4GenericIon.hh"
39	#include "CLHEP/Random/Random.h"
40
41	G4InclAblaCascadeInterface::G4InclAblaCascadeInterface(const G4String& nam)
42	:G4VIntraNuclearTransportModel(nam)
43	{
44	hazard = new G4Hazard();
45	const G4long* table_entry = CLHEP::HepRandom::getTheSeeds(); // Get random seed from CLHEP.
46	hazard->ial = (*table_entry);
47
48	varntp = new G4VarNtp();
49	calincl = 0;
50	ws = new G4Ws();
51	mat = new G4Mat();
52	incl = new G4Incl(hazard, calincl, ws, mat, varntp);
53	if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
54	incl->setUseFermiBreakUp(true);
55	}
56
57	verboseLevel = 0;
58	}
59
60	G4InclAblaCascadeInterface::~G4InclAblaCascadeInterface()
61	{
62	delete hazard;
63	delete varntp;
64	delete ws;
65	delete mat;
66	delete incl;
67	}
68
69	G4HadFinalState* G4InclAblaCascadeInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
70	{
71	G4int maxTries = 200;
72
73	G4int particleI;
74	G4int bulletType = 0;
75
76	// Print diagnostic messages: 0 = silent, 1 and 2 = verbose
77	verboseLevel = 0;
78
79	// Increase the event number:
80	eventNumber++;
81
82	if (verboseLevel > 1) {
83	G4cout << " >>> G4InclAblaCascadeInterface::ApplyYourself called" << G4endl;
84	}
85
86	if(verboseLevel > 1) {
87	G4cout <<"G4InclAblaCascadeInterface: Now processing INCL4 event number:" << eventNumber << G4endl;
88	}
89
90	#ifdef DEBUGINCL
91	G4cout <<"Bullet energy = " << bulletE / MeV << G4endl;
92	#endif
93
94	G4double eKin;
95	G4double momx = 0.0, momy = 0.0, momz = 0.0;
96	G4DynamicParticle *cascadeParticle = 0;
97	G4ParticleDefinition *aParticleDefinition = 0;
98
99	G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp();
100	G4FragmentVector *theFermiBreakupResult = 0;
101	G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
102
103	// INCL assumes the projectile particle is going in the direction of
104	// the Z-axis. Here we construct proper rotation to convert the
105	// momentum vectors of the outcoming particles to the original
106	// coordinate system.
107	G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
108	G4LorentzRotation toZ;
109	toZ.rotateZ(-projectileMomentum.phi());
110	toZ.rotateY(-projectileMomentum.theta());
111	G4LorentzRotation toLabFrame = toZ.inverse();
112
113	theResult.Clear(); // Make sure the output data structure is clean.
114
115	calincl = new G4InclInput(aTrack, theNucleus, false);
116	incl->setInput(calincl);
117
118	// G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
119	// calincl->printInfo();
120
121	#ifdef DEBUGINCL
122	G4int baryonBullet = 0, chargeBullet = 0;
123	if(bulletType == proton \|\| bulletType == neutron) baryonBullet = 1;
124	if(bulletType == proton \|\| bulletType == pionPlus) chargeBullet = 1;
125	if(bulletType == pionMinus) chargeBullet = -1;
126	G4int baryonNumber = int(std::floor(targetA)) + baryonBullet;
127	G4int chargeNumber = int(std::floor(targetZ)) + chargeBullet;
128	G4double mass = aTrack.GetDefinition()->GetPDGMass();
129	G4double amass = theNucleus.AtomicMass(targetA, targetZ);
130	G4double eKinSum = bulletE;
131	G4LorentzVector labv = G4LorentzVector(0.0, 0.0, std::sqrt(bulletE(bulletE + 2.mass)), bulletE + mass + amass);
132	G4LorentzVector labvA = G4LorentzVector(0.0, 0.0, 0.0, 0.0);
133	G4cout <<"Energy in the beginning = " << labv.e() / MeV << G4endl;
134	#endif
135
136	// Check wheter the input is acceptable.
137	if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
138	ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon
139	ws->xfoisa = 8;
140	ws->npaulstr = 0;
141
142	int nTries = 0;
143	varntp->ntrack = 0;
144
145	mat->nbmat = 1;
146	mat->amat[0] = int(calincl->targetA());
147	mat->zmat[0] = int(calincl->targetZ());
148
149	incl->initIncl(true);
150
151	while((varntp->ntrack <= 0) && (nTries < maxTries)) { // Loop until we produce real cascade
152	nTries++;
153	if(verboseLevel > 1) {
154	G4cout <<"G4InclAblaCascadeInterface: Try number = " << nTries << G4endl;
155	}
156	incl->processEventInclAbla(calincl, eventNumber);
157
158	if(verboseLevel > 1) {
159	G4cout <<"G4InclAblaCascadeInterface: number of tracks = " << varntp->ntrack <<G4endl;
160	}
161	}
162
163	if(verboseLevel > 1) {
164	/**
165	* Diagnostic output
166	*/
167	G4cout <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
168	G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
169
170	G4cout <<"G4InclAblaCascadeInterface: Target A: " << calincl->targetA() << G4endl;
171	G4cout <<"G4InclAblaCascadeInterface: Target Z: " << calincl->targetZ() << G4endl;
172
173	if(verboseLevel > 3) {
174	diagdata <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
175	diagdata <<"G4InclAblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
176
177	diagdata <<"G4InclAblaCascadeInterface: Target A: " << calincl->targetA() << G4endl;
178	diagdata <<"G4InclAblaCascadeInterface: Target Z: " << calincl->targetZ() << G4endl;
179	}
180
181	}
182
183	// Check whether a valid cascade was produced.
184	// If not return the original bullet particle with the same momentum.
185	if(varntp->ntrack <= 0) {
186	if(verboseLevel > 1) {
187	G4cout <<"WARNING G4InclAblaCascadeInterface: No cascade. Returning original particle with original momentum." << G4endl;
188	G4cout <<"\t Reached maximum trials of 200 to produce inelastic scattering." << G4endl;
189	}
190
191	theResult.SetStatusChange(stopAndKill);
192
193	G4int bulletType = calincl->bulletType();
194	aParticleDefinition = G4InclInput::getParticleDefinition(bulletType);
195
196	if(aParticleDefinition != 0) {
197	cascadeParticle = new G4DynamicParticle();
198	cascadeParticle->SetDefinition(aParticleDefinition);
199	cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
200	theResult.AddSecondary(cascadeParticle);
201	}
202	}
203
204	// Convert INCL4 output to Geant4 compatible data structures.
205	// Elementary particles are converted to G4DynamicParticle.
206	theResult.SetStatusChange(stopAndKill);
207
208	#ifdef DEBUGINCL
209	G4cout << "E [MeV]" << std::setw(12)
210	<< " Ekin [MeV]" << std::setw(12)
211	<< "Px [MeV]" << std::setw(12)
212	<< " Py [MeV]" << std::setw(12)
213	<< "Pz [MeV]" << std::setw(12)
214	<< "Pt [MeV]" << std::setw(12)
215	<< "A" << std::setw(12)
216	<< "Z" << G4endl;
217	#endif
218
219	for(particleI = 0; particleI < varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.
220	// Get energy/momentum and construct momentum vector in INCL4 coordinates.
221	momx = varntp->plab[particleI]std::sin(varntp->tetlab[particleI]CLHEP::pi/180.0)std::cos(varntp->philab[particleI]CLHEP::pi/180.0)*MeV;
222	momy = varntp->plab[particleI]std::sin(varntp->tetlab[particleI]CLHEP::pi/180.0)std::sin(varntp->philab[particleI]CLHEP::pi/180.0)*MeV;
223	momz = varntp->plab[particleI]std::cos(varntp->tetlab[particleI]CLHEP::pi/180.0)*MeV;
224
225	eKin = varntp->enerj[particleI] * MeV;
226
227	G4ThreeVector momDirection(momx, momy, momz); // Direction of the particle.
228	momDirection = momDirection.unit();
229	if(verboseLevel > 2) {
230	G4cout <<"G4InclAblaCascadeInterface: " << G4endl;
231	G4cout <<"A = " << varntp->avv[particleI] << " Z = " << varntp->zvv[particleI] << G4endl;
232	G4cout <<"eKin = " << eKin << " MeV" << G4endl;
233	G4cout <<"px = " << momDirection.x() << " py = " << momDirection.y() <<" pz = " << momDirection.z() << G4endl;
234	}
235
236	G4int particleIdentified = 0; // Check particle ID.
237
238	if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) { // Proton
239	cascadeParticle =
240	new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
241	particleIdentified++;
242	#ifdef DEBUGINCL
243	baryonNumber--;
244	chargeNumber--;
245	#endif
246	}
247
248	if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) { // Neutron
249	cascadeParticle =
250	new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
251	particleIdentified++;
252	#ifdef DEBUGINCL
253	baryonNumber--;
254	#endif
255	}
256
257	if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) { // PionPlus
258	cascadeParticle =
259	new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
260	particleIdentified++;
261	#ifdef DEBUGINCL
262	chargeNumber--;
263	#endif
264	}
265
266	if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) { // PionZero
267	cascadeParticle =
268	new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
269	particleIdentified++;
270	}
271
272	if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) { // PionMinus
273	cascadeParticle =
274	new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
275	particleIdentified++;
276	#ifdef DEBUGINCL
277	chargeNumber++;
278	#endif
279	}
280
281	if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
282	G4ParticleDefinition * aIonDef = 0;
283
284	G4int A = G4int(varntp->avv[particleI]);
285	G4int Z = G4int(varntp->zvv[particleI]);
286	G4double excitationE = G4double(varntp->exini) * MeV;
287
288	if(verboseLevel > 1) {
289	G4cout <<"Finding ion: A = " << A << " Z = " << Z << " E* = " << excitationE/MeV << G4endl;
290	}
291	aIonDef = theTableOfParticles->GetIon(Z, A, excitationE);
292
293	if(aIonDef == 0) {
294	if(verboseLevel > 1) {
295	G4cout <<"G4InclAblaCascadeInterface: " << G4endl;
296	G4cout <<"FATAL ERROR: aIonDef = 0" << G4endl;
297	G4cout <<"A = " << A << " Z = " << Z << " E* = " << excitationE << G4endl;
298	}
299	}
300
301	if(aIonDef != 0) { // If the ion was identified add it to output.
302	cascadeParticle =
303	new G4DynamicParticle(aIonDef, momDirection, eKin);
304	particleIdentified++;
305	#ifdef DEBUGINCL
306	baryonNumber = baryonNumber - A;
307	chargeNumber = chargeNumber - Z;
308	#endif
309	}
310	}
311
312	if(particleIdentified == 1) { // Particle identified properly.
313	cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
314	#ifdef DEBUGINCL
315	G4ParticleDefinition *pd = cascadeParticle->GetDefinition();
316	G4LorentzVector fm = cascadeParticle->Get4Momentum();
317	G4ThreeVector mom = cascadeParticle->GetMomentum();
318	G4double m = pd->GetPDGMass();
319	G4double p = mom.mag();
320	labv -= fm;
321	if(varntp->avv[particleI] > 1) {
322	labvA += fm;
323	}
324	G4double px = mom.x() * MeV;
325	G4double py = mom.y() * MeV;
326	G4double pz = mom.z() * MeV;
327	G4double pt = std::sqrt(pxpx+pypy);
328	G4double e = fm.e();
329	eKinSum -= cascadeParticle->GetKineticEnergy() * MeV;
330	G4double exE;
331	if(varntp->avv[particleI] > 1) {
332	exE = varntp->exini;
333	}
334	else {
335	exE = 0.0;
336	}
337	G4cout << fm.e() / MeV
338	<< std::setw(12) << cascadeParticle->GetKineticEnergy() / MeV
339	<< std::setw(12) << mom.x() / MeV
340	<< std::setw(12) << mom.y() / MeV
341	<< std::setw(12) << mom.z() / MeV
342	<< std::setw(12) << pt / MeV
343	<< std::setw(12) << varntp->avv[particleI]
344	<< std::setw(12) << varntp->zvv[particleI] << G4endl;
345	#endif
346	theResult.AddSecondary(cascadeParticle); // Put data into G4HadFinalState.
347	}
348	else { // Particle identification failed.
349	if(particleIdentified > 1) { // Particle was identified as more than one particle type.
350	if(verboseLevel > 1) {
351	G4cout <<"G4InclAblaCascadeInterface: One outcoming particle was identified as";
352	G4cout <<"more than one particle type. This is probably due to a bug in the interface." << G4endl;
353	G4cout <<"Particle A:" << varntp->avv[particleI] << "Z: " << varntp->zvv[particleI] << G4endl;
354	G4cout << "(particleIdentified =" << particleIdentified << ")" << G4endl;
355	}
356	}
357	}
358	}
359
360	// Finally do Fermi break-up if needed
361	if(varntp->needsFermiBreakup) {
362	// baryonNumberBalanceInINCL -= varntp->massini;
363	// chargeNumberBalanceInINCL -= varntp->mzini;
364	// Call Fermi Break-up
365	G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV;
366	G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV,
367	varntp->erecrem * MeV + nuclearMass);
368	G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini),
369	varntp->exini,
370	varntp->erecrem,
371	varntp->pxrem,
372	varntp->pyrem,
373	varntp->pzrem);
374	G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV,
375	momentumScaling * varntp->pzrem * MeV,
376	varntp->erecrem + nuclearMass);
377
378	// For four-momentum, baryon number and charge conservation check:
379	G4LorentzVector fourMomentumBalance = p4;
380	G4int baryonNumberBalance = G4int(varntp->massini);
381	G4int chargeBalance = G4int(varntp->mzini);
382
383	G4LorentzRotation toFragmentZ;
384	toFragmentZ.rotateZ(-p4.theta());
385	toFragmentZ.rotateY(-p4.phi());
386	G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
387	// p4 *= toFragmentZ;
388
389	G4LorentzVector p4rest = p4;
390	// p4rest.boost(-p4.boostVector());
391	if(verboseLevel > 0) {
392	G4cout <<"Cascade remnant nucleus:" << G4endl;
393	G4cout <<"p4: " << G4endl;
394	G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
395	G4cout <<" E = " << p4.e() << G4endl;
396
397	G4cout <<"p4rest: " << G4endl;
398	G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
399	G4cout <<" E = " << p4rest.e() << G4endl;
400	}
401
402	G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest);
403	theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant);
404	if(theFermiBreakupResult != 0) {
405	G4FragmentVector::iterator fragment;
406	for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) {
407	G4ParticleDefinition *theFragmentDefinition = 0;
408	if((fragment)->GetA_asInt() == 1 && (fragment)->GetZ_asInt() == 0) { // Neutron
409	theFragmentDefinition = G4Neutron::NeutronDefinition();
410	} else if ((fragment)->GetA_asInt() == 1 && (fragment)->GetZ_asInt() == 1) {
411	theFragmentDefinition = G4Proton::ProtonDefinition();
412	} else {
413	theFragmentDefinition = theTableOfParticles->GetIon((fragment)->GetZ_asInt(), (fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
414	}
415
416	if(theFragmentDefinition != 0) {
417	G4DynamicParticle theFragment = new G4DynamicParticle(theFragmentDefinition, (fragment)->GetMomentum());
418	G4LorentzVector labMomentum = theFragment->Get4Momentum();
419	// labMomentum.boost(p4.boostVector());
420	// labMomentum *= toFragmentLab;
421	// labMomentum *= toLabFrame;
422	theFragment->Set4Momentum(labMomentum);
423	fourMomentumBalance -= theFragment->Get4Momentum();
424	baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
425	chargeBalance -= theFragmentDefinition->GetAtomicNumber();
426	if(verboseLevel > 0) {
427	G4cout <<"Resulting fragment: " << G4endl;
428	G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
429	G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
430	}
431	theResult.AddSecondary(theFragment);
432	} else {
433	G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
434	G4cout <<"Resulting fragment: " << G4endl;
435	G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
436	G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
437	G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
438	G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
439	}
440	}
441	delete theFermiBreakupResult;
442	theFermiBreakupResult = 0;
443
444	if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) {
445	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
446	G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
447	G4cout <<"Vector components (px, py, pz, E) = ("
448	<< fourMomentumBalance.px() << ", "
449	<< fourMomentumBalance.py() << ", "
450	<< fourMomentumBalance.pz() << ", "
451	<< fourMomentumBalance.e() << ")" << G4endl;
452	}
453	if(baryonNumberBalance != 0) {
454	G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
455	}
456	if(chargeBalance != 0) {
457	G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
458	}
459	}
460	}
461
462	#ifdef DEBUGINCL
463	G4cout <<"--------------------------------------------------------------------------------" << G4endl;
464	G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2));
465	G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2));
466	G4cout << labv.e() / MeV << std::setw(12)
467	<< eKinSum / MeV << std::setw(12)
468	<< labv.x() / MeV << std::setw(12)
469	<< labv.y() / MeV << std::setw(12)
470	<< labv.z() / MeV << std::setw(12)
471	<< pt / MeV << std::setw(12)
472	<< baryonNumber << std::setw(12)
473	<< chargeNumber << " totals" << G4endl;
474	G4cout << " - " << std::setw(12)
475	<< " - " << std::setw(12)
476	<< labvA.x() / MeV << std::setw(12)
477	<< labvA.y() / MeV << std::setw(12)
478	<< labvA.z() / MeV << std::setw(12)
479	<< ptA / MeV << std::setw(12)
480	<< " - " << std::setw(12) << " - " << " totals ABLA" << G4endl;
481	G4cout << G4endl;
482
483	if(verboseLevel > 3) {
484	if(baryonNumber != 0) {
485	G4cout <<"WARNING G4InclCascadeInterface: Baryon number conservation violated." << G4endl;
486	G4cout <<"Baryon number balance after the event: " << baryonNumber << G4endl;
487	if(baryonNumber < 0) {
488	G4cout <<"Too many baryons produced." << G4endl;
489	} else {
490	G4cout <<"Too few baryons produced." << G4endl;
491	}
492	}
493	}
494	#endif
495
496	varntp->ntrack = 0; // Clean up the number of generated particles in the event.
497	}
498	/**
499	* Report unsupported features.
500	* (Check bullet, target, energy range)
501	*/
502	else { // If the bullet type was not recognized by the interface, it will be returned back without any interaction.
503	theResult.SetStatusChange(stopAndKill);
504
505	cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
506
507	theResult.AddSecondary(cascadeParticle);
508
509	if(verboseLevel > 1) {
510	G4cout <<"ERROR G4InclAblaCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;
511	}
512	if(verboseLevel > 3) {
513	diagdata <<"ERROR G4InclAblaCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl;
514	}
515
516	if(bulletType == 0) {
517	if(verboseLevel > 1) {
518	G4cout <<"G4InclAblaCascadeInterface: Unknown bullet type" << G4endl;
519	G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
520	}
521	if(verboseLevel > 3) {
522	diagdata <<"G4InclAblaCascadeInterface: Unknown bullet type" << G4endl;
523	diagdata <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
524	}
525	}
526
527	if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
528	if(verboseLevel > 1) {
529	G4cout <<"Unsupported target: " << G4endl;
530	G4cout <<"Target A: " << calincl->targetA() << G4endl;
531	G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
532	}
533	if(verboseLevel > 3) {
534	diagdata <<"Unsupported target: " << G4endl;
535	diagdata <<"Target A: " << calincl->targetA() << G4endl;
536	diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
537	}
538	}
539
540	if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
541	if(verboseLevel > 1) {
542	G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
543	G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
544	}
545	if(verboseLevel > 3) {
546	diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
547	}
548	}
549
550	if(verboseLevel > 3) {
551	diagdata <<"WARNING: returning the original bullet with original energy back to Geant4." << G4endl;
552	}
553	}
554
555	delete fermiBreakUp;
556	delete calincl;
557	calincl = 0;
558	return &theResult;
559	}
560
561	G4ReactionProductVector* G4InclAblaCascadeInterface::Propagate(G4KineticTrackVector* , G4V3DNucleus* ) {
562	return 0;
563	}
564
565

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

Context Navigation

source: trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc @ 1358

Download in other formats: