Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4CascadeElasticInterface.cc @ 1315

Last change on this file since 1315 was 962, checked in by garnier, 15 years ago
update processes
File size: 12.8 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	#include "G4CascadeElasticInterface.hh"
27	#include "globals.hh"
28	#include "G4DynamicParticleVector.hh"
29	#include "G4IonTable.hh"
30	#include "G4InuclCollider.hh"
31	#include "G4IntraNucleiCascader.hh"
32	#include "G4ElementaryParticleCollider.hh"
33	#include "G4NonEquilibriumEvaporator.hh"
34	#include "G4EquilibriumEvaporator.hh"
35	#include "G4Fissioner.hh"
36	#include "G4BigBanger.hh"
37	#include "G4InuclElementaryParticle.hh"
38	#include "G4InuclNuclei.hh"
39	#include "G4InuclParticle.hh"
40	#include "G4CollisionOutput.hh"
41	#include "G4V3DNucleus.hh"
42	#include "G4Track.hh"
43	#include "G4Nucleus.hh"
44	#include "G4NucleiModel.hh"
45	#include "G4LorentzRotation.hh"
46
47
48	typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
49	typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
50
51	G4CascadeElasticInterface::G4CascadeElasticInterface()
52	:verboseLevel(0) {
53
54	if (verboseLevel > 3) {
55	G4cout << " >>> G4CascadeElasticInterface::G4CascadeElasticInterface" << G4endl;
56	}
57	}
58
59	G4ReactionProductVector* G4CascadeElasticInterface::Propagate(G4KineticTrackVector* ,
60	G4V3DNucleus* ) {
61	return 0;
62	}
63
64	// #define debug_G4CascadeElasticInterface
65
66	G4HadFinalState* G4CascadeElasticInterface::ApplyYourself(const G4HadProjectile& aTrack,
67	G4Nucleus& theNucleus) {
68	#ifdef debug_G4CascadeElasticInterface
69	static G4int counter(0);
70	counter++;
71	G4cerr << "Reaction number "<< counter << " "<<aTrack.GetDynamicParticle()->GetDefinition()->GetParticleName()<<" "<< aTrack.GetDynamicParticle()->GetKineticEnergy()<<G4endl;
72	#endif
73
74	theResult.Clear();
75
76	if (verboseLevel > 3) {
77	G4cout << " >>> G4CascadeElasticInterface::ApplyYourself" << G4endl;
78	};
79
80	G4double eInit = 0.0;
81	G4double eTot = 0.0;
82	G4double sumBaryon = 0.0;
83	G4double sumEnergy = 0.0;
84
85	// Make conversion between native Geant4 and Bertini cascade classes.
86	// NOTE: Geant4 units are MeV = 1 and GeV = 1000. Cascade code by default use GeV = 1.
87
88	enum particleType { nuclei = 0, proton = 1, neutron = 2, pionPlus = 3, pionMinus = 5, pionZero = 7, photon = 10 };
89
90	G4int bulletType = 0;
91
92	// Coding particles
93	if (aTrack.GetDefinition() == G4Proton::Proton() ) bulletType = proton;
94	if (aTrack.GetDefinition() == G4Neutron::Neutron() ) bulletType = neutron;
95	if (aTrack.GetDefinition() == G4PionPlus::PionPlus() ) bulletType = pionPlus;
96	if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;
97	if (aTrack.GetDefinition() == G4PionZero::PionZero() ) bulletType = pionZero;
98	if (aTrack.GetDefinition() == G4Gamma::Gamma() ) bulletType = photon;
99
100	// Code momentum and energy.
101	G4double px,py,pz;
102	px=aTrack.Get4Momentum().px() / GeV;
103	py=aTrack.Get4Momentum().py() / GeV;
104	pz=aTrack.Get4Momentum().pz() / GeV;
105
106	G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
107	G4LorentzRotation toZ;
108	toZ.rotateZ(-projectileMomentum.phi());
109	toZ.rotateY(-projectileMomentum.theta());
110	G4LorentzRotation toLabFrame = toZ.inverse();
111
112	G4CascadeMomentum momentumBullet;
113	momentumBullet[0] =0.;
114	momentumBullet[1] =0;
115	momentumBullet[2] =0;
116	momentumBullet[3] =std::sqrt(pxpx+pypy+pz*pz);
117
118	G4InuclElementaryParticle * bullet = new G4InuclElementaryParticle(momentumBullet, bulletType);
119
120	sumEnergy = bullet->getKineticEnergy(); // In GeV
121	if (bulletType == proton \|\| bulletType == neutron) {
122	sumBaryon += 1;
123	}
124
125	// Set target
126	G4InuclNuclei* target = 0;
127	G4InuclParticle* targetH = 0;
128	// and outcoming particles
129	G4DynamicParticle* cascadeParticle = 0;
130
131	G4CascadeMomentum targetMomentum;
132
133	G4double theNucleusA = theNucleus.GetN();
134
135	if ( !(G4int(theNucleusA) == 1) ) {
136	target = new G4InuclNuclei(targetMomentum,
137	theNucleusA,
138	theNucleus.GetZ());
139	target->setEnergy();
140
141	const G4CascadeMomentum& bmom = bullet->getMomentum();
142	eInit = std::sqrt(bmom[0] * bmom[0]);
143	const G4CascadeMomentum& tmom = target->getMomentum();
144	eInit += std::sqrt(tmom[0] * tmom[0]);
145
146	sumBaryon += theNucleusA;
147
148	if (verboseLevel > 2) {
149	G4cout << "Bullet: " << G4endl;
150	bullet->printParticle();
151	}
152	if (verboseLevel > 2) {
153	G4cout << "Target: " << G4endl;
154	target->printParticle();
155	}
156	}
157
158	G4CollisionOutput output;
159
160	// Colliders initialisation
161	G4ElementaryParticleCollider* colep = new G4ElementaryParticleCollider;
162	G4IntraNucleiCascader* inc = new G4IntraNucleiCascader; // the actual cascade
163	inc->setInteractionCase(1); // Interaction type is particle with nuclei.
164
165	G4NonEquilibriumEvaporator* noneq = new G4NonEquilibriumEvaporator;
166	G4EquilibriumEvaporator* eqil = new G4EquilibriumEvaporator;
167	G4Fissioner* fiss = new G4Fissioner;
168	G4BigBanger* bigb = new G4BigBanger;
169	G4InuclCollider* collider = new G4InuclCollider(colep, inc, noneq, eqil, fiss, bigb);
170
171	G4int maxTries = 10; // maximum tries for inelastic collision to avoid infinite loop
172	G4int nTries = 0; // try counter
173
174	if (G4int(theNucleusA) == 1) { // special treatment for target H(1,1) (proton)
175
176	targetH = new G4InuclElementaryParticle(targetMomentum, 1);
177
178	G4float cutElastic[8];
179	cutElastic[proton ] = 1.0; // GeV
180	cutElastic[neutron ] = 1.0;
181	cutElastic[pionPlus ] = 0.6;
182	cutElastic[pionMinus] = 0.2;
183	cutElastic[pionZero ] = 0.2;
184
185	if (momentumBullet[3] > cutElastic[bulletType]) { // inelastic collision possible
186
187	do { // we try to create inelastic interaction
188	output = collider->collide(bullet, targetH);
189	nTries++;
190	} while(
191	(nTries < maxTries) &&
192	(output.getOutgoingParticles().size() == 2 && // elastic: bullet + p = H(1,1) coming out
193	(output.getOutgoingParticles().begin()->type() == bulletType \|\|
194	output.getOutgoingParticles().begin()->type() == proton)
195	)
196	);
197
198	} else { // only elastic collision is energetically possible
199	output = collider->collide(bullet, targetH);
200	}
201
202	sumBaryon += 1;
203
204	const G4CascadeMomentum& bmom = bullet->getMomentum();
205	eInit = std::sqrt(bmom[0] * bmom[0]);
206	const G4CascadeMomentum& tmom = targetH->getMomentum();
207	eInit += std::sqrt(tmom[0] * tmom[0]);
208
209	if (verboseLevel > 2) {
210	G4cout << "Target: " << G4endl;
211	targetH->printParticle();
212	}
213
214	} else { // treat all other targets excepet H(1,1)
215
216	do // we try to create inelastic interaction
217	{
218	output = collider->collide(bullet, target );
219	nTries++;
220	} while(
221	(nTries < maxTries) &&
222	(output.getOutgoingParticles().size() + output.getNucleiFragments().size() < 2.5) &&
223	(output.getOutgoingParticles().size()!=0) &&
224	(output.getOutgoingParticles().begin()->type()==bullet->type())
225	);
226
227	}
228
229	if (verboseLevel > 1)
230	{
231	G4cout << " Cascade output: " << G4endl;
232	output.printCollisionOutput();
233	}
234
235	// Convert cascade data to use hadronics interface
236	std::vector<G4InuclNuclei> nucleiFragments = output.getNucleiFragments();
237	std::vector<G4InuclElementaryParticle> particles = output.getOutgoingParticles();
238
239	theResult.SetStatusChange(stopAndKill);
240
241	if (!particles.empty()) {
242	particleIterator ipart;
243	G4int outgoingParticle;
244
245	for (ipart = particles.begin(); ipart != particles.end(); ipart++) {
246	outgoingParticle = ipart->type();
247	const G4CascadeMomentum& mom = ipart->getMomentum();
248	eTot += std::sqrt(mom[0] * mom[0]);
249
250	G4double ekin = ipart->getKineticEnergy() * GeV;
251	G4ThreeVector aMom(mom[1], mom[2], mom[3]);
252	aMom = aMom.unit();
253
254	if (outgoingParticle == proton \|\| outgoingParticle == neutron) {
255	sumBaryon -= 1;
256	}
257
258	sumEnergy -= ekin / GeV;
259
260	switch(outgoingParticle) {
261
262	case proton:
263	#ifdef debug_G4CascadeElasticInterface
264	G4cerr << "proton " << counter << " " << aMom << " " << ekin << G4endl;
265	#endif
266	cascadeParticle =
267	new G4DynamicParticle(G4Proton::ProtonDefinition(), aMom, ekin);
268	break;
269
270	case neutron:
271
272	#ifdef debug_G4CascadeElasticInterface
273	G4cerr << "neutron "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
274	#endif
275	cascadeParticle =
276	new G4DynamicParticle(G4Neutron::NeutronDefinition(), aMom, ekin);
277	break;
278
279	case pionPlus:
280	cascadeParticle =
281	new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), aMom, ekin);
282
283	#ifdef debug_G4CascadeElasticInterface
284	G4cerr << "pionPlus "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
285	#endif
286	break;
287
288	case pionMinus:
289	cascadeParticle =
290	new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), aMom, ekin);
291
292	#ifdef debug_G4CascadeElasticInterface
293	G4cerr << "pionMinus "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
294	#endif
295	break;
296
297	case pionZero:
298	cascadeParticle =
299	new G4DynamicParticle(G4PionZero::PionZeroDefinition(), aMom, ekin);
300
301	#ifdef debug_G4CascadeElasticInterface
302	G4cerr << "pionZero "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
303	#endif
304	break;
305
306	case photon:
307	cascadeParticle =
308	new G4DynamicParticle(G4Gamma::Gamma(), aMom, ekin);
309
310	#ifdef debug_G4CascadeElasticInterface
311	G4cerr << "photon "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
312	#endif
313	break;
314
315	default:
316	G4cout << " ERROR: G4CascadeElasticInterface::Propagate undefined particle type"
317	<< G4endl;
318	}
319
320	cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
321	theResult.AddSecondary(cascadeParticle);
322	}
323	}
324
325	// get nuclei fragments
326	G4DynamicParticle * aFragment = 0;
327	G4ParticleDefinition * aIonDef = 0;
328	G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
329
330	if (!nucleiFragments.empty()) {
331	nucleiIterator ifrag;
332
333	for (ifrag = nucleiFragments.begin(); ifrag != nucleiFragments.end(); ifrag++)
334	{
335	G4double eKin = ifrag->getKineticEnergy() * GeV;
336	const G4CascadeMomentum& mom = ifrag->getMomentum();
337	eTot += std::sqrt(mom[0] * mom[0]);
338
339	G4ThreeVector aMom(mom[1], mom[2], mom[3]);
340	aMom = aMom.unit();
341
342	// hpw @@@ ==> Should be zero: G4double fragmentExitation = ifrag->getExitationEnergyInGeV();
343
344	if (verboseLevel > 2) {
345	G4cout << " Nuclei fragment: " << G4endl;
346	ifrag->printParticle();
347	}
348
349	G4int A = G4int(ifrag->getA());
350	G4int Z = G4int(ifrag->getZ());
351	aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
352
353	aFragment = new G4DynamicParticle(aIonDef, aMom, eKin);
354
355	sumBaryon -= A;
356	sumEnergy -= eKin / GeV;
357
358	aFragment->Set4Momentum(aFragment->Get4Momentum()*=toLabFrame);
359	theResult.AddSecondary(aFragment);
360	}
361	}
362
363	if (verboseLevel > 2) {
364	if (sumBaryon != 0) {
365	G4cout << "ERROR: no baryon number conservation, sum of baryons = "
366	<< sumBaryon << G4endl;
367	}
368
369	if (sumEnergy > 0.01 ) {
370	G4cout << "Kinetic energy conservation violated by "
371	<< sumEnergy << " GeV" << G4endl;
372	}
373
374	G4cout << "Total energy conservation at level ~"
375	<< (eInit - eTot) * GeV << " MeV" << G4endl;
376
377	if (sumEnergy < -5.0e-5 ) { // 0.05 MeV
378	G4cout << "FATAL ERROR: energy created "
379	<< sumEnergy * GeV << " MeV" << G4endl;
380	}
381	}
382
383	delete bullet;
384	delete colep;
385	delete inc;
386	delete noneq;
387	delete fiss;
388	delete eqil;
389	delete bigb;
390	delete collider;
391
392	if(target != 0) delete target;
393	if(targetH != 0) delete targetH;
394	// if(cascadeParticle != 0) delete cascadeParticle;
395	// if(aFragment != 0) delete aFragment;
396
397	return &theResult;
398	}

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

Download in other formats: