Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4CascadeInterface.cc @ 1196

Last change on this file since 1196 was 1196, checked in by garnier, 15 years ago
update CVS release candidate geant4.9.3.01
File size: 17.3 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	#include "G4CascadeInterface.hh"
28	#include "globals.hh"
29	#include "G4DynamicParticleVector.hh"
30	#include "G4IonTable.hh"
31	#include "G4InuclCollider.hh"
32	// #include "G4IntraNucleiCascader.hh"
33	#include "G4InuclElementaryParticle.hh"
34	#include "G4InuclNuclei.hh"
35	#include "G4InuclParticle.hh"
36	#include "G4CollisionOutput.hh"
37	#include "G4V3DNucleus.hh"
38	#include "G4Track.hh"
39	#include "G4Nucleus.hh"
40	#include "G4NucleiModel.hh"
41	#include "G4LorentzRotation.hh"
42
43
44	//#define BERTDEV 1 // A flag to activate a development version of Bertini cascade
45
46	typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
47	typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
48
49	G4CascadeInterface::G4CascadeInterface(const G4String& nam)
50	:G4VIntraNuclearTransportModel(nam), verboseLevel(0)
51	{
52	if (verboseLevel > 3) {
53	G4cout << " >>> G4CascadeInterface::G4CascadeInterface" << G4endl;
54	}
55
56	}
57
58
59	G4CascadeInterface::~G4CascadeInterface()
60	{}
61
62	G4ReactionProductVector* G4CascadeInterface::Propagate(G4KineticTrackVector* ,
63	G4V3DNucleus* ) {
64	return 0;
65	}
66
67	// #define debug_G4CascadeInterface
68
69	G4HadFinalState* G4CascadeInterface::ApplyYourself(const G4HadProjectile& aTrack,
70	G4Nucleus& theNucleus) {
71	#ifdef debug_G4CascadeInterface
72	static G4int counter(0);
73	counter++;
74	G4cerr << "Reaction number "<< counter << " "<<aTrack.GetDynamicParticle()->GetDefinition()->GetParticleName()<<" "<< aTrack.GetDynamicParticle()->GetKineticEnergy()<<G4endl;
75	#endif
76
77	theResult.Clear();
78
79	if (verboseLevel > 3) {
80	G4cout << " >>> G4CascadeInterface::ApplyYourself" << G4endl;
81	};
82
83	G4double eInit = 0.0;
84	G4double eTot = 0.0;
85	G4double sumBaryon = 0.0;
86	G4double sumEnergy = 0.0;
87
88	// Make conversion between native Geant4 and Bertini cascade classes.
89	// NOTE: Geant4 units are MeV = 1 and GeV = 1000. Cascade code by default use GeV = 1.
90
91	enum particleType { nuclei = 0, proton = 1, neutron = 2, pionPlus = 3,
92	pionMinus = 5, pionZero = 7, photon = 10,
93	kaonPlus = 11, kaonMinus = 13, kaonZero = 15,
94	kaonZeroBar = 17, lambda = 21, sigmaPlus = 23,
95	sigmaZero = 25, sigmaMinus = 27, xiZero = 29, xiMinus = 31 };
96
97	G4int bulletType = 0;
98
99	// Coding particles
100	if (aTrack.GetDefinition() == G4Proton::Proton() ) bulletType = proton;
101	if (aTrack.GetDefinition() == G4Neutron::Neutron() ) bulletType = neutron;
102	if (aTrack.GetDefinition() == G4PionPlus::PionPlus() ) bulletType = pionPlus;
103	if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;
104	if (aTrack.GetDefinition() == G4PionZero::PionZero() ) bulletType = pionZero;
105	if (aTrack.GetDefinition() == G4Gamma::Gamma() ) bulletType = photon;
106	if (aTrack.GetDefinition() == G4KaonPlus::KaonPlus() ) bulletType = kaonPlus;
107	if (aTrack.GetDefinition() == G4KaonMinus::KaonMinus() ) bulletType = kaonMinus;
108	if (aTrack.GetDefinition() == G4Lambda::Lambda() ) bulletType = lambda;
109	if (aTrack.GetDefinition() == G4SigmaPlus::SigmaPlus() ) bulletType = sigmaPlus;
110	if (aTrack.GetDefinition() == G4SigmaZero::SigmaZero() ) bulletType = sigmaZero;
111	if (aTrack.GetDefinition() == G4SigmaMinus::SigmaMinus() ) bulletType = sigmaMinus;
112	if (aTrack.GetDefinition() == G4XiZero::XiZero() ) bulletType = xiZero;
113	if (aTrack.GetDefinition() == G4XiMinus::XiMinus() ) bulletType = xiMinus;
114
115	if (aTrack.GetDefinition() == G4KaonZeroLong::KaonZeroLong() \|\|
116	aTrack.GetDefinition() == G4KaonZeroShort::KaonZeroShort() ) {
117	if (G4UniformRand() > 0.5) {
118	bulletType = kaonZero;
119	} else {
120	bulletType = kaonZeroBar;
121	}
122	}
123
124	// Code momentum and energy.
125	G4double px,py,pz;
126	px=aTrack.Get4Momentum().px() / GeV;
127	py=aTrack.Get4Momentum().py() / GeV;
128	pz=aTrack.Get4Momentum().pz() / GeV;
129
130	G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
131	G4LorentzRotation toZ;
132	toZ.rotateZ(-projectileMomentum.phi());
133	toZ.rotateY(-projectileMomentum.theta());
134	G4LorentzRotation toLabFrame = toZ.inverse();
135
136	G4CascadeMomentum momentumBullet;
137	momentumBullet[0] =0.;
138	momentumBullet[1] =0;
139	momentumBullet[2] =0;
140	momentumBullet[3] =std::sqrt(pxpx+pypy+pz*pz);
141
142	G4InuclElementaryParticle * bullet = new G4InuclElementaryParticle(momentumBullet, bulletType);
143
144	sumEnergy = bullet->getKineticEnergy(); // In GeV
145
146	if (bulletType == proton \|\| bulletType == neutron \|\| bulletType == lambda \|\|
147	bulletType == sigmaPlus \|\| bulletType == sigmaZero \|\| bulletType == sigmaMinus \|\|
148	bulletType == xiZero \|\| bulletType == xiMinus) {
149
150	sumBaryon += 1;
151	}
152
153	// Set target
154	G4InuclNuclei* target = 0;
155	G4InuclParticle* targetH = 0;
156	// and outcoming particles
157	G4DynamicParticle* cascadeParticle = 0;
158
159	G4CascadeMomentum targetMomentum;
160
161	G4double theNucleusA = theNucleus.GetN();
162
163	if ( !(G4int(theNucleusA) == 1) ) {
164	target = new G4InuclNuclei(targetMomentum,
165	theNucleusA,
166	theNucleus.GetZ());
167	target->setEnergy();
168
169	const G4CascadeMomentum& bmom = bullet->getMomentum();
170	eInit = std::sqrt(bmom[0] * bmom[0]);
171	const G4CascadeMomentum& tmom = target->getMomentum();
172	eInit += std::sqrt(tmom[0] * tmom[0]);
173
174	sumBaryon += theNucleusA;
175
176	if (verboseLevel > 2) {
177	G4cout << "Bullet: " << G4endl;
178	bullet->printParticle();
179	}
180	if (verboseLevel > 2) {
181	G4cout << "Target: " << G4endl;
182	target->printParticle();
183	}
184	}
185
186	G4CollisionOutput output;
187
188	// Colliders initialisation
189	// G4ElementaryParticleCollider* colep = new G4ElementaryParticleCollider;
190	// G4IntraNucleiCascader* inc = new G4IntraNucleiCascader; // the actual cascade
191	// inc->setInteractionCase(1); // Interaction type is particle with nuclei.
192	inc.setInteractionCase(1); // Interaction type is particle with nuclei.
193
194	// G4NonEquilibriumEvaporator* noneq = new G4NonEquilibriumEvaporator;
195	// G4EquilibriumEvaporator* eqil = new G4EquilibriumEvaporator;
196	// G4Fissioner* fiss = new G4Fissioner;
197	// G4BigBanger* bigb = new G4BigBanger;
198	G4InuclCollider* collider = new G4InuclCollider(&colep, &inc, &noneq, &eqil, &fiss, &bigb);
199
200	G4int maxTries = 100; // maximum tries for inelastic collision to avoid infinite loop
201	G4int nTries = 0; // try counter
202
203	#ifdef BERTDEV
204	G4int coulombOK =0; // flag for correct Coulomb barrier
205	#endif
206	if (G4int(theNucleusA) == 1) { // special treatment for target H(1,1) (proton)
207
208	targetH = new G4InuclElementaryParticle(targetMomentum, 1);
209
210	G4float cutElastic[32];
211
212	cutElastic[proton ] = 1.0; // GeV
213	cutElastic[neutron ] = 1.0;
214	cutElastic[pionPlus ] = 0.6;
215	cutElastic[pionMinus] = 0.2;
216	cutElastic[pionZero ] = 0.2;
217	cutElastic[kaonPlus ] = 0.5;
218	cutElastic[kaonMinus] = 0.5;
219	cutElastic[kaonMinus] = 0.5;
220	cutElastic[kaonZero] = 0.5;
221	cutElastic[kaonZeroBar] = 0.5;
222	cutElastic[lambda] = 1.0;
223	cutElastic[sigmaPlus] = 1.0;
224	cutElastic[sigmaZero] = 1.0;
225	cutElastic[sigmaMinus] = 1.0;
226	cutElastic[xiZero] = 1.0;
227	cutElastic[xiMinus] = 1.0;
228
229	if (momentumBullet[3] > cutElastic[bulletType]) { // inelastic collision possible
230
231	do { // we try to create inelastic interaction
232	output = collider->collide(bullet, targetH);
233	nTries++;
234	} while(
235	(nTries < maxTries) &&
236	(output.getOutgoingParticles().size() == 2 && // elastic: bullet + p = H(1,1) coming out
237	(output.getOutgoingParticles().begin()->type() == bulletType \|\|
238	output.getOutgoingParticles().begin()->type() == proton)
239	)
240	);
241	} else { // only elastic collision is energetically possible
242	output = collider->collide(bullet, targetH);
243	}
244
245	sumBaryon += 1;
246
247	const G4CascadeMomentum& bmom = bullet->getMomentum();
248	eInit = std::sqrt(bmom[0] * bmom[0]);
249	const G4CascadeMomentum& tmom = targetH->getMomentum();
250	eInit += std::sqrt(tmom[0] * tmom[0]);
251
252	if (verboseLevel > 2) {
253	G4cout << "Target: " << G4endl;
254	targetH->printParticle();
255	}
256
257	} else { // treat all other targets excepet H(1,1)
258
259	do // we try to create inelastic interaction
260	{
261	#ifdef BERTDEV
262	coulombOK=0; // by default coulomb analysis is OK
263	#endif
264	output = collider->collide(bullet, target );
265	nTries++;
266
267	#ifdef BERTDEV
268	G4double coulumbBarrier = 8.7 * MeV;
269	std::vector<G4InuclElementaryParticle> p= output.getOutgoingParticles();
270	if(!p.empty()) {
271	for( particleIterator ipart = p.begin(); ipart != p.end(); ipart++) {
272	if (ipart->type() == proton) {
273	G4double e = ipart->getKineticEnergy()*GeV;
274	if (e < coulumbBarrier) coulombOK= 1; // If event with coulomb barrier violation detected -> retry
275	// G4cout << "///AH "<< e << "" << coulumbBarrier <<G4endl;
276	}
277	}
278	}
279	} while(
280	(nTries < maxTries) && // conditions for next try
281	(coulombOK==1) &&
282	((output.getOutgoingParticles().size() + output.getNucleiFragments().size()) > 2.5) &&
283	(output.getOutgoingParticles().size()!=0)
284	);
285	#else
286
287	} while(
288	(nTries < maxTries) &&
289	(output.getOutgoingParticles().size() + output.getNucleiFragments().size() < 2.5) &&
290	(output.getOutgoingParticles().size()!=0) &&
291	(output.getOutgoingParticles().begin()->type()==bullet->type())
292	);
293	#endif
294	}
295
296	if (verboseLevel > 1)
297	{
298	G4cout << " Cascade output: " << G4endl;
299	output.printCollisionOutput();
300	}
301
302	// Convert cascade data to use hadronics interface
303	std::vector<G4InuclNuclei> nucleiFragments = output.getNucleiFragments();
304	std::vector<G4InuclElementaryParticle> particles = output.getOutgoingParticles();
305
306	theResult.SetStatusChange(stopAndKill);
307
308	if (!particles.empty()) {
309	particleIterator ipart;
310	G4int outgoingParticle;
311
312	for (ipart = particles.begin(); ipart != particles.end(); ipart++) {
313	outgoingParticle = ipart->type();
314	const G4CascadeMomentum& mom = ipart->getMomentum();
315	eTot += std::sqrt(mom[0] * mom[0]);
316
317	G4double ekin = ipart->getKineticEnergy() * GeV;
318	G4ThreeVector aMom(mom[1], mom[2], mom[3]);
319	aMom = aMom.unit();
320
321	if (ipart->baryon() ) {
322	sumBaryon -= 1;
323	}
324
325	sumEnergy -= ekin / GeV;
326
327	switch(outgoingParticle) {
328
329	case proton:
330	#ifdef debug_G4CascadeInterface
331	G4cerr << "proton " << counter << " " << aMom << " " << ekin << G4endl;
332	#endif
333	cascadeParticle =
334	new G4DynamicParticle(G4Proton::ProtonDefinition(), aMom, ekin);
335	break;
336
337	case neutron:
338
339	#ifdef debug_G4CascadeInterface
340	G4cerr << "neutron "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
341	#endif
342	cascadeParticle =
343	new G4DynamicParticle(G4Neutron::NeutronDefinition(), aMom, ekin);
344	break;
345
346	case pionPlus:
347	cascadeParticle =
348	new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), aMom, ekin);
349
350	#ifdef debug_G4CascadeInterface
351	G4cerr << "pionPlus "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
352	#endif
353	break;
354
355	case pionMinus:
356	cascadeParticle =
357	new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), aMom, ekin);
358
359	#ifdef debug_G4CascadeInterface
360	G4cerr << "pionMinus "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
361	#endif
362	break;
363
364	case pionZero:
365	cascadeParticle =
366	new G4DynamicParticle(G4PionZero::PionZeroDefinition(), aMom, ekin);
367
368	#ifdef debug_G4CascadeInterface
369	G4cerr << "pionZero "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
370	#endif
371	break;
372
373	case photon:
374	cascadeParticle =
375	new G4DynamicParticle(G4Gamma::Gamma(), aMom, ekin);
376
377	#ifdef debug_G4CascadeInterface
378	G4cerr << "photon "<< counter<<" "<<aMom<<" "<< ekin<<G4endl;
379	#endif
380	break;
381
382
383	case kaonPlus:
384	cascadeParticle =
385	new G4DynamicParticle(G4KaonPlus::KaonPlusDefinition(), aMom, ekin);
386	break;
387
388	case kaonMinus:
389	cascadeParticle =
390	new G4DynamicParticle(G4KaonMinus::KaonMinusDefinition(), aMom, ekin);
391	break;
392
393	case kaonZero:
394	if (G4UniformRand() > 0.5) {
395	cascadeParticle = new G4DynamicParticle(
396	G4KaonZeroLong::KaonZeroLongDefinition(),
397	aMom, ekin);
398	} else {
399	cascadeParticle = new G4DynamicParticle(
400	G4KaonZeroShort::KaonZeroShortDefinition(),
401	aMom, ekin);
402	}
403	break;
404
405	case kaonZeroBar:
406	if (G4UniformRand() > 0.5) {
407	cascadeParticle = new G4DynamicParticle(
408	G4KaonZeroLong::KaonZeroLongDefinition(),
409	aMom, ekin);
410	} else {
411	cascadeParticle = new G4DynamicParticle(
412	G4KaonZeroShort::KaonZeroShortDefinition(),
413	aMom, ekin);
414	}
415	break;
416
417	case lambda:
418	cascadeParticle =
419	new G4DynamicParticle(G4Lambda::LambdaDefinition(), aMom, ekin);
420	break;
421
422	case sigmaPlus:
423	cascadeParticle =
424	new G4DynamicParticle(G4SigmaPlus::SigmaPlusDefinition(), aMom, ekin);
425	break;
426
427	case sigmaZero:
428	cascadeParticle =
429	new G4DynamicParticle(G4SigmaZero::SigmaZeroDefinition(), aMom, ekin);
430	break;
431
432	case sigmaMinus:
433	cascadeParticle =
434	new G4DynamicParticle(G4SigmaMinus::SigmaMinusDefinition(), aMom, ekin);
435	break;
436
437	case xiZero:
438	cascadeParticle =
439	new G4DynamicParticle(G4XiZero::XiZeroDefinition(), aMom, ekin);
440	break;
441
442	case xiMinus:
443	cascadeParticle =
444	new G4DynamicParticle(G4XiMinus::XiMinusDefinition(), aMom, ekin);
445	break;
446
447	default:
448	G4cout << " ERROR: G4CascadeInterface::Propagate undefined particle type"
449	<< G4endl;
450	}
451
452	cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
453	theResult.AddSecondary(cascadeParticle);
454	}
455	}
456
457	// get nuclei fragments
458	G4DynamicParticle * aFragment = 0;
459	G4ParticleDefinition * aIonDef = 0;
460	G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
461
462	if (!nucleiFragments.empty()) {
463	nucleiIterator ifrag;
464
465	for (ifrag = nucleiFragments.begin(); ifrag != nucleiFragments.end(); ifrag++)
466	{
467	G4double eKin = ifrag->getKineticEnergy() * GeV;
468	const G4CascadeMomentum& mom = ifrag->getMomentum();
469	eTot += std::sqrt(mom[0] * mom[0]);
470
471	G4ThreeVector aMom(mom[1], mom[2], mom[3]);
472	aMom = aMom.unit();
473
474	// hpw @@@ ==> Should be zero: G4double fragmentExitation = ifrag->getExitationEnergyInGeV();
475
476	if (verboseLevel > 2) {
477	G4cout << " Nuclei fragment: " << G4endl;
478	ifrag->printParticle();
479	}
480
481	G4int A = G4int(ifrag->getA());
482	G4int Z = G4int(ifrag->getZ());
483	aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
484
485	aFragment = new G4DynamicParticle(aIonDef, aMom, eKin);
486
487	sumBaryon -= A;
488	sumEnergy -= eKin / GeV;
489
490	aFragment->Set4Momentum(aFragment->Get4Momentum()*=toLabFrame);
491	theResult.AddSecondary(aFragment);
492	}
493	}
494
495	if (verboseLevel > 2) {
496	if (sumBaryon != 0) {
497	G4cout << "ERROR: no baryon number conservation, sum of baryons = "
498	<< sumBaryon << G4endl;
499	}
500
501	if (sumEnergy > 0.01 ) {
502	G4cout << "Kinetic energy conservation violated by "
503	<< sumEnergy << " GeV" << G4endl;
504	}
505
506	G4cout << "Total energy conservation at level ~"
507	<< (eInit - eTot) * GeV << " MeV" << G4endl;
508
509	if (sumEnergy < -5.0e-5 ) { // 0.05 MeV
510	G4cout << "FATAL ERROR: energy created "
511	<< sumEnergy * GeV << " MeV" << G4endl;
512	}
513	}
514
515	delete bullet;
516	// delete inc;
517	delete collider;
518
519	if(target != 0) delete target;
520	if(targetH != 0) delete targetH;
521	// if(cascadeParticle != 0) delete cascadeParticle;
522	// if(aFragment != 0) delete aFragment;
523
524	return &theResult;
525	}

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

Download in other formats: