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

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 22.0 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.16 2010/10/29 06:48:43 gunter 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
41G4InclAblaCascadeInterface::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
60G4InclAblaCascadeInterface::~G4InclAblaCascadeInterface()
61{
62 delete hazard;
63 delete varntp;
64 delete ws;
65 delete mat;
66 delete incl;
67}
68
69G4HadFinalState* 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(px*px+py*py);
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 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) {
442 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
443 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
444 G4cout <<"Vector components (px, py, pz, E) = ("
445 << fourMomentumBalance.px() << ", "
446 << fourMomentumBalance.py() << ", "
447 << fourMomentumBalance.pz() << ", "
448 << fourMomentumBalance.e() << ")" << G4endl;
449 }
450 if(baryonNumberBalance != 0) {
451 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
452 }
453 if(chargeBalance != 0) {
454 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
455 }
456 }
457 }
458
459#ifdef DEBUGINCL
460 G4cout <<"--------------------------------------------------------------------------------" << G4endl;
461 G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2));
462 G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2));
463 G4cout << labv.e() / MeV << std::setw(12)
464 << eKinSum / MeV << std::setw(12)
465 << labv.x() / MeV << std::setw(12)
466 << labv.y() / MeV << std::setw(12)
467 << labv.z() / MeV << std::setw(12)
468 << pt / MeV << std::setw(12)
469 << baryonNumber << std::setw(12)
470 << chargeNumber << " totals" << G4endl;
471 G4cout << " - " << std::setw(12)
472 << " - " << std::setw(12)
473 << labvA.x() / MeV << std::setw(12)
474 << labvA.y() / MeV << std::setw(12)
475 << labvA.z() / MeV << std::setw(12)
476 << ptA / MeV << std::setw(12)
477 << " - " << std::setw(12) << " - " << " totals ABLA" << G4endl;
478 G4cout << G4endl;
479
480 if(verboseLevel > 3) {
481 if(baryonNumber != 0) {
482 G4cout <<"WARNING G4InclCascadeInterface: Baryon number conservation violated." << G4endl;
483 G4cout <<"Baryon number balance after the event: " << baryonNumber << G4endl;
484 if(baryonNumber < 0) {
485 G4cout <<"Too many baryons produced." << G4endl;
486 } else {
487 G4cout <<"Too few baryons produced." << G4endl;
488 }
489 }
490 }
491#endif
492
493 varntp->ntrack = 0; // Clean up the number of generated particles in the event.
494 }
495 /**
496 * Report unsupported features.
497 * (Check bullet, target, energy range)
498 */
499 else { // If the bullet type was not recognized by the interface, it will be returned back without any interaction.
500 theResult.SetStatusChange(stopAndKill);
501
502 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
503
504 theResult.AddSecondary(cascadeParticle);
505
506 if(verboseLevel > 1) {
507 G4cout <<"ERROR G4InclAblaCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;
508 }
509 if(verboseLevel > 3) {
510 diagdata <<"ERROR G4InclAblaCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl;
511 }
512
513 if(bulletType == 0) {
514 if(verboseLevel > 1) {
515 G4cout <<"G4InclAblaCascadeInterface: Unknown bullet type" << G4endl;
516 G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
517 }
518 if(verboseLevel > 3) {
519 diagdata <<"G4InclAblaCascadeInterface: Unknown bullet type" << G4endl;
520 diagdata <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
521 }
522 }
523
524 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
525 if(verboseLevel > 1) {
526 G4cout <<"Unsupported target: " << G4endl;
527 G4cout <<"Target A: " << calincl->targetA() << G4endl;
528 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
529 }
530 if(verboseLevel > 3) {
531 diagdata <<"Unsupported target: " << G4endl;
532 diagdata <<"Target A: " << calincl->targetA() << G4endl;
533 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
534 }
535 }
536
537 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
538 if(verboseLevel > 1) {
539 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
540 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
541 }
542 if(verboseLevel > 3) {
543 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
544 }
545 }
546
547 if(verboseLevel > 3) {
548 diagdata <<"WARNING: returning the original bullet with original energy back to Geant4." << G4endl;
549 }
550 }
551
552 delete fermiBreakUp;
553 delete calincl;
554 calincl = 0;
555 return &theResult;
556}
557
558G4ReactionProductVector* G4InclAblaCascadeInterface::Propagate(G4KineticTrackVector* , G4V3DNucleus* ) {
559 return 0;
560}
561
562
Note: See TracBrowser for help on using the repository browser.