Context Navigation

G4InclCascadeInterface.cc@ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 15 years ago
geant4 tag 9.4
File size: 21.2 KB

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

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

Context Navigation

source: trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc@ 1347

Download in other formats: