Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4IBertini.cc @ 962

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

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

Download in other formats: