Context Navigation

source: trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago
update ti head
File size: 101.4 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 "globals.hh"
28	#include "G4Proton.hh"
29	#include "G4Neutron.hh"
30	#include "G4LorentzRotation.hh"
31	#include "G4BinaryCascade.hh"
32	#include "G4KineticTrackVector.hh"
33	#include "G4ReactionProductVector.hh"
34	#include "G4Track.hh"
35	#include "G4V3DNucleus.hh"
36	#include "G4Fancy3DNucleus.hh"
37	#include "G4Scatterer.hh"
38	#include "G4MesonAbsorption.hh"
39	#include "G4ping.hh"
40	#include "G4Delete.hh"
41
42	#include "G4CollisionManager.hh"
43	#include "G4Absorber.hh"
44
45	#include "G4CollisionInitialState.hh"
46	#include "G4ListOfCollisions.hh"
47	#include "G4Fragment.hh"
48	#include "G4RKPropagation.hh"
49
50	#include "G4NuclearShellModelDensity.hh"
51	#include "G4NuclearFermiDensity.hh"
52	#include "G4FermiMomentum.hh"
53
54	#include "G4PreCompoundModel.hh"
55	#include "G4ExcitationHandler.hh"
56
57	#include "G4FermiPhaseSpaceDecay.hh"
58
59	#include <algorithm>
60	#include "G4ShortLivedConstructor.hh"
61	#include <typeinfo>
62
63	// turn on general debugging info, and consistency checks
64	//#define debug_G4BinaryCascade 1
65
66	// more detailed debugging -- deprecated
67	//#define debug_1_BinaryCascade 1
68
69	// specific debuuging info per method or functionality
70	//#define debug_BIC_ApplyCollision 1
71	//#define debug_BIC_CheckPauli 1
72	//#define debug_BIC_CorrectFinalPandE 1
73	//#define debug_BIC_Propagate 1
74	//#define debug_BIC_Propagate_Excitation 1
75	//#define debug_BIC_Propagate_Collisions 1
76	//#define debug_BIC_Propagate_finals 1
77	//#define debug_BIC_DoTimeStep 1
78	//#define debug_BIC_CorrectBarionsOnBoundary 1
79	//#define debug_BIC_GetExcitationEnergy 1
80	//#define debug_BIC_FinalNucleusMomentum 1
81	//#define debug_BIC_FindFragments 1
82	//#define debug_BIC_BuildTargetList 1
83	//#define debug_BIC_FindCollision 1
84
85	// C O N S T R U C T O R S A N D D E S T R U C T O R S
86	//
87
88	G4BinaryCascade::G4BinaryCascade() :
89	G4VIntraNuclearTransportModel("Binary Cascade")
90	{
91	// initialise the resonance sector
92	G4ShortLivedConstructor ShortLived;
93	ShortLived.ConstructParticle();
94
95	theCollisionMgr = new G4CollisionManager;
96	theDecay=new G4BCDecay;
97	theLateParticle= new G4BCLateParticle;
98	theImR.push_back(theDecay);
99	G4MesonAbsorption * aAb=new G4MesonAbsorption;
100	theImR.push_back(aAb);
101	G4Scatterer * aSc=new G4Scatterer;
102	theH1Scatterer = new G4Scatterer;
103	theImR.push_back(aSc);
104
105	thePropagator = new G4RKPropagation;
106	theCurrentTime = 0.;
107	theBCminP = 45*MeV;
108	theCutOnP = 90*MeV;
109	theCutOnPAbsorb= 0*MeV; // No Absorption of slow Mesons, other than above G4MesonAbsorption
110
111	theExcitationHandler = new G4ExcitationHandler;
112	SetDeExcitation(new G4PreCompoundModel(theExcitationHandler));
113	SetMinEnergy(0.0*GeV);
114	SetMaxEnergy(10.1*GeV);
115	//PrintWelcomeMessage();
116	thePrimaryEscape = true;
117	thePrimaryType = 0;
118	}
119
120
121	G4BinaryCascade::G4BinaryCascade(const G4BinaryCascade& )
122	: G4VIntraNuclearTransportModel("Binary Cascade")
123	{
124	}
125
126
127	G4BinaryCascade::~G4BinaryCascade()
128	{
129	ClearAndDestroy(&theTargetList);
130	ClearAndDestroy(&theSecondaryList);
131	ClearAndDestroy(&theCapturedList);
132	ClearAndDestroy(&theProjectileList);
133	delete thePropagator;
134	delete theCollisionMgr;
135	std::for_each(theImR.begin(), theImR.end(), Delete<G4BCAction>());
136	delete theLateParticle;
137	delete theExcitationHandler;
138	delete theH1Scatterer;
139	}
140
141	//----------------------------------------------------------------------------
142
143	//
144	// I M P L E M E N T A T I O N
145	//
146
147
148	//----------------------------------------------------------------------------
149	G4HadFinalState * G4BinaryCascade::ApplyYourself(const G4HadProjectile & aTrack,
150	//----------------------------------------------------------------------------
151	G4Nucleus & aNucleus)
152	{
153	static G4int eventcounter=0;
154	//if(eventcounter == 100*(eventcounter/100) )
155	eventcounter++;
156	if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl;
157
158	G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
159	G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition());
160
161	if(initial4Momentum.e()-initial4Momentum.m()<theBCminP &&
162	( definition==G4Neutron::NeutronDefinition() \|\| definition==G4Proton::ProtonDefinition() ) )
163	{
164	return theDeExcitation->ApplyYourself(aTrack, aNucleus);
165	}
166
167	theParticleChange.Clear();
168	// initialize the G4V3DNucleus from G4Nucleus
169	the3DNucleus = new G4Fancy3DNucleus;
170
171	// Build a KineticTrackVector with the G4Track
172	G4KineticTrackVector * secondaries;// = new G4KineticTrackVector;
173	G4ThreeVector initialPosition(0., 0., 0.); // will be set later
174
175	if(!getenv("I_Am_G4BinaryCascade_Developer") )
176	{
177	if(definition!=G4Neutron::NeutronDefinition() &&
178	definition!=G4Proton::ProtonDefinition() &&
179	definition!=G4PionPlus::PionPlusDefinition() &&
180	definition!=G4PionMinus::PionMinusDefinition() )
181	{
182	G4cerr << "You are using G4BinaryCascade for projectiles other than nucleons or pions."<<G4endl;
183	G4cerr << "If you want to continue, please switch on the developer environment: "<<G4endl;
184	G4cerr << "setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<G4endl;
185	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade - used for unvalid particle type - Fatal");
186	}
187	}
188
189	// keep primary
190	thePrimaryType = definition;
191	thePrimaryEscape = false;
192
193	// try until an interaction will happen
194	G4ReactionProductVector * products = 0;
195	G4int interactionCounter = 0;
196	do
197	{
198	// reset status that could be changed in previous loop event
199	theCollisionMgr->ClearAndDestroy();
200
201	if(products != 0)
202	{ // free memory from previous loop event
203	ClearAndDestroy(products);
204	delete products;
205	products=0;
206	}
207
208	the3DNucleus->Init(aNucleus.GetA_asInt(), aNucleus.GetZ_asInt());
209	thePropagator->Init(the3DNucleus);
210	// GF Leak on kt??? but where to delete?
211	G4KineticTrack * kt;// = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
212	do // sample impact parameter until collisions are found
213	{
214	theCurrentTime=0;
215	G4double radius = the3DNucleus->GetOuterRadius()+3*fermi;
216	initialPosition=GetSpherePoint(1.1*radius, initial4Momentum); // get random position
217	kt = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
218	kt->SetState(G4KineticTrack::outside);
219	// secondaries has been cleared by Propagate() in the previous loop event
220	secondaries= new G4KineticTrackVector;
221	secondaries->push_back(kt);
222	products = Propagate(secondaries, the3DNucleus);
223	} while(! products ); // until we FIND a collision...
224
225	if(++interactionCounter>99) break;
226	} while(products->size() == 0); // ...untill we find an ALLOWED collision
227
228	if(products->size()==0)
229	{
230	if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number void ######### "<<eventcounter<<G4endl;
231	theParticleChange.SetStatusChange(isAlive);
232	theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
233	theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
234	return &theParticleChange;
235	}
236	// G4cout << "BIC Applyyourself: number of products " << products->size() << G4endl;
237
238	// Fill the G4ParticleChange * with products
239	theParticleChange.SetStatusChange(stopAndKill);
240	G4ReactionProductVector::iterator iter;
241
242	for(iter = products->begin(); iter != products->end(); ++iter)
243	{
244	G4DynamicParticle * aNew =
245	new G4DynamicParticle((*iter)->GetDefinition(),
246	(*iter)->GetTotalEnergy(),
247	(*iter)->GetMomentum());
248	// FixMe: should I use "position" or "time" specifyed AddSecondary() methods?
249	theParticleChange.AddSecondary(aNew);
250
251	}
252
253	// DebugEpConservation(track, products);
254
255	ClearAndDestroy(products);
256	delete products;
257
258	delete the3DNucleus;
259	the3DNucleus = NULL; // protect from wrong usage...
260
261	if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number ends ######### "<<eventcounter<<G4endl;
262	if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
263	{
264	G4cout <<" BIC-fin-weight change " << theParticleChange.GetWeightChange()<< G4endl;
265	}
266	return &theParticleChange;
267	}
268
269	//----------------------------------------------------------------------------
270	G4ReactionProductVector * G4BinaryCascade::Propagate(
271	//----------------------------------------------------------------------------
272	G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus)
273	{
274	G4ping debug("debug_G4BinaryCascade");
275	debug.push_back("trial");
276	debug.dump();
277	#ifdef debug_BIC_Propagate
278	G4cout << "G4BinaryCascade Propagate starting -------------------------------------------------------" <<G4endl;
279	#endif
280	G4ReactionProductVector * products = new G4ReactionProductVector;
281	the3DNucleus = nucleus;
282	theOuterRadius = the3DNucleus->GetOuterRadius();
283	theCurrentTime=0;
284	// build theSecondaryList, theProjectileList and theCapturedList
285	ClearAndDestroy(&theCapturedList);
286	ClearAndDestroy(&theSecondaryList);
287	theSecondaryList.clear();
288	ClearAndDestroy(&theProjectileList);
289	ClearAndDestroy(&theFinalState);
290	std::vector<G4KineticTrack *>::iterator iter;
291
292	// GF FIXME ? in propagate mode this test is wrong! Could be in Apply....
293	if(nucleus->GetMassNumber() == 1) // 1H1 is special case
294	{
295	#ifdef debug_BIC_Propagate
296	G4cout << " special case 1H1.... " << G4endl;
297	#endif
298	return Propagate1H1(secondaries,nucleus);
299	}
300
301	BuildTargetList();
302
303	#ifdef debug_BIC_GetExcitationEnergy
304	G4cout << "ExcitationEnergy0 " << GetExcitationEnergy() << G4endl;
305	#endif
306
307	thePropagator->Init(the3DNucleus);
308
309
310	theCutOnP=90*MeV;
311	if(nucleus->GetMass()>30) theCutOnP = 70*MeV;
312	if(nucleus->GetMass()>60) theCutOnP = 50*MeV;
313	if(nucleus->GetMass()>120) theCutOnP = 45*MeV;
314
315	G4double StartingTime=DBL_MAX; // Search for minimal formation time Uzhi
316	for(iter = secondaries->begin(); iter != secondaries->end(); ++iter) // Uzhi
317	{ // Uzhi
318	if((*iter)->GetFormationTime() < StartingTime) // Uzhi
319	StartingTime = (*iter)->GetFormationTime(); // Uzhi
320	} // Uzhi
321
322	for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
323	{
324	// G4cout << " Formation time : " << (*iter)->GetDefinition()->GetParticleName() << " "
325	// << (*iter)->GetFormationTime() << G4endl;
326	G4double FormTime = (*iter)->GetFormationTime() - StartingTime; // Uzhi
327	(*iter)->SetFormationTime(FormTime); // Uzhi
328	if( (*iter)->GetState() == G4KineticTrack::undefined )
329	{
330	FindLateParticleCollision(*iter);
331	} else
332	{
333	theSecondaryList.push_back(*iter);
334	#ifdef debug_BIC_Propagate
335	G4cout << " Adding initial secondary " << *iter
336	<< " time" << (*iter)->GetFormationTime()
337	<< ", state " << (*iter)->GetState() << G4endl;
338	#endif
339	}
340	// theCollisionMgr->Print();
341	theProjectileList.push_back(new G4KineticTrack((iter)));
342	}
343	FindCollisions(&theSecondaryList);
344	secondaries->clear(); // Don't leave "G4KineticTrack *"s in two vectors
345	delete secondaries;
346
347	// if called stand alone, build theTargetList and find first collisions
348
349	if(theCollisionMgr->Entries() == 0 ) //late particles ALWAYS create Entries
350	{
351	//G4cout << " no collsions -> return 0" << G4endl;
352	delete products;
353	return 0;
354	}
355
356	// end of initialization: do the job now
357	// loop untill there are no more collisions
358
359
360	G4bool haveProducts = false;
361	G4int collisionCount=0;
362	theMomentumTransfer=G4ThreeVector(0,0,0);
363	while(theCollisionMgr->Entries() > 0)
364	{
365	// G4cout << "Propagate Captured size: " << theCapturedList.size() << G4endl;
366	// FixMe: at the moment I absorb pi with mom < theCutOnp, and I capture
367	// p and n with mom < theCutOnP. What about antip, k and sigmas
368	// with mom < theCutOnP?
369	if(Absorb()) { // absorb secondaries
370	haveProducts = true;
371	//G4cout << "Absorb sucess " << G4endl;
372	}
373
374	if(Capture()) { // capture secondaries (nucleons)
375	haveProducts = true;
376	//G4cout << "Capture sucess " << G4endl;
377	}
378
379	// propagate to the next collision if any (collisions could have been deleted
380	// by previous absorption or capture)
381	if(theCollisionMgr->Entries() > 0)
382	{
383	G4CollisionInitialState *
384	nextCollision = theCollisionMgr->GetNextCollision();
385	#ifdef debug_BIC_Propagate_Collisions
386	G4cout << " NextCollision * , Time, curtime = " << nextCollision << " "
387	<<nextCollision->GetCollisionTime()<< " " <<
388	theCurrentTime<< G4endl;
389	#endif
390	debug.push_back("======> test 1"); debug.dump();
391	if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) )
392	{
393	// Check if nextCollision is still valid, ie. particle did not leave nucleus
394	if (theCollisionMgr->GetNextCollision() != nextCollision )
395	{
396	nextCollision = 0;
397	}
398	}
399	debug.push_back("======> test 2"); debug.dump();
400	// G4cerr <<"post- DoTimeStep 1"<<G4endl;
401
402	if( nextCollision )
403	{
404	debug.push_back("======> test 3"); debug.dump();
405	if (ApplyCollision(nextCollision))
406	{
407	//G4cerr << "ApplyCollision sucess " << G4endl;
408	haveProducts = true;
409	collisionCount++;
410	debug.push_back("======> test 4.1"); debug.dump();
411	} else {
412	//G4cerr << "ApplyCollision failure " << G4endl;
413	theCollisionMgr->RemoveCollision(nextCollision);
414	debug.push_back("======> test 4.2"); debug.dump();
415	}
416	}
417	debug.push_back("======> test 5"); debug.dump();
418	// G4cerr <<"post-post- DoTimeStep 1"<<G4endl;
419	}
420	}
421
422	//--------- end of while on Collsions
423
424	// No more collisions: absorb, capture and propagate the secondaries out of the nucleus
425	if(Absorb()) {
426	haveProducts = true;
427	// G4cout << "Absorb sucess " << G4endl;
428	}
429
430	if(Capture()) {
431	haveProducts = true;
432	// G4cout << "Capture sucess " << G4endl;
433	}
434
435	if(!haveProducts) // no collisions happened. Return an empty vector.
436	{
437	//G4cout << " no products return 0" << G4endl;
438	return products;
439	}
440
441
442	#ifdef debug_BIC_Propagate
443	G4cout << " Momentum transfer to Nucleus " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
444	G4cout << " Stepping particles out...... " << G4endl;
445	#endif
446
447	StepParticlesOut();
448
449
450	if ( theSecondaryList.size() > 0 )
451	{
452	#ifdef debug_G4BinaryCascade
453	G4cerr << "G4BinaryCascade: Warning, have active particles at end" << G4endl;
454	#endif
455	// add left secondaries to FinalSate
456	std::vector<G4KineticTrack *>::iterator iter;
457	for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
458	{
459	theFinalState.push_back(*iter);
460	}
461	theSecondaryList.clear();
462
463	}
464	while ( theCollisionMgr->Entries() > 0 )
465	{
466	#ifdef debug_G4BinaryCascade
467	G4cerr << " Warning: remove left over collision " << G4endl;
468	#endif
469	theCollisionMgr->RemoveCollision(theCollisionMgr->GetNextCollision());
470	}
471
472	#ifdef debug_BIC_Propagate_Excitation
473
474	PrintKTVector(&theProjectileList,std::string(" theProjectileList"));
475	PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
476	G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
477	// PrintKTVector(&theTargetList,std::string(" theTargetList"));
478	PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
479
480	G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
481	G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
482	PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
483	#endif
484
485	//
486
487
488	G4double ExcitationEnergy=GetExcitationEnergy();
489
490	#ifdef debug_BIC_Propagate_finals
491	PrintKTVector(&theFinalState,std::string(" FinalState be4 corr"));
492	G4cout << " Excitation Energy prefinal, #collisions:, out, captured "
493	<< ExcitationEnergy << " "
494	<< collisionCount << " "
495	<< theFinalState.size() << " "
496	<< theCapturedList.size()<<G4endl;
497	#endif
498
499	if (ExcitationEnergy < 0 )
500	{
501	G4int maxtry=5, ntry=0;
502	do {
503	CorrectFinalPandE();
504	ExcitationEnergy=GetExcitationEnergy();
505	} while ( ++ntry < maxtry && ExcitationEnergy < 0 );
506	}
507
508	#ifdef debug_BIC_Propagate_finals
509	PrintKTVector(&theFinalState,std::string(" FinalState corrected"));
510	G4cout << " Excitation Energy final, #collisions:, out, captured "
511	<< ExcitationEnergy << " "
512	<< collisionCount << " "
513	<< theFinalState.size() << " "
514	<< theCapturedList.size()<<G4endl;
515	#endif
516
517
518	if ( ExcitationEnergy < 0. )
519	{
520	// if ( ExcitationEnergy < 0. )
521	{
522	//#ifdef debug_G4BinaryCascade
523	// G4cerr << "G4BinaryCascade-Warning: negative excitation energy ";
524	// G4cerr <<ExcitationEnergy<<G4endl;
525	// PrintKTVector(&theFinalState,std::string("FinalState"));
526	// PrintKTVector(&theCapturedList,std::string("captured"));
527	// G4cout << "negative ExE:Final 4Momentum .mag: " << GetFinal4Momentum()
528	// << " "<< GetFinal4Momentum().mag()<< G4endl
529	// << "negative ExE:FinalNucleusMom .mag: " << GetFinalNucleusMomentum()
530	// << " "<< GetFinalNucleusMomentum().mag()<< G4endl;
531	//#endif
532	}
533	ClearAndDestroy(products);
534	//G4cout << " negative Excitation E return empty products " << products << G4endl;
535	return products; // return empty products
536	}
537
538	// find a fragment and call the precompound model.
539	G4Fragment * fragment = 0;
540	G4ReactionProductVector * precompoundProducts = 0;
541
542	G4LorentzVector pFragment(0);
543	// G4cout << " final4mon " << GetFinal4Momentum() /MeV << G4endl;
544
545	// if ( ExcitationEnergy >= 0 ) // closed by Uzhi
546	// { // closed by Uzhi
547	fragment = FindFragments();
548	if(fragment) // Uzhi
549	{ // Uzhi
550	if(fragment->GetA() >1.5) // Uzhi
551	{
552	if (theDeExcitation) // pre-compound
553	{
554	// G4cout << " going to preco with fragment 4 mom " << fragment->GetMomentum() << G4endl;
555	pFragment=fragment->GetMomentum();
556	precompoundProducts= theDeExcitation->DeExcite(*fragment);
557	delete fragment;
558	fragment=0;
559	} else if (theExcitationHandler) // de-excitation
560	{
561	// G4cout << " going to de-excit with fragment 4 mom " << fragment->GetMomentum() << G4endl;
562	pFragment=fragment->GetMomentum();
563	precompoundProducts=theExcitationHandler->BreakItUp(*fragment);
564	delete fragment;
565	fragment=0;
566	}
567	} else
568	{ // fragment->GetA() < 1.5
569	precompoundProducts = new G4ReactionProductVector();
570	std::vector<G4KineticTrack *>::iterator i;
571	if ( theTargetList.size() == 1 )
572	{
573	i=theTargetList.begin();
574	G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition());
575	aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass());
576	aNew->SetMomentum(G4ThreeVector(0));// see boost for preCompoundProducts below..
577	precompoundProducts->push_back(aNew);
578	}
579
580	if ( theCapturedList.size() == 1 ) // Uzhi
581	{ // Uzhi
582	i=theCapturedList.begin(); // Uzhi
583	G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition()); // Uzhi
584	aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass()); // Uzhi
585	aNew->SetMomentum(G4ThreeVector(0));// see boost below.. // Uzhi
586	precompoundProducts->push_back(aNew); // Uzhi
587	} // Uzhi
588	} // End of fragment->GetA() < 1.5
589	} else // End of if(fragment)
590	{ // No fragment, can be neutrons only // Uzhi
591	precompoundProducts = new G4ReactionProductVector();
592
593	if ( (theTargetList.size()+theCapturedList.size()) > 0 )
594	{
595	std::vector<G4KineticTrack *>::iterator aNuc;
596	G4LorentzVector aVec;
597	std::vector<G4double> masses;
598	G4double sumMass(0);
599
600	if ( theTargetList.size() != 0) // Uzhi
601	{
602	for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
603	{
604	G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
605	masses.push_back(mass);
606	sumMass += mass;
607	}
608	} // Uzhi
609
610	if ( theCapturedList.size() != 0) // Uzhi
611	{ // Uzhi
612	for(aNuc = theCapturedList.begin(); // Uzhi
613	aNuc != theCapturedList.end(); aNuc++) // Uzhi
614	{ // Uzhi
615	G4double mass=(*aNuc)->GetDefinition()->GetPDGMass(); // Uzhi
616	masses.push_back(mass); // Uzhi
617	sumMass += mass; // Uzhi
618	}
619	}
620
621	G4LorentzVector finalP=GetFinal4Momentum();
622	G4FermiPhaseSpaceDecay decay;
623	// G4cout << " some neutrons? " << masses.size() <<" " ;
624	// G4cout<< theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl;
625
626	G4double eCMS=finalP.mag();
627	if ( eCMS < sumMass ) // @@GF --- Cheat!!
628	{
629	eCMS=sumMass + (2MeVmasses.size());
630	finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS)));
631	}
632
633	precompoundLorentzboost.set(finalP.boostVector());
634	std::vector<G4LorentzVector> momenta=decay.Decay(eCMS,masses);
635	std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
636
637	if ( theTargetList.size() != 0)
638	{
639	for ( aNuc=theTargetList.begin();
640	(aNuc != theTargetList.end()) && (aMom!=momenta->end());
641	aNuc++, aMom++ )
642	{
643	G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition());
644	aNew->SetTotalEnergy((*aMom)->e());
645	aNew->SetMomentum((*aMom)->vect());
646	precompoundProducts->push_back(aNew);
647
648	delete *aMom;
649	}
650	}
651
652	if ( theCapturedList.size() != 0) // Uzhi
653	{ // Uzhi
654	for ( aNuc=theCapturedList.begin(); // Uzhi
655	(aNuc != theCapturedList.end()) && (aMom!=momenta->end()); // Uzhi
656	aNuc++, aMom++ ) // Uzhi
657	{ // Uzhi
658	G4ReactionProduct * aNew = new G4ReactionProduct( // Uzhi
659	(*aNuc)->GetDefinition()); // Uzhi
660	aNew->SetTotalEnergy((*aMom)->e()); // Uzhi
661	aNew->SetMomentum((*aMom)->vect()); // Uzhi
662	precompoundProducts->push_back(aNew); // Uzhi
663	delete *aMom; // Uzhi
664	} // Uzhi
665	} // Uzhi
666
667	if(momenta) delete momenta;
668	}
669	} // End if(!fragment)
670
671
672	{
673	// fill in products the outgoing particles
674	G4double Ekinout=0;
675	G4LorentzVector pSumBic(0);
676	size_t i(0);
677	for(i = 0; i< theFinalState.size(); i++)
678	{
679	G4KineticTrack * kt = theFinalState[i];
680	if(kt->GetDefinition()->IsShortLived())
681	{
682	G4KineticTrackVector * dec = kt->Decay();
683	G4KineticTrackVector::const_iterator jter;
684	for(jter=dec->begin(); jter != dec->end(); jter++)
685	{
686	theFinalState.push_back(*jter);
687	}
688	delete dec;
689	}
690	else
691	{
692	G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
693	aNew->SetMomentum(kt->Get4Momentum().vect());
694	aNew->SetTotalEnergy(kt->Get4Momentum().e());
695	Ekinout += aNew->GetKineticEnergy();
696	pSumBic += kt->Get4Momentum();
697	if(kt->IsParticipant())
698	{
699	aNew->SetNewlyAdded(true);
700	}
701	else
702	{
703	aNew->SetNewlyAdded(false);
704	}
705	//G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
706	products->push_back(aNew);
707
708	#ifdef debug_BIC_Propagate_finals
709	if (! kt->GetDefinition()->GetPDGStable() )
710	{
711	if (kt->GetDefinition()->IsShortLived())
712	{
713	G4cout << "final shortlived : ";
714	} else
715	{
716	G4cout << "final non stable : ";
717	}
718	G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
719	}
720	#endif
721	}
722
723	}
724	//G4cout << " Total Ekin " << Ekinout << G4endl;
725	}
726	// add precompound products to products
727	G4LorentzVector pSumPreco(0), pPreco(0);
728	if ( precompoundProducts )
729	{
730	std::vector<G4ReactionProduct *>::iterator j;
731	for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
732	{
733	// boost back to system of moving nucleus
734	G4LorentzVector pProduct((j)->GetMomentum(),(j)->GetTotalEnergy());
735	pPreco+= pProduct;
736	#ifdef debug_BIC_Propagate_finals
737	G4cout << " pProduct be4 boost " <<pProduct << G4endl;
738	#endif
739	pProduct *= precompoundLorentzboost;
740	#ifdef debug_BIC_Propagate_finals
741	G4cout << " pProduct aft boost " <<pProduct << G4endl;
742	#endif
743	pSumPreco += pProduct;
744	(*j)->SetTotalEnergy(pProduct.e());
745	(*j)->SetMomentum(pProduct.vect());
746	(*j)->SetNewlyAdded(true);
747	products->push_back(*j);
748	}
749	// G4cout << " unboosted preco result mom " << pPreco / MeV << " ..- fragmentMom " << (pPreco - pFragment)/MeV<< G4endl;
750	// G4cout << " preco result mom " << pSumPreco / MeV << " ..-file4Mom " << (pSumPreco - GetFinal4Momentum())/MeV<< G4endl;
751	precompoundProducts->clear();
752	delete precompoundProducts;
753	}
754	else
755	{
756	G4ReactionProduct * aNew=0;
757	// return nucleus e and p
758	if (fragment != 0 ) {
759	aNew = new G4ReactionProduct(G4Gamma::GammaDefinition()); // we only want to pass e/p
760	aNew->SetMomentum(fragment->GetMomentum().vect());
761	aNew->SetTotalEnergy(fragment->GetMomentum().e());
762	delete fragment;
763	fragment=0;
764	} else if (products->size() == 0) {
765	// FixMe GF: for testing without precompound, return 1 gamma of 0.01 MeV in +x
766	#include "G4Gamma.hh"
767	aNew = new G4ReactionProduct(G4Gamma::GammaDefinition());
768	aNew->SetMomentum(G4ThreeVector(0.01*MeV,0,0));
769	aNew->SetTotalEnergy(0.01*MeV);
770	}
771	if ( aNew != 0 ) products->push_back(aNew);
772	}
773
774	// G4cerr <<"mon - end of event "<<G4endl;
775	thePrimaryEscape = true;
776	//G4cout << " return products " << products << G4endl;
777	return products;
778	}
779
780
781	//----------------------------------------------------------------------------
782	G4double G4BinaryCascade::GetExcitationEnergy()
783	//----------------------------------------------------------------------------
784	{
785
786	G4ping debug("debug_ExcitationEnergy");
787	// get A and Z for the residual nucleus
788	#if defined(debug_G4BinaryCascade) \|\| defined(debug_BIC_GetExcitationEnergy)
789	G4int finalA = theTargetList.size()+theCapturedList.size();
790	G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
791	if ( (currentA - finalA) != 0 \|\| (currentZ - finalZ) != 0 )
792	{
793	G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} "
794	<< currentA << " " << finalA << " "<< currentZ << " " << finalZ << G4endl;
795	}
796
797	#endif
798
799	G4double excitationE(0);
800	G4double nucleusMass(0);
801	if(currentZ>.5)
802	{
803	nucleusMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA);
804	}
805	else if (currentZ==0 ) // Uzhi && currentA==1 ) // Uzhi
806	{ // Uzhi
807	if(currentA == 1) {nucleusMass = G4Neutron::Neutron()->GetPDGMass();}// Uzhi
808	else {nucleusMass = GetFinalNucleusMomentum().mag() // Uzhi
809	- 3.MeVcurrentA;} // Uzhi
810	} // Uzhi
811	else
812	{
813	#ifdef debug_G4BinaryCascade
814	G4cout << "G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=("
815	<< currentA << "," << currentZ << ")" << G4endl;
816	#endif
817	return 0;
818	}
819
820	#ifdef debug_BIC_GetExcitationEnergy
821	debug.push_back("====> current A, Z");
822	debug.push_back(currentZ);
823	debug.push_back(currentA);
824	debug.push_back(finalZ);
825	debug.push_back(finalA);
826	debug.push_back(nucleusMass);
827	debug.push_back(GetFinalNucleusMomentum().mag());
828	debug.dump();
829	// PrintKTVector(&theTargetList, std::string(" current target list info"));
830	PrintKTVector(&theCapturedList, std::string(" current captured list info"));
831	#endif
832
833	excitationE = GetFinalNucleusMomentum().mag() - nucleusMass;
834
835	#ifdef debug_BIC_GetExcitationEnergy
836	// ------ debug
837	if ( excitationE < 0 )
838	{
839	G4cout << "negative ExE final Ion mass " <<nucleusMass<< G4endl;
840	G4LorentzVector Nucl_mom=GetFinalNucleusMomentum();
841	if(finalZ>.5) G4cout << " Final nuclmom/mass " << Nucl_mom << " " << Nucl_mom.mag()
842	<< " (A,Z)=("<< finalA <<","<<finalZ <<")"
843	<< " mass " << nucleusMass << " "
844	<< " excitE " << excitationE << G4endl;
845
846
847	G4int A = the3DNucleus->GetMassNumber();
848	G4int Z = the3DNucleus->GetCharge();
849	G4double initialExc(0);
850	if(Z>.5)
851	{
852	initialExc = theInitial4Mom.mag()-
853	G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z, A);
854	G4cout << "GetExcitationEnergy: Initial nucleus A Z " << A << " " << Z << " " << initialExc << G4endl;
855	}
856	}
857
858	#endif
859
860	return excitationE;
861	}
862
863
864	//----------------------------------------------------------------------------
865	//
866	// P R I V A T E M E M B E R F U N C T I O N S
867	//
868	//----------------------------------------------------------------------------
869
870	//----------------------------------------------------------------------------
871	void G4BinaryCascade::BuildTargetList()
872	//----------------------------------------------------------------------------
873	{
874
875	if(!the3DNucleus->StartLoop())
876	{
877	// G4cerr << "G4BinaryCascade::BuildTargetList(): StartLoop() error!"
878	// << G4endl;
879	return;
880	}
881
882	ClearAndDestroy(&theTargetList); // clear theTargetList before rebuilding
883
884	G4Nucleon * nucleon;
885	G4ParticleDefinition * definition;
886	G4ThreeVector pos;
887	G4LorentzVector mom;
888	// if there are nucleon hit by higher energy models, then SUM(momenta) != 0
889	theInitial4Mom = G4LorentzVector(0,0,0,0);
890	currentA=0;
891	currentZ=0;
892	// G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
893	while((nucleon = the3DNucleus->GetNextNucleon()) != NULL)
894	{
895	// check if nucleon is hit by higher energy model.
896	if ( ! nucleon->AreYouHit() )
897	{
898	definition = nucleon->GetDefinition();
899	pos = nucleon->GetPosition();
900	mom = nucleon->GetMomentum();
901	// G4cout << "Nucleus " << pos.mag()/fermi << " " << mom.e() << G4endl;
902	theInitial4Mom += mom;
903	// the potential inside the nucleus is taken into account, and nucleons are on mass shell.
904	mom.setE( std::sqrt( mom.vect().mag2() + sqr(definition->GetPDGMass()) ) );
905	G4KineticTrack * kt = new G4KineticTrack(definition, 0., pos, mom);
906	kt->SetState(G4KineticTrack::inside);
907	kt->SetNucleon(nucleon);
908	theTargetList.push_back(kt);
909	++currentA;
910	if (definition->GetPDGCharge() > .5 ) ++currentZ;
911	}
912	#ifdef debug_BIC_BuildTargetList
913	else { G4cout << "nucleon is hit" << nucleon << G4endl;}
914	#endif
915	}
916	massInNucleus = 0;
917	if(currentZ>.5)
918	{
919	massInNucleus = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA);
920	} else if (currentZ==0 && currentA>=1 )
921	{
922	massInNucleus = currentA * G4Neutron::Neutron()->GetPDGMass();
923	} else
924	{
925	G4cerr << "G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=("
926	<< currentA << "," << currentZ << ")" << G4endl;
927	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::BuildTargetList()");
928	}
929	currentInitialEnergy= theInitial4Mom.e();
930	G4KineticTrackVector::iterator i;
931	for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
932	{
933	currentInitialEnergy+= (*i)->GetTrackingMomentum().e();
934	}
935
936	#ifdef debug_BIC_BuildTargetList
937	G4cout << "G4BinaryCascade::BuildTargetList(): nucleus (A,Z)=("
938	<< currentA << "," << currentZ << ") mass: " << massInNucleus <<
939	", theInitial4Mom " << theInitial4Mom <<
940	", currentInitialEnergy " << currentInitialEnergy << G4endl;
941	#endif
942
943	}
944
945
946	//----------------------------------------------------------------------------
947	void G4BinaryCascade::FindCollisions(G4KineticTrackVector * secondaries)
948	//----------------------------------------------------------------------------
949	{
950	for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin();
951	i != secondaries->end(); ++i)
952	{
953	#ifdef debug_G4BinaryCascade
954	if ( (i)->GetTrackingMomentum().mag2() < -1.eV )
955	{
956	G4cout << "G4BinaryCascade::FindCollisions(): negative m2:" << (*i)->GetTrackingMomentum().mag2() << G4endl;
957	}
958	#endif
959
960	for(std::vector<G4BCAction *>::iterator j = theImR.begin();
961	j!=theImR.end(); j++)
962	{
963	// G4cout << "G4BinaryCascade::FindCollisions: at action " << *j << G4endl;
964	const std::vector<G4CollisionInitialState *> & aCandList
965	= (j)->GetCollisions(i, theTargetList, theCurrentTime);
966	for(size_t count=0; count<aCandList.size(); count++)
967	{
968	theCollisionMgr->AddCollision(aCandList[count]);
969	// G4cout << "====================== New Collision ================="<<G4endl;
970	// theCollisionMgr->Print();
971	}
972	}
973	}
974	}
975
976	//----------------------------------------------------------------------------
977	void G4BinaryCascade::FindDecayCollision(G4KineticTrack * secondary)
978	//----------------------------------------------------------------------------
979	{
980	if ( secondary->GetTrackingMomentum().mag2() < -1.*eV )
981	{
982	G4cout << "G4BinaryCascade::FindDecayCollision(): negative m2:" << secondary->GetTrackingMomentum().mag2() << G4endl;
983	}
984	const std::vector<G4CollisionInitialState *> & aCandList
985	= theDecay->GetCollisions(secondary, theTargetList, theCurrentTime);
986	for(size_t count=0; count<aCandList.size(); count++)
987	{
988	theCollisionMgr->AddCollision(aCandList[count]);
989	}
990	}
991
992	//----------------------------------------------------------------------------
993	void G4BinaryCascade::FindLateParticleCollision(G4KineticTrack * secondary)
994	//----------------------------------------------------------------------------
995	{
996	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
997	if ( secondary->GetTrackingMomentum().mag2() < -1.*eV )
998	{
999	G4cout << "G4BinaryCascade::FindDecayCollision(): negative m2:" << secondary->GetTrackingMomentum().mag2() << G4endl;
1000	}
1001
1002	G4double tin=0., tout=0.;
1003	if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(secondary,tin,tout))
1004	{
1005	if ( tin > 0 )
1006	{
1007	secondary->SetState(G4KineticTrack::outside);
1008	} else if ( tout > 0 )
1009	{
1010	secondary->SetState(G4KineticTrack::inside);
1011	} else {
1012	//G4cout << "G4BC set miss , tin, tout " << tin << " , " << tout <<G4endl;
1013	secondary->SetState(G4KineticTrack::miss_nucleus);
1014	}
1015	} else {
1016	secondary->SetState(G4KineticTrack::miss_nucleus);
1017	//G4cout << "G4BC set miss ,no intersect tin, tout " << tin << " , " << tout <<G4endl;
1018	}
1019
1020	// for barions, correct for fermi energy
1021	G4int PDGcode=std::abs(secondary->GetDefinition()->GetPDGEncoding());
1022	if ( PDGcode > 1000 )
1023	{
1024	G4double initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
1025	secondary->GetPosition());
1026	secondary->Update4Momentum(secondary->Get4Momentum().e() - initial_Efermi);
1027	}
1028
1029	#ifdef debug_BIC_FindCollision
1030	G4cout << "FindLateP Particle, 4-mom, times newState "
1031	<< secondary->GetDefinition()->GetParticleName() << " "
1032	<< secondary->Get4Momentum()
1033	<< " times " << tin << " " << tout << " "
1034	<< secondary->GetState() << G4endl;
1035	#endif
1036
1037	const std::vector<G4CollisionInitialState *> & aCandList
1038	= theLateParticle->GetCollisions(secondary, theTargetList, theCurrentTime);
1039	for(size_t count=0; count<aCandList.size(); count++)
1040	{
1041	#ifdef debug_BIC_FindCollision
1042	G4cout << " Adding a late Col : " << aCandList[count] << G4endl;
1043	#endif
1044	theCollisionMgr->AddCollision(aCandList[count]);
1045	}
1046	}
1047
1048
1049	//----------------------------------------------------------------------------
1050	G4bool G4BinaryCascade::ApplyCollision(G4CollisionInitialState * collision)
1051	//----------------------------------------------------------------------------
1052	{
1053	G4ping debug("debug_ApplyCollision");
1054	#ifdef debug_BIC_ApplyCollision
1055	G4cerr << "G4BinaryCascade::ApplyCollision start"<<G4endl;
1056	theCollisionMgr->Print();
1057	#endif
1058	G4KineticTrack * primary = collision->GetPrimary();
1059	G4KineticTrackVector target_collection=collision->GetTargetCollection();
1060	G4bool haveTarget=target_collection.size()>0;
1061	if( haveTarget && (primary->GetState() != G4KineticTrack::inside) )
1062	{
1063	#ifdef debug_G4BinaryCascade
1064	G4cout << "G4BinaryCasacde::ApplyCollision(): StateError " << primary << G4endl;
1065	G4KineticTrackVector debug;
1066	debug.push_back(primary);
1067	PrintKTVector(&debug,std::string("primay- ..."));
1068	PrintKTVector(&target_collection,std::string("... targets"));
1069	//GF throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
1070	#else
1071	return false;
1072	#endif
1073	}
1074
1075	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
1076
1077	#ifdef debug_BIC_ApplyCollision
1078	G4cout << "ApplyCollisions : projte 4mom " << primary->GetTrackingMomentum()<< G4endl;
1079	#endif
1080
1081	G4int initialBaryon(0);
1082	G4int initialCharge(0);
1083	G4LorentzVector mom4Primary(0);
1084
1085	if (primary->GetState() == G4KineticTrack::inside)
1086	{
1087	initialBaryon = primary->GetDefinition()->GetBaryonNumber();
1088	initialCharge = G4lrint(primary->GetDefinition()->GetPDGCharge());
1089	}
1090
1091	G4int PDGcode=std::abs(primary->GetDefinition()->GetPDGEncoding());
1092	mom4Primary=primary->Get4Momentum();
1093
1094	// for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
1095	G4double initial_Efermi(0);
1096	if (primary->GetState() == G4KineticTrack::inside ) {
1097	initial_Efermi=RKprop->GetField(primary->GetDefinition()->GetPDGEncoding(),primary->GetPosition());
1098
1099	if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1100	{
1101	initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
1102	primary->GetPosition());
1103	primary->Update4Momentum(mom4Primary.e() - initial_Efermi);
1104	}
1105
1106	std::vector<G4KineticTrack *>::iterator titer;
1107	for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
1108	{
1109	G4ParticleDefinition * aDef=(*titer)->GetDefinition();
1110	G4int aCode=aDef->GetPDGEncoding();
1111	G4ThreeVector aPos=(*titer)->GetPosition();
1112	initial_Efermi+= RKprop->GetField(aCode, aPos);
1113	}
1114	}
1115	//****************************************
1116
1117	G4KineticTrackVector * products=0;
1118	products = collision->GetFinalState();
1119
1120	#ifdef debug_BIC_ApplyCollision
1121	G4bool havePion=false;
1122	for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1123	{
1124	G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
1125	if (std::abs(PDGcode)==211 \|\| PDGcode==111 ) havePion=true;
1126	}
1127	if ( !products \|\| havePion)
1128	{
1129	G4cout << " Collision " << collision << ", type: "<< typeid(*collision->GetGenerator()).name()
1130	<< ", with NO products! " <<G4endl;
1131	G4cout << G4endl<<"Initial condition are these:"<<G4endl;
1132	G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
1133	PrintKTVector(collision->GetPrimary());
1134	for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
1135	{
1136	G4cout << "targ: "
1137	<<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
1138	}
1139	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
1140	}
1141	#endif
1142	G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
1143	// if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl;
1144	// if ( lateParticleCollision && products ) PrintKTVector(products, " reaction products");
1145	//****************************************
1146
1147	// reset primary to initial state
1148	primary->Set4Momentum(mom4Primary);
1149
1150
1151	G4int lateBaryon(0), lateCharge(0);
1152
1153	if ( lateParticleCollision )
1154	{ // for late particles, reset charges
1155	//G4cout << "lateP, initial B C state " << initialBaryon << " "
1156	// << initialCharge<< " " << primary->GetState() << " "<< primary->GetDefinition()->GetParticleName()<< G4endl;
1157	lateBaryon = initialBaryon;
1158	lateCharge = initialCharge;
1159	initialBaryon=initialCharge=0;
1160	}
1161
1162	initialBaryon += collision->GetTargetBaryonNumber();
1163	initialCharge+=G4lrint(collision->GetTargetCharge());
1164	if(!products \|\| products->size()==0 \|\| !CheckPauliPrinciple(products))
1165	{
1166	#ifdef debug_BIC_ApplyCollision
1167	if (products) G4cout << " ======Failed Pauli =====" << G4endl;
1168	G4cerr << "G4BinaryCascade::ApplyCollision blocked"<<G4endl;
1169	#endif
1170	if (products) ClearAndDestroy(products);
1171	if ( ! haveTarget ) FindDecayCollision(primary); // for decay, sample new decay
1172	delete products;
1173	return false;
1174	}
1175
1176	if (primary->GetState() == G4KineticTrack::inside ) { // if the primary was outside, nothing to correct
1177	G4double final_Efermi(0);
1178	G4KineticTrackVector resonances;
1179	for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1180	{
1181	G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
1182	// G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
1183	final_Efermi+=RKprop->GetField(PDGcode,(*i)->GetPosition());
1184	if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1185	{
1186	resonances.push_back(*i);
1187	}
1188	}
1189	if ( resonances.size() > 0 )
1190	{
1191	G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
1192	for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
1193	{
1194	G4LorentzVector mom=(*res)->Get4Momentum();
1195	G4double mass2=mom.mag2();
1196	G4double newEnergy=mom.e() + delta_Fermi;
1197	G4double newEnergy2= newEnergy*newEnergy;
1198	//G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
1199	if ( newEnergy2 < mass2 )
1200	{
1201	ClearAndDestroy(products);
1202	if (target_collection.size() == 0 ) FindDecayCollision(primary); // for decay, sample new decay
1203	delete products;
1204	return false;
1205	}
1206	// G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl;
1207	G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
1208	(*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
1209	}
1210	}
1211	}
1212
1213	#ifdef debug_BIC_ApplyCollision
1214	DebugApplyCollision(collision, products);
1215	#endif
1216
1217	G4int finalBaryon(0);
1218	G4int finalCharge(0);
1219	G4KineticTrackVector toFinalState;
1220	for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1221	{
1222	if ( ! lateParticleCollision )
1223	{
1224	(*i)->SetState(primary->GetState()); // decay may be anywhere!
1225	if ( (*i)->GetState() == G4KineticTrack::inside ){
1226	finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1227	finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
1228	}
1229	} else {
1230	G4double tin=0., tout=0.;
1231	if (((G4RKPropagation)thePropagator)->GetSphereIntersectionTimes((i),tin,tout))
1232	{
1233	if ( tin > 0 )
1234	{
1235	(*i)->SetState(G4KineticTrack::outside);
1236	}
1237	else if ( tout > 0 )
1238	{
1239	(*i)->SetState(G4KineticTrack::inside);
1240	finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1241	finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
1242	}
1243	else
1244	{
1245	(*i)->SetState(G4KineticTrack::gone_out);
1246	toFinalState.push_back((*i));
1247	}
1248	} else
1249	{
1250	(*i)->SetState(G4KineticTrack::miss_nucleus);
1251	//G4cout << " G4BC - miss -late Part- no intersection found " << G4endl;
1252	toFinalState.push_back((*i));
1253	}
1254
1255	//G4cout << " PDGcode, state " << (i)->GetDefinition()->GetPDGEncoding() << " " << (i)->GetState()<<G4endl;
1256
1257	}
1258	}
1259	if(!toFinalState.empty())
1260	{
1261	theFinalState.insert(theFinalState.end(),
1262	toFinalState.begin(),toFinalState.end());
1263	std::vector<G4KineticTrack *>::iterator iter1, iter2;
1264	for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
1265	++iter1)
1266	{
1267	iter2 = std::find(products->begin(), products->end(),
1268	*iter1);
1269	if ( iter2 != products->end() ) products->erase(iter2);
1270	}
1271	theCollisionMgr->RemoveTracksCollisions(&toFinalState);
1272	}
1273
1274	//G4cout << " currentA, Z be4: " << currentA << " " << currentZ << G4endl;
1275	currentA += finalBaryon-initialBaryon;
1276	currentZ += finalCharge-initialCharge;
1277	//G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl;
1278
1279	G4KineticTrackVector oldSecondaries;
1280	if (primary)
1281	{
1282	oldSecondaries.push_back(primary);
1283	primary->Hit();
1284	}
1285
1286	#ifdef debug_G4BinaryCascade
1287	if ( (finalBaryon-initialBaryon-lateBaryon) != 0 \|\| (finalCharge-initialCharge-lateCharge) != 0 )
1288	{
1289	G4cout << "G4BinaryCascade: Error in Balancing: " << G4endl;
1290	G4cout << "initial/final baryon number, initial/final Charge "
1291	<< initialBaryon <<" "<< finalBaryon <<" "
1292	<< initialCharge <<" "<< finalCharge <<" "
1293	<< " in Collision type: "<< typeid(*collision->GetGenerator()).name()
1294	<< ", with number of products: "<< products->size() <<G4endl;
1295	G4cout << G4endl<<"Initial condition are these:"<<G4endl;
1296	G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
1297	for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
1298	{
1299	G4cout << "targ: "
1300	<<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
1301	}
1302	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
1303	G4cout << G4endl<<G4endl;
1304	}
1305	#endif
1306
1307	G4KineticTrackVector oldTarget = collision->GetTargetCollection();
1308	for(size_t ii=0; ii< oldTarget.size(); ii++)
1309	{
1310	oldTarget[ii]->Hit();
1311	}
1312
1313	debug.push_back("=======> we have hit nucleons <=======");
1314
1315	UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products);
1316	std::for_each(oldSecondaries.begin(), oldSecondaries.end(), Delete<G4KineticTrack>());
1317	std::for_each(oldTarget.begin(), oldTarget.end(), Delete<G4KineticTrack>());
1318
1319	delete products;
1320	return true;
1321	}
1322
1323
1324	//----------------------------------------------------------------------------
1325	G4bool G4BinaryCascade::Absorb()
1326	//----------------------------------------------------------------------------
1327	{
1328	// Do it in two step: cannot change theSecondaryList inside the first loop.
1329	G4Absorber absorber(theCutOnPAbsorb);
1330
1331	// Build the vector of G4KineticTracks that must be absorbed
1332	G4KineticTrackVector absorbList;
1333	std::vector<G4KineticTrack *>::iterator iter;
1334	// PrintKTVector(&theSecondaryList, " testing for Absorb" );
1335	for(iter = theSecondaryList.begin();
1336	iter != theSecondaryList.end(); ++iter)
1337	{
1338	G4KineticTrack * kt = *iter;
1339	if(kt->GetState() == G4KineticTrack::inside)// absorption happens only inside the nucleus
1340	{
1341	if(absorber.WillBeAbsorbed(*kt))
1342	{
1343	absorbList.push_back(kt);
1344	}
1345	}
1346	}
1347
1348	if(absorbList.empty())
1349	return false;
1350
1351	G4KineticTrackVector toDelete;
1352	for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
1353	{
1354	G4KineticTrack * kt = *iter;
1355	if(!absorber.FindAbsorbers(*kt, theTargetList))
1356	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1357
1358	if(!absorber.FindProducts(*kt))
1359	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1360
1361	G4KineticTrackVector * products = absorber.GetProducts();
1362	// ------ debug
1363	G4int count = 0;
1364	// ------ end debug
1365	while(!CheckPauliPrinciple(products))
1366	{
1367	// ------ debug
1368	count++;
1369	// ------ end debug
1370	ClearAndDestroy(products);
1371	if(!absorber.FindProducts(*kt))
1372	throw G4HadronicException(__FILE__, __LINE__,
1373	"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1374	}
1375	// ------ debug
1376	// G4cerr << "Absorb CheckPauliPrinciple count= " << count << G4endl;
1377	// ------ end debug
1378	G4KineticTrackVector toRemove; // build a vector for UpdateTrack...
1379	toRemove.push_back(kt);
1380	toDelete.push_back(kt); // delete the track later
1381	G4KineticTrackVector * absorbers = absorber.GetAbsorbers();
1382	UpdateTracksAndCollisions(&toRemove, absorbers, products);
1383	ClearAndDestroy(absorbers);
1384	}
1385	ClearAndDestroy(&toDelete);
1386	return true;
1387	}
1388
1389
1390
1391	// Capture all p and n with Energy < theCutOnP
1392	//----------------------------------------------------------------------------
1393	G4bool G4BinaryCascade::Capture(G4bool verbose)
1394	//----------------------------------------------------------------------------
1395	{
1396	G4KineticTrackVector captured;
1397	G4bool capture = false;
1398	std::vector<G4KineticTrack *>::iterator i;
1399
1400	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
1401
1402	G4double capturedEnergy = 0;
1403	G4int particlesAboveCut=0;
1404	G4int particlesBelowCut=0;
1405	if ( verbose ) G4cout << " Capture: secondaries " << theSecondaryList.size() << G4endl;
1406	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1407	{
1408	G4KineticTrack * kt = *i;
1409	if (verbose) G4cout << "Capture position, radius, state " <<kt->GetPosition().mag()<<" "<<theOuterRadius<<" "<<kt->GetState()<<G4endl;
1410	if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1411	{
1412	if((kt->GetDefinition() == G4Proton::Proton()) \|\|
1413	(kt->GetDefinition() == G4Neutron::Neutron()))
1414	{
1415	//GF cut on kinetic energy if(kt->Get4Momentum().vect().mag() >= theCutOnP)
1416	G4double field=RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
1417	-RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
1418	G4double energy= kt->Get4Momentum().e() - kt->GetActualMass() + field;
1419	if (verbose ) G4cout << "Capture: .e(), mass, field, energy" << kt->Get4Momentum().e() <<" "<<kt->GetActualMass()<<" "<<field<<" "<<energy<< G4endl;
1420	// if( energy < theCutOnP )
1421	// {
1422	capturedEnergy+=energy;
1423	++particlesBelowCut;
1424	// } else
1425	// {
1426	// ++particlesAboveCut;
1427	// }
1428	}
1429	}
1430	}
1431	if (verbose) G4cout << "Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
1432	<< particlesAboveCut << " " << particlesBelowCut << " " << capturedEnergy
1433	<< " " << capturedEnergy/particlesBelowCut << " " << 0.2*theCutOnP << G4endl;
1434	// if(particlesAboveCut==0 && particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1435	if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1436	{
1437	capture=true;
1438	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1439	{
1440	G4KineticTrack * kt = *i;
1441	if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1442	{
1443	if((kt->GetDefinition() == G4Proton::Proton()) \|\|
1444	(kt->GetDefinition() == G4Neutron::Neutron()))
1445	{
1446	captured.push_back(kt);
1447	kt->Hit(); //
1448	theCapturedList.push_back(kt);
1449	}
1450	}
1451	}
1452	UpdateTracksAndCollisions(&captured, NULL, NULL);
1453	}
1454
1455	return capture;
1456	}
1457
1458
1459	//----------------------------------------------------------------------------
1460	G4bool G4BinaryCascade::CheckPauliPrinciple(G4KineticTrackVector * products)
1461	//----------------------------------------------------------------------------
1462	{
1463	G4int A = the3DNucleus->GetMassNumber();
1464	G4int Z = the3DNucleus->GetCharge();
1465
1466	G4FermiMomentum fermiMom;
1467	fermiMom.Init(A, Z);
1468
1469	const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
1470
1471	G4KineticTrackVector::iterator i;
1472	G4ParticleDefinition * definition;
1473
1474	// ------ debug
1475	G4bool myflag = true;
1476	// ------ end debug
1477	// G4ThreeVector xpos(0);
1478	for(i = products->begin(); i != products->end(); ++i)
1479	{
1480	definition = (*i)->GetDefinition();
1481	if((definition == G4Proton::Proton()) \|\|
1482	(definition == G4Neutron::Neutron()))
1483	{
1484	G4ThreeVector pos = (*i)->GetPosition();
1485	G4double d = density->GetDensity(pos);
1486	// energy correspondiing to fermi momentum
1487	G4double eFermi = std::sqrt( sqr(fermiMom.GetFermiMomentum(d)) + (*i)->Get4Momentum().mag2() );
1488	if( definition == G4Proton::Proton() )
1489	{
1490	eFermi -= the3DNucleus->CoulombBarrier();
1491	}
1492	G4LorentzVector mom = (*i)->Get4Momentum();
1493	// ------ debug
1494	/*
1495	* G4cout << "p:[" << (1/MeV)mom.x() << " " << (1/MeV)mom.y() << " "
1496	* << (1/MeV)*mom.z() << "] \|p3\|:"
1497	* << (1/MeV)mom.vect().mag() << " E:" << (1/MeV)mom.t() << " ] m: "
1498	* << (1/MeV)*mom.mag() << " pos["
1499	* << (1/fermi)pos.x() << " "<< (1/fermi)pos.y() << " "
1500	* << (1/fermi)*pos.z() << "] \|Dpos\|: "
1501	* << (1/fermi)*(pos-xpos).mag() << " Pfermi:"
1502	* << (1/MeV)*p << G4endl;
1503	* xpos=pos;
1504	*/
1505	// ------ end debug
1506	if(mom.e() < eFermi )
1507	{
1508	// ------ debug
1509	myflag = false;
1510	// ------ end debug
1511	// return false;
1512	}
1513	}
1514	}
1515	#ifdef debug_BIC_CheckPauli
1516	if ( myflag )
1517	{
1518	for(i = products->begin(); i != products->end(); ++i)
1519	{
1520	definition = (*i)->GetDefinition();
1521	if((definition == G4Proton::Proton()) \|\|
1522	(definition == G4Neutron::Neutron()))
1523	{
1524	G4ThreeVector pos = (*i)->GetPosition();
1525	G4double d = density->GetDensity(pos);
1526	G4double pFermi = fermiMom.GetFermiMomentum(d);
1527	G4LorentzVector mom = (*i)->Get4Momentum();
1528	G4double field =((G4RKPropagation*)thePropagator)->GetField(definition->GetPDGEncoding(),pos);
1529	if ( mom.e()-mom.mag()+field > 160*MeV )
1530	{
1531	G4cout << "momentum problem pFermi=" << pFermi
1532	<< " mom, mom.m " << mom << " " << mom.mag()
1533	<< " field " << field << G4endl;
1534	}
1535	}
1536	}
1537	}
1538	#endif
1539
1540	return myflag;
1541	}
1542
1543	//----------------------------------------------------------------------------
1544	void G4BinaryCascade::StepParticlesOut()
1545	//----------------------------------------------------------------------------
1546	{
1547	G4int counter=0;
1548	G4int countreset=0;
1549	//G4cout << " nucl. Radius " << radius << G4endl;
1550	// G4cerr <<"pre-while- theSecondaryList "<<G4endl;
1551	while( theSecondaryList.size() > 0 )
1552	{
1553	G4int nsec=0;
1554	G4double minTimeStep = 1.e-12ns; // about 30fermi/(0.1c_light);1.e-12ns
1555	// i.e. a big step
1556	std::vector<G4KineticTrack *>::iterator i;
1557	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1558	{
1559	G4KineticTrack * kt = *i;
1560	if( kt->GetState() == G4KineticTrack::inside )
1561	{
1562	nsec++;
1563	G4double tStep(0), tdummy(0);
1564	((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep);
1565	#ifdef debug_BIC_StepParticlesOut
1566	G4cout << " minTimeStep, tStep Particle " <<minTimeStep << " " <<tStep
1567	<< " " <<kt->GetDefinition()->GetParticleName()
1568	<< " 4mom " << kt->GetTrackingMomentum()<<G4endl;
1569	#endif
1570	if(tStep<minTimeStep && tStep> 0 )
1571	{
1572	minTimeStep = tStep;
1573	}
1574	} else if ( kt->GetState() != G4KineticTrack::outside ){
1575	PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
1576	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1577	}
1578	}
1579	minTimeStep *= 1.2;
1580	// G4cerr << "CaptureCount = "<<counter<<" "<<nsec<<" "<<minTimeStep<<" "<<1*ns<<G4endl;
1581	G4double timeToCollision=DBL_MAX;
1582	G4CollisionInitialState * nextCollision=0;
1583	if(theCollisionMgr->Entries() > 0)
1584	{
1585	nextCollision = theCollisionMgr->GetNextCollision();
1586	timeToCollision = nextCollision->GetCollisionTime()-theCurrentTime;
1587	G4cout << " NextCollision * , Time= " << nextCollision << " "
1588	<<timeToCollision<< G4endl;
1589	}
1590	if ( timeToCollision > minTimeStep )
1591	{
1592	DoTimeStep(minTimeStep);
1593	++counter;
1594	} else
1595	{
1596	if (!DoTimeStep(timeToCollision) )
1597	{
1598	// Check if nextCollision is still valid, ie. partcile did not leave nucleus
1599	if (theCollisionMgr->GetNextCollision() != nextCollision )
1600	{
1601	nextCollision = 0;
1602	}
1603	}
1604	// G4cerr <<"post- DoTimeStep 3"<<G4endl;
1605
1606	if(nextCollision)
1607	{
1608	if ( ApplyCollision(nextCollision))
1609	{
1610	// G4cout << "ApplyCollision sucess " << G4endl;
1611	} else
1612	{
1613	theCollisionMgr->RemoveCollision(nextCollision);
1614	}
1615	}
1616	}
1617
1618	if(countreset>100)
1619	{
1620	#ifdef debug_G4BinaryCascade
1621	G4cerr << "G4BinaryCascade.cc: Warning - aborting looping particle(s)" << G4endl;
1622	#endif
1623
1624	// add left secondaries to FinalSate
1625	std::vector<G4KineticTrack *>::iterator iter;
1626	for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
1627	{
1628	theFinalState.push_back(*iter);
1629	}
1630	theSecondaryList.clear();
1631
1632	break;
1633	}
1634
1635	if(Absorb())
1636	{
1637	// haveProducts = true;
1638	// G4cout << "Absorb sucess " << G4endl;
1639	}
1640
1641	if(Capture(false))
1642	{
1643	// haveProducts = true;
1644	#ifdef debug_BIC_StepParticlesOut
1645	G4cout << "Capture sucess " << G4endl;
1646	#endif
1647	}
1648	if ( counter > 100 && theCollisionMgr->Entries() == 0) // no collision, and stepping a while....
1649	{
1650	#ifdef debug_1_BinaryCascade
1651	PrintKTVector(&theSecondaryList,std::string("stepping 100 steps"));
1652	#endif
1653	FindCollisions(&theSecondaryList);
1654	counter=0;
1655	++countreset;
1656	}
1657	}
1658	// G4cerr <<"Finished capture loop "<<G4endl;
1659
1660	//G4cerr <<"pre- DoTimeStep 4"<<G4endl;
1661	DoTimeStep(DBL_MAX);
1662	//G4cerr <<"post- DoTimeStep 4"<<G4endl;
1663
1664
1665	}
1666
1667	//----------------------------------------------------------------------------
1668	void G4BinaryCascade::CorrectFinalPandE()
1669	//----------------------------------------------------------------------------
1670	//
1671	// Modify momenta of outgoing particles.
1672	// Assume two body decay, nucleus(@nominal mass) + sum of final state particles(SFSP).
1673	// momentum of SFSP shall be less than momentum for two body decay.
1674	//
1675	{
1676	#ifdef debug_BIC_CorrectFinalPandE
1677	G4cerr << "BIC: -CorrectFinalPandE called" << G4endl;
1678	#endif
1679
1680	if ( theFinalState.size() == 0 ) return;
1681
1682	G4KineticTrackVector::iterator i;
1683	G4LorentzVector pNucleus=GetFinal4Momentum();
1684	if ( pNucleus.e() == 0 ) return; // check against explicit 0 from GetNucleus4Momentum()
1685	#ifdef debug_BIC_CorrectFinalPandE
1686	G4cerr << " -CorrectFinalPandE 3" << G4endl;
1687	#endif
1688	G4LorentzVector pFinals(0);
1689	G4int nFinals(0);
1690	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1691	{
1692	pFinals += (*i)->Get4Momentum();
1693	++nFinals;
1694	#ifdef debug_BIC_CorrectFinalPandE
1695	G4cout <<"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
1696	<< " 4mom " << (*i)->Get4Momentum()<< G4endl;
1697	#endif
1698	}
1699	#ifdef debug_BIC_CorrectFinalPandE
1700	G4cout << "CorrectFinalPandE pN pF: " <<pNucleus << " " <<pFinals << G4endl;
1701	#endif
1702	G4LorentzVector pCM=pNucleus + pFinals;
1703
1704	G4LorentzRotation toCMS(-pCM.boostVector());
1705	pFinals *=toCMS;
1706
1707	#ifdef debug_BIC_CorrectFinalPandE
1708	G4cout << "CorrectFinalPandE pCM, CMS pCM " << pCM << " " <<toCMS*pCM<< G4endl;
1709	G4cout << "CorrectFinal CMS pN pF " <<toCMS*pNucleus << " "
1710	<<pFinals << G4endl
1711	<< " nucleus initial mass : " <<GetFinal4Momentum().mag()
1712	<<" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus << " " <<pNucleus.mag()<< " "
1713	<< pFinals.mag() << " " << pCM.mag() << G4endl;
1714	#endif
1715
1716	G4LorentzRotation toLab = toCMS.inverse();
1717
1718	G4double s = pCM.mag2();
1719	// G4double m10 = massInNucleus; //pNucleus.mag();
1720	G4double m10 = GetIonMass(currentZ,currentA);
1721	G4double m20 = pFinals.mag();
1722	if( s-(m10+m20)*(m10+m20) < 0 )
1723	{
1724	#ifdef debug_BIC_CorrectFinalPandE
1725	G4cout << "G4BinaryCascade::CorrectFinalPandE() : error! " << G4endl;
1726
1727	G4cout << "not enough mass to correct: mass, A,Z, mass(nucl), mass(finals) "
1728	<< std::sqrt(-s+(m10+m20)*(m10+m20)) << " "
1729	<< currentA << " " << currentZ << " "
1730	<< m10 << " " << m20
1731	<< G4endl;
1732	G4cerr << " -CorrectFinalPandE 4" << G4endl;
1733
1734	PrintKTVector(&theFinalState," mass problem");
1735	#endif
1736	return;
1737	}
1738
1739	// Three momentum in cm system
1740	G4double pInCM = std::sqrt((s-(m10+m20)(m10+m20))(s-(m10-m20)(m10-m20))/(4.s));
1741	#ifdef debug_BIC_CorrectFinalPandE
1742	G4cout <<" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM
1743	<< " " << (pFinals).vect().mag()<< " " << pInCM/(pFinals).vect().mag() << G4endl;
1744	#endif
1745	if ( pFinals.vect().mag() > pInCM )
1746	{
1747	G4ThreeVector p3finals=pInCM*pFinals.vect().unit();
1748
1749	// G4ThreeVector deltap=(p3finals - pFinals.vect() ) / nFinals;
1750	G4double factor=std::max(0.98,pInCM/pFinals.vect().mag()); // small correction
1751	G4LorentzVector qFinals(0);
1752	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1753	{
1754	// G4ThreeVector p3((toCMS(i)->Get4Momentum()).vect() + deltap);
1755	G4ThreeVector p3(factor(toCMS(*i)->Get4Momentum()).vect());
1756	G4LorentzVector p(p3,std::sqrt((*i)->Get4Momentum().mag2() + p3.mag2()));
1757	qFinals += p;
1758	p *= toLab;
1759	#ifdef debug_BIC_CorrectFinalPandE
1760	G4cout << " final p corrected: " << p << G4endl;
1761	#endif
1762	(*i)->Set4Momentum(p);
1763	}
1764	#ifdef debug_BIC_CorrectFinalPandE
1765	G4cout << "CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() << " "
1766	<<GetFinal4Momentum().mag() << G4endl
1767	<< " CMS pFinals , mag, 3.mag : " << qFinals << " " << qFinals.mag() << " " << qFinals.vect().mag()<< G4endl;
1768	G4cerr << " -CorrectFinalPandE 5 " << factor << G4endl;
1769	#endif
1770	}
1771	#ifdef debug_BIC_CorrectFinalPandE
1772	else { G4cerr << " -CorrectFinalPandE 6 - no correction done" << G4endl; }
1773	#endif
1774
1775	}
1776
1777	//----------------------------------------------------------------------------
1778	void G4BinaryCascade::UpdateTracksAndCollisions(
1779	//----------------------------------------------------------------------------
1780	G4KineticTrackVector * oldSecondaries,
1781	G4KineticTrackVector * oldTarget,
1782	G4KineticTrackVector * newSecondaries)
1783	{
1784	// G4cout << "Entering ... "<<oldTarget<<G4endl;
1785	std::vector<G4KineticTrack *>::iterator iter1, iter2;
1786
1787	// remove old secondaries from the secondary list
1788	if(oldSecondaries)
1789	{
1790	if(!oldSecondaries->empty())
1791	{
1792	for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
1793	++iter1)
1794	{
1795	iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(),
1796	*iter1);
1797	if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2);
1798	}
1799	theCollisionMgr->RemoveTracksCollisions(oldSecondaries);
1800	}
1801	}
1802
1803	// remove old target from the target list
1804	if(oldTarget)
1805	{
1806	// G4cout << "################## Debugging 0 "<<G4endl;
1807	if(oldTarget->size()!=0)
1808	{
1809
1810	// G4cout << "################## Debugging 1 "<<oldTarget->size()<<G4endl;
1811	for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
1812	{
1813	iter2 = std::find(theTargetList.begin(), theTargetList.end(),
1814	*iter1);
1815	theTargetList.erase(iter2);
1816	}
1817	theCollisionMgr->RemoveTracksCollisions(oldTarget);
1818	}
1819	}
1820
1821	if(newSecondaries)
1822	{
1823	if(!newSecondaries->empty())
1824	{
1825	// insert new secondaries in the secondary list
1826	for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
1827	++iter1)
1828	{
1829	theSecondaryList.push_back(*iter1);
1830	}
1831	// look for collisions of new secondaries
1832	FindCollisions(newSecondaries);
1833	}
1834	}
1835	// G4cout << "Exiting ... "<<oldTarget<<G4endl;
1836	}
1837
1838
1839	class SelectFromKTV
1840	{
1841	private:
1842	G4KineticTrackVector * ktv;
1843	G4KineticTrack::CascadeState wanted_state;
1844	public:
1845	SelectFromKTV(G4KineticTrackVector * out, G4KineticTrack::CascadeState astate)
1846	:
1847	ktv(out), wanted_state(astate)
1848	{};
1849	void operator() (G4KineticTrack *& kt) const
1850	{
1851	if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt);
1852	};
1853	};
1854
1855
1856
1857	//----------------------------------------------------------------------------
1858	G4bool G4BinaryCascade::DoTimeStep(G4double theTimeStep)
1859	//----------------------------------------------------------------------------
1860	{
1861
1862	#ifdef debug_BIC_DoTimeStep
1863	G4ping debug("debug_G4BinaryCascade");
1864	debug.push_back("======> DoTimeStep 1"); debug.dump();
1865	G4cerr <<"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
1866	<< " , time="<<theCurrentTime << G4endl;
1867	PrintKTVector(&theSecondaryList, std::string("DoTimeStep - theSecondaryList"));
1868	//PrintKTVector(&theTargetList, std::string("DoTimeStep - theTargetList"));
1869	#endif
1870
1871	G4bool success=true;
1872	std::vector<G4KineticTrack *>::iterator iter;
1873
1874	G4KineticTrackVector * kt_outside = new G4KineticTrackVector;
1875	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1876	SelectFromKTV(kt_outside,G4KineticTrack::outside));
1877	//PrintKTVector(kt_outside, std::string("DoTimeStep - found outside"));
1878
1879	G4KineticTrackVector * kt_inside = new G4KineticTrackVector;
1880	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1881	SelectFromKTV(kt_inside, G4KineticTrack::inside));
1882	// PrintKTVector(kt_inside, std::string("DoTimeStep - found inside"));
1883	//-----
1884	G4KineticTrackVector dummy; // needed for re-usability
1885	#ifdef debug_BIC_DoTimeStep
1886	G4cout << "NOW WE ARE ENTERING THE TRANSPORT"<<G4endl;
1887	#endif
1888
1889	// =================== Here we move the particles ===================
1890
1891	thePropagator->Transport(theSecondaryList, dummy, theTimeStep);
1892
1893	// =================== Here we move the particles ===================
1894
1895	//------
1896
1897	theMomentumTransfer += thePropagator->GetMomentumTransfer();
1898	#ifdef debug_BIC_DoTimeStep
1899	G4cout << "DoTimeStep : theMomentumTransfer = " << theMomentumTransfer << G4endl;
1900	PrintKTVector(&theSecondaryList, std::string("DoTimeStep - secondaries aft trsprt"));
1901	#endif
1902
1903	// Partclies which went INTO nucleus
1904
1905	G4KineticTrackVector * kt_gone_in = new G4KineticTrackVector;
1906	std::for_each( kt_outside->begin(),kt_outside->end(),
1907	SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
1908	// PrintKTVector(kt_gone_in, std::string("DoTimeStep - gone in"));
1909
1910
1911	// Partclies which went OUT OF nucleus
1912	G4KineticTrackVector * kt_gone_out = new G4KineticTrackVector;
1913	std::for_each( kt_inside->begin(),kt_inside->end(),
1914	SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
1915
1916	// PrintKTVector(kt_gone_out, std::string("DoTimeStep - gone out"));
1917
1918	G4KineticTrackVector *fail=CorrectBarionsOnBoundary(kt_gone_in,kt_gone_out);
1919
1920	if ( fail )
1921	{
1922	// some particle(s) supposed to enter/leave were miss_nucleus/captured by the correction
1923	kt_gone_in->clear();
1924	std::for_each( kt_outside->begin(),kt_outside->end(),
1925	SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
1926
1927	kt_gone_out->clear();
1928	std::for_each( kt_inside->begin(),kt_inside->end(),
1929	SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
1930
1931	#ifdef debug_BIC_DoTimeStep
1932	PrintKTVector(fail,std::string(" Failed to go in/out -> miss_nucleus/captured"));
1933	PrintKTVector(kt_gone_in, std::string("recreated kt_gone_in"));
1934	PrintKTVector(kt_gone_out, std::string("recreated kt_gone_out"));
1935	#endif
1936	delete fail;
1937	}
1938
1939	// Add tracks missing nucleus and tracks going straight though to addFinals
1940	std::for_each( kt_outside->begin(),kt_outside->end(),
1941	SelectFromKTV(kt_gone_out,G4KineticTrack::miss_nucleus));
1942	//PrintKTVector(kt_gone_out, std::string("miss to append to final state.."));
1943	std::for_each( kt_outside->begin(),kt_outside->end(),
1944	SelectFromKTV(kt_gone_out,G4KineticTrack::gone_out));
1945
1946	#ifdef debug_BIC_DoTimeStep
1947	PrintKTVector(kt_gone_out, std::string("append to final state.."));
1948	#endif
1949
1950	theFinalState.insert(theFinalState.end(),
1951	kt_gone_out->begin(),kt_gone_out->end());
1952
1953	// Partclies which could not leave nucleus, captured...
1954	G4KineticTrackVector * kt_captured = new G4KineticTrackVector;
1955	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1956	SelectFromKTV(kt_captured, G4KineticTrack::captured));
1957
1958	// Check no track is part in next collision, ie.
1959	// this step was to far, and collisions should not occur any more
1960
1961	if ( theCollisionMgr->Entries()> 0 )
1962	{
1963	if (kt_gone_out->size() )
1964	{
1965	G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
1966	iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
1967	if ( iter != kt_gone_out->end() )
1968	{
1969	success=false;
1970	#ifdef debug_BIC_DoTimeStep
1971	G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
1972	#endif
1973	}
1974	}
1975	if ( kt_captured->size() )
1976	{
1977	G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
1978	iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
1979	if ( iter != kt_captured->end() )
1980	{
1981	success=false;
1982	#ifdef debug_BIC_DoTimeStep
1983	G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
1984	#endif
1985	}
1986	}
1987
1988	}
1989	// PrintKTVector(kt_gone_out," kt_gone_out be4 updatetrack...");
1990	UpdateTracksAndCollisions(kt_gone_out,0 ,0);
1991
1992
1993	if ( kt_captured->size() )
1994	{
1995	theCapturedList.insert(theCapturedList.end(),
1996	kt_captured->begin(),kt_captured->end());
1997	//should be std::for_each(kt_captured->begin(),kt_captured->end(),
1998	// std::mem_fun(&G4KineticTrack::Hit));
1999	// but VC 6 requires:
2000	std::vector<G4KineticTrack *>::iterator i_captured;
2001	for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
2002	{
2003	(*i_captured)->Hit();
2004	}
2005	// PrintKTVector(kt_captured," kt_captured be4 updatetrack...");
2006	UpdateTracksAndCollisions(kt_captured, NULL, NULL);
2007	}
2008
2009	#ifdef debug_G4BinaryCascade
2010	delete kt_inside;
2011	kt_inside = new G4KineticTrackVector;
2012	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2013	SelectFromKTV(kt_inside, G4KineticTrack::inside));
2014	if ( currentZ != (GetTotalCharge(theTargetList)
2015	+ GetTotalCharge(theCapturedList)
2016	+ GetTotalCharge(*kt_inside)) )
2017	{
2018	G4cout << " error-DoTimeStep aft, A, Z: " << currentA << " " << currentZ
2019	<< " sum(tgt,capt,active) "
2020	<< GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside)
2021	<< " targets: " << GetTotalCharge(theTargetList)
2022	<< " captured: " << GetTotalCharge(theCapturedList)
2023	<< " active: " << GetTotalCharge(*kt_inside)
2024	<< G4endl;
2025	}
2026	#endif
2027
2028	delete kt_inside;
2029	delete kt_outside;
2030	delete kt_captured;
2031	delete kt_gone_in;
2032	delete kt_gone_out;
2033
2034	// G4cerr <<"G4BinaryCascade::DoTimeStep: exit "<<G4endl;
2035	theCurrentTime += theTimeStep;
2036
2037	//debug.push_back("======> DoTimeStep 2"); debug.dump();
2038	return success;
2039
2040	}
2041
2042	//----------------------------------------------------------------------------
2043	G4KineticTrackVector* G4BinaryCascade::CorrectBarionsOnBoundary(
2044	G4KineticTrackVector *in,
2045	G4KineticTrackVector *out)
2046	//----------------------------------------------------------------------------
2047	{
2048	G4KineticTrackVector * kt_fail(0);
2049	std::vector<G4KineticTrack *>::iterator iter;
2050	// G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2051	// << currentZ << " "<< currentA << G4endl;
2052	if (in->size())
2053	{
2054	G4int secondaries_in(0);
2055	G4int secondaryBarions_in(0);
2056	G4int secondaryCharge_in(0);
2057	G4double secondaryMass_in(0);
2058
2059	for ( iter =in->begin(); iter != in->end(); ++iter)
2060	{
2061	++secondaries_in;
2062	secondaryCharge_in += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2063	if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
2064	{
2065	secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
2066	if((*iter)->GetDefinition() == G4Neutron::Neutron() \|\|
2067	(*iter)->GetDefinition() == G4Proton::Proton() )
2068	{
2069	secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
2070	} else {
2071	secondaryMass_in += G4Proton::Proton()->GetPDGMass();
2072	}
2073	}
2074	}
2075	G4double mass_initial= GetIonMass(currentZ,currentA);
2076
2077	currentZ += secondaryCharge_in;
2078	currentA += secondaryBarions_in;
2079
2080	// G4cout << "CorrectBarionsOnBoundary,secondaryCharge_in, secondaryBarions_in "
2081	// << secondaryCharge_in << " "<< secondaryBarions_in << G4endl;
2082
2083	G4double mass_final= GetIonMass(currentZ,currentA);
2084
2085	G4double correction= secondaryMass_in + mass_initial - mass_final;
2086	if (secondaries_in>1)
2087	{correction /= secondaries_in;}
2088
2089	#ifdef debug_BIC_CorrectBarionsOnBoundary
2090	G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2091	<< "secondaryCharge_in,secondaryBarions_in,"
2092	<< "energy correction,m_secondry,m_nucl_init,m_nucl_final "
2093	<< currentZ << " "<< currentA <<" "
2094	<< secondaryCharge_in<<" "<<secondaryBarions_in<<" "
2095	<< correction << " "
2096	<< secondaryMass_in << " "
2097	<< mass_initial << " "
2098	<< mass_final << " "
2099	<< G4endl;
2100	PrintKTVector(in,std::string("in be4 correction"));
2101	#endif
2102
2103	for ( iter = in->begin(); iter != in->end(); ++iter)
2104	{
2105	if ((iter)->GetTrackingMomentum().e()+correction > (iter)->GetActualMass())
2106	{
2107	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + correction);
2108	} else {
2109	//particle cannot go in, put to miss_nucleus
2110	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2111	(*iter)->SetState(G4KineticTrack::miss_nucleus);
2112	// Undo correction for Colomb Barrier
2113	G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2114	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + barrier);
2115	if ( ! kt_fail ) kt_fail=new G4KineticTrackVector;
2116	kt_fail->push_back(*iter);
2117	currentZ -= G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2118	currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
2119
2120	}
2121
2122	}
2123	#ifdef debug_BIC_CorrectBarionsOnBoundary
2124	G4cout << " CorrectBarionsOnBoundary, aft, A, Z, sec-Z,A,m,m_in_nucleus "
2125	<< currentA << " " << currentZ << " "
2126	<< secondaryCharge_in << " " << secondaryBarions_in << " "
2127	<< secondaryMass_in << " "
2128	<< G4endl;
2129	PrintKTVector(in,std::string("in AFT correction"));
2130	#endif
2131
2132	}
2133	//----------------------------------------------
2134	if (out->size())
2135	{
2136	G4int secondaries_out(0);
2137	G4int secondaryBarions_out(0);
2138	G4int secondaryCharge_out(0);
2139	G4double secondaryMass_out(0);
2140
2141	for ( iter =out->begin(); iter != out->end(); ++iter)
2142	{
2143	++secondaries_out;
2144	secondaryCharge_out += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2145	if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
2146	{
2147	secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
2148	if((*iter)->GetDefinition() == G4Neutron::Neutron() \|\|
2149	(*iter)->GetDefinition() == G4Proton::Proton() )
2150	{
2151	secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
2152	} else {
2153	secondaryMass_out += G4Neutron::Neutron()->GetPDGMass();
2154	}
2155	}
2156	}
2157
2158	G4double mass_initial= GetIonMass(currentZ,currentA);
2159	currentA -=secondaryBarions_out;
2160	currentZ -=secondaryCharge_out;
2161
2162	// G4cout << "CorrectBarionsOnBoundary,secondaryCharge_out, secondaryBarions_out"
2163	// << secondaryCharge_out << " "<< secondaryBarions_out << G4endl;
2164
2165	// a delta minus will do currentZ < 0 in light nuclei
2166	// if (currentA < 0 \|\| currentZ < 0 )
2167	if (currentA < 0 )
2168	{
2169	G4cerr << "G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
2170	secondaryBarions_out << " " << secondaryCharge_out << G4endl;
2171	PrintKTVector(&theTargetList,"CorrectBarionsOnBoundary Target");
2172	PrintKTVector(&theCapturedList,"CorrectBarionsOnBoundary Captured");
2173	PrintKTVector(&theSecondaryList,"CorrectBarionsOnBoundary Secondaries");
2174	G4cerr << "G4BinaryCascade - currentA, currentZ " << currentA << " " << currentZ << G4endl;
2175	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
2176	}
2177	G4double mass_final=GetIonMass(currentZ,currentA);
2178	G4double correction= mass_initial - mass_final - secondaryMass_out;
2179
2180	if (secondaries_out>1) correction /= secondaries_out;
2181	#ifdef debug_BIC_CorrectBarionsOnBoundary
2182	G4cout << "DoTimeStep,currentZ,currentA,"
2183	<< "secondaries_out,"
2184	<<"secondaryCharge_out,secondaryBarions_out,"
2185	<<"energy correction,m_secondry,m_nucl_init,m_nucl_final "
2186	<< " "<< currentZ << " "<< currentA <<" "
2187	<< secondaries_out << " "
2188	<< secondaryCharge_out<<" "<<secondaryBarions_out<<" "
2189	<< correction << " "
2190	<< secondaryMass_out << " "
2191	<< mass_initial << " "
2192	<< mass_final << " "
2193	<< G4endl;
2194	PrintKTVector(out,std::string("out be4 correction"));
2195	#endif
2196
2197	for ( iter = out->begin(); iter != out->end(); ++iter)
2198	{
2199	if ((iter)->GetTrackingMomentum().e()+correction > (iter)->GetActualMass())
2200	{
2201	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + correction);
2202	} else
2203	{
2204	// particle cannot go out due to change of nuclear potential!
2205	// capture protons and neutrons;
2206	if(((*iter)->GetDefinition() == G4Proton::Proton()) \|\|
2207	((*iter)->GetDefinition() == G4Neutron::Neutron()))
2208	{
2209	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2210	(*iter)->SetState(G4KineticTrack::captured);
2211	// Undo correction for Colomb Barrier
2212	G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2213	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() - barrier);
2214	if ( kt_fail == 0 ) kt_fail=new G4KineticTrackVector;
2215	kt_fail->push_back(*iter);
2216	currentZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2217	currentA += (*iter)->GetDefinition()->GetBaryonNumber();
2218	}
2219	#ifdef debug_BIC_CorrectBarionsOnBoundary
2220	else
2221	{
2222	G4cout << "Not correcting outgoing " << *iter << " "
2223	<< (*iter)->GetDefinition()->GetPDGEncoding() << " "
2224	<< (*iter)->GetDefinition()->GetParticleName() << G4endl;
2225	PrintKTVector(out,std::string("outgoing, one not corrected"));
2226	}
2227	#endif
2228	}
2229	}
2230
2231	#ifdef debug_BIC_CorrectBarionsOnBoundary
2232	PrintKTVector(out,std::string("out AFTER correction"));
2233	G4cout << " DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
2234	<< currentA << " "<< currentZ << " "
2235	<< secondaryCharge_out << " "<< secondaryBarions_out << " "<<
2236	secondaryMass_out << " "
2237	<< massInNucleus << " "
2238	<< G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2239	<< " " << massInNucleus -G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2240	<< G4endl;
2241	#endif
2242	}
2243
2244	return kt_fail;
2245	}
2246
2247
2248	//----------------------------------------------------------------------------
2249
2250	G4Fragment * G4BinaryCascade::FindFragments()
2251	//----------------------------------------------------------------------------
2252	{
2253
2254	G4int a = theTargetList.size()+theCapturedList.size();
2255	#ifdef debug_BIC_FindFragments
2256	G4cout << "target, captured, secondary: "
2257	<< theTargetList.size() << " "
2258	<< theCapturedList.size()<< " "
2259	<< theSecondaryList.size()
2260	<< G4endl;
2261	#endif
2262	G4int zTarget = 0;
2263	G4KineticTrackVector::iterator i;
2264	for(i = theTargetList.begin(); i != theTargetList.end(); ++i)
2265	{
2266	if((*i)->GetDefinition()->GetPDGCharge() == eplus)
2267	{
2268	zTarget++;
2269	}
2270	}
2271
2272	G4int zCaptured = 0;
2273	G4LorentzVector CapturedMomentum=0;
2274	for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2275	{
2276	CapturedMomentum += (*i)->Get4Momentum();
2277	if((*i)->GetDefinition()->GetPDGCharge() == eplus)
2278	{
2279	zCaptured++;
2280	}
2281	}
2282
2283	G4int z = zTarget+zCaptured;
2284
2285	#ifdef debug_G4BinaryCascade
2286	if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) )
2287	{
2288	G4cout << " FindFragment Counting error z a " << z << " " <<a << " "
2289	<< GetTotalCharge(theTargetList) << " " << GetTotalCharge(theCapturedList)<<
2290	G4endl;
2291	PrintKTVector(&theTargetList, std::string("theTargetList"));
2292	PrintKTVector(&theCapturedList, std::string("theCapturedList"));
2293	}
2294	#endif
2295	//debug
2296	/*
2297	* G4cout << " Fragment mass table / real "
2298	* << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
2299	* << " / " << GetFinal4Momentum().mag()
2300	* << " difference "
2301	* << GetFinal4Momentum().mag() - G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
2302	* << G4endl;
2303	*/
2304	//
2305	// if(getenv("BCDEBUG") ) G4cerr << "Fragment A, Z "<< a <<" "<< z<<G4endl;
2306	if ( z < 1 ) return 0;
2307
2308	G4int holes = the3DNucleus->GetMassNumber() - theTargetList.size();
2309	G4int excitons = theCapturedList.size();
2310	#ifdef debug_BIC_FindFragments
2311	G4cout << "Fragment: a= " << a
2312	<< " z= " << z
2313	<< " particles= " << excitons
2314	<< " Charged= " << zCaptured
2315	<< " holes= " << holes
2316	<< " excitE= " <<GetExcitationEnergy()
2317	<< " Final4Momentum= " << GetFinalNucleusMomentum()
2318	<< " capturMomentum= " << CapturedMomentum
2319	<< G4endl;
2320	#endif
2321
2322	G4Fragment * fragment = new G4Fragment(a,z,GetFinalNucleusMomentum());
2323	fragment->SetNumberOfHoles(holes);
2324
2325	//GF fragment->SetNumberOfParticles(excitons-holes);
2326	fragment->SetNumberOfParticles(excitons);
2327	fragment->SetNumberOfCharged(zCaptured);
2328	G4ParticleDefinition * aIonDefinition =
2329	G4ParticleTable::GetParticleTable()->FindIon(a,z,0,z);
2330	fragment->SetParticleDefinition(aIonDefinition);
2331
2332	return fragment;
2333	}
2334
2335	//----------------------------------------------------------------------------
2336
2337	G4LorentzVector G4BinaryCascade::GetFinal4Momentum()
2338	//----------------------------------------------------------------------------
2339	{
2340	// the initial 3-momentum will differ from 0, if nucleus created by string model.
2341	G4LorentzVector final4Momentum = theInitial4Mom;
2342	G4KineticTrackVector::iterator i;
2343	for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
2344	{
2345	final4Momentum += (*i)->GetTrackingMomentum();
2346	//G4cout << "Initial state: "<<(*i)->Get4Momentum()<<G4endl;
2347	}
2348	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2349	{
2350	final4Momentum -= (*i)->Get4Momentum();
2351	}
2352
2353	if((final4Momentum.vect()/final4Momentum.e()).mag()>1.0 && currentA > 0)
2354	{
2355	# ifdef debug_BIC_Final4Momentum
2356	G4cerr << G4endl;
2357	G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl;
2358	G4KineticTrackVector::iterator i;
2359	G4cerr <<" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<G4endl;
2360	for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
2361	{
2362	G4cerr << " Initial state (get4M), (trackingM): "
2363	<<(i)->Get4Momentum()<<", " << (i)->GetTrackingMomentum() <<","
2364	<<(*i)->GetDefinition()->GetParticleName()<<G4endl;
2365	}
2366	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2367	{
2368	G4cerr <<" Final state: "<<(i)->Get4Momentum()<<(i)->GetDefinition()->GetParticleName()<<G4endl;
2369	}
2370	G4cerr<< " Final4Momentum = "<<final4Momentum <<" "<<final4Momentum.m()<<G4endl;
2371	G4cerr <<" current A, Z = "<< currentA<<", "<<currentZ<<G4endl;
2372	G4cerr << G4endl;
2373	# endif
2374
2375	final4Momentum=G4LorentzVector(0);
2376	}
2377	return final4Momentum;
2378	}
2379
2380	//----------------------------------------------------------------------------
2381	G4LorentzVector G4BinaryCascade::GetFinalNucleusMomentum()
2382	//----------------------------------------------------------------------------
2383	{
2384	// return momentum of nucleus for use with precompound model; also keep transformation to
2385	// apply to precompoud products.
2386
2387	G4LorentzVector CapturedMomentum=0;
2388	G4KineticTrackVector::iterator i;
2389	// G4cout << "GetFinalNucleusMomentum Captured size: " <<theCapturedList.size() << G4endl;
2390	for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2391	{
2392	CapturedMomentum += (*i)->Get4Momentum();
2393	}
2394	//G4cout << "GetFinalNucleusMomentum CapturedMomentum= " <<CapturedMomentum << G4endl;
2395	// G4cerr << "it 9"<<G4endl;
2396
2397	G4LorentzVector NucleusMomentum = GetFinal4Momentum();
2398	if ( NucleusMomentum.e() > 0 )
2399	{
2400	// G4cout << "GetFinalNucleusMomentum GetFinal4Momentum= " <<NucleusMomentum <<" "<<NucleusMomentum.mag()<<G4endl;
2401	// boost nucleus to a frame such that the momentum of nucleus == momentum of Captured
2402	G4ThreeVector boost= (NucleusMomentum.vect() -CapturedMomentum.vect())/NucleusMomentum.e();
2403	if(boost.mag2()>1.0)
2404	{
2405	# ifdef debug_BIC_FinalNucleusMomentum
2406	G4cerr << "G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<G4endl;
2407	G4cerr << "it 0"<<boost <<G4endl;
2408	G4cerr << "it 01"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2409	G4cout <<" testing boost "<<boost<<" "<<boost.mag()<<G4endl;
2410	# endif
2411	boost=G4ThreeVector(0);
2412	NucleusMomentum=G4LorentzVector(0);
2413	}
2414	G4LorentzRotation nucleusBoost( -boost );
2415	precompoundLorentzboost.set( boost );
2416	#ifdef debug_debug_BIC_FinalNucleusMomentum
2417	G4cout << "GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2418	#endif
2419	NucleusMomentum *= nucleusBoost;
2420	#ifdef debug_BIC_FinalNucleusMomentum
2421	G4cout << "GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<G4endl;
2422	#endif
2423	}
2424	return NucleusMomentum;
2425	}
2426
2427	//----------------------------------------------------------------------------
2428	G4ReactionProductVector * G4BinaryCascade::Propagate1H1(
2429	//----------------------------------------------------------------------------
2430	G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus)
2431	{
2432	G4ReactionProductVector * products = new G4ReactionProductVector;
2433	G4ParticleDefinition * aHTarg = G4Proton::ProtonDefinition();
2434	G4double mass = aHTarg->GetPDGMass();
2435	if (nucleus->GetCharge() == 0) aHTarg = G4Neutron::NeutronDefinition();
2436	mass = aHTarg->GetPDGMass();
2437	G4KineticTrackVector * secs = 0;
2438	G4ThreeVector pos(0,0,0);
2439	G4LorentzVector mom(mass);
2440	G4KineticTrack aTarget(aHTarg, 0., pos, mom);
2441	G4bool done(false);
2442	std::vector<G4KineticTrack *>::iterator iter, jter;
2443	// data member static G4Scatterer theH1Scatterer;
2444	//G4cout << " start 1H1 for " << (*secondaries).front()->GetDefinition()->GetParticleName()
2445	// << " on " << aHTarg->GetParticleName() << G4endl;
2446	G4int tryCount(0);
2447	while(!done && tryCount++ <200)
2448	{
2449	if(secs)
2450	{
2451	std::for_each(secs->begin(), secs->end(), DeleteKineticTrack());
2452	delete secs;
2453	}
2454	secs = theH1Scatterer->Scatter((secondaries).front(), aTarget);
2455	for(size_t ss=0; secs && ss<secs->size(); ss++)
2456	{
2457	// G4cout << "1H1 " << (*secs)[ss]->GetDefinition()->GetParticleName()
2458	// << ", shortlived? "<< (*secs)[ss]->GetDefinition()->IsShortLived()<< G4endl;
2459	if((*secs)[ss]->GetDefinition()->IsShortLived()) done = true;
2460	}
2461	// G4cout << G4endl;
2462	}
2463	size_t current(0);
2464	for(current=0; current<secs->size(); current++)
2465	{
2466	if((*secs)[current]->GetDefinition()->IsShortLived())
2467	{
2468	G4KineticTrackVector * dec = (*secs)[current]->Decay();
2469	for(jter=dec->begin(); jter != dec->end(); jter++)
2470	{
2471	//G4cout << "Decay"<<G4endl;
2472	secs->push_back(*jter);
2473	//G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<<G4endl;
2474	}
2475	delete (*secs)[current];
2476	delete dec;
2477	}
2478	else
2479	{
2480	//G4cout << "FS "<<G4endl;
2481	//G4cout << "FS "<<(*secs)[current]->GetDefinition()->GetParticleName()<<G4endl;
2482	theFinalState.push_back((*secs)[current]);
2483	}
2484	}
2485	//G4cout << "Through loop"<<G4endl;
2486	delete secs;
2487	for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2488	{
2489	G4KineticTrack * kt = *iter;
2490	G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
2491	aNew->SetMomentum(kt->Get4Momentum().vect());
2492	aNew->SetTotalEnergy(kt->Get4Momentum().e());
2493	products->push_back(aNew);
2494	#ifdef debug_1_BinaryCascade
2495	if (! kt->GetDefinition()->GetPDGStable() )
2496	{
2497	if (kt->GetDefinition()->IsShortLived())
2498	{
2499	G4cout << "final shortlived : ";
2500	} else
2501	{
2502	G4cout << "final un stable : ";
2503	}
2504	G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
2505	}
2506	#endif
2507	delete kt;
2508	}
2509	theFinalState.clear();
2510	return products;
2511
2512	}
2513
2514	//----------------------------------------------------------------------------
2515	G4ThreeVector G4BinaryCascade::GetSpherePoint(
2516	G4double r, const G4LorentzVector & mom4)
2517	//----------------------------------------------------------------------------
2518	{
2519	// Get a point outside radius.
2520	// point is random in plane (circle of radius r) orthogonal to mom,
2521	// plus -1rmom->vect()->unit();
2522	G4ThreeVector o1, o2;
2523	G4ThreeVector mom = mom4.vect();
2524
2525	o1= mom.orthogonal(); // we simply need any vector non parallel
2526	o2= mom.cross(o1); // o2 is now orthogonal to mom and o1, ie. o1 and o2 define plane.
2527
2528	G4double x2, x1;
2529
2530	do
2531	{
2532	x1=(G4UniformRand()-.5)*2;
2533	x2=(G4UniformRand()-.5)*2;
2534	} while (sqr(x1) +sqr(x2) > 1.);
2535
2536	return G4ThreeVector(r(x1o1.unit() + x2o2.unit() - 1.5 mom.unit()));
2537
2538
2539
2540	/*
2541	* // Get a point uniformly distributed on the surface of a sphere,
2542	* // with z < 0.
2543	* G4double b = r*G4UniformRand(); // impact parameter
2544	* G4double phi = G4UniformRand()2pi;
2545	* G4double x = b*std::cos(phi);
2546	* G4double y = b*std::sin(phi);
2547	* G4double z = -std::sqrt(rr-bb);
2548	* z *= 1.001; // Get position a little bit out of the sphere...
2549	* point.setX(x);
2550	* point.setY(y);
2551	* point.setZ(z);
2552	*/
2553	}
2554
2555	//----------------------------------------------------------------------------
2556
2557	void G4BinaryCascade::ClearAndDestroy(G4KineticTrackVector * ktv)
2558	//----------------------------------------------------------------------------
2559	{
2560	std::vector<G4KineticTrack *>::iterator i;
2561	for(i = ktv->begin(); i != ktv->end(); ++i)
2562	delete (*i);
2563	ktv->clear();
2564	}
2565
2566	//----------------------------------------------------------------------------
2567	void G4BinaryCascade::ClearAndDestroy(G4ReactionProductVector * rpv)
2568	//----------------------------------------------------------------------------
2569	{
2570	std::vector<G4ReactionProduct *>::iterator i;
2571	for(i = rpv->begin(); i != rpv->end(); ++i)
2572	delete (*i);
2573	rpv->clear();
2574	}
2575
2576	//----------------------------------------------------------------------------
2577	void G4BinaryCascade::PrintKTVector(G4KineticTrackVector * ktv, std::string comment)
2578	//----------------------------------------------------------------------------
2579	{
2580	if (comment.size() > 0 ) G4cout << comment << G4endl;
2581	G4cout << " vector: " << ktv << ", number of tracks: " << ktv->size()
2582	<< G4endl;
2583	std::vector<G4KineticTrack *>::iterator i;
2584	G4int count;
2585
2586	for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
2587	{
2588	G4KineticTrack * kt = *i;
2589	G4cout << " track n. " << count;
2590	PrintKTVector(kt);
2591	}
2592	}
2593	//----------------------------------------------------------------------------
2594	void G4BinaryCascade::PrintKTVector(G4KineticTrack * kt, std::string comment)
2595	//----------------------------------------------------------------------------
2596	{
2597	if (comment.size() > 0 ) G4cout << comment << G4endl;
2598	G4cout << ", id: " << kt << G4endl;
2599	G4ThreeVector pos = kt->GetPosition();
2600	G4LorentzVector mom = kt->Get4Momentum();
2601	G4LorentzVector tmom = kt->GetTrackingMomentum();
2602	G4ParticleDefinition * definition = kt->GetDefinition();
2603	G4cout << " definition: " << definition->GetPDGEncoding() << " pos: "
2604	<< 1/fermipos << " R: " << 1/fermipos.mag() << " 4mom: "
2605	<< 1/MeVmom <<"Tr_mom" << 1/MeVtmom << " P: " << 1/MeV*mom.vect().mag()
2606	<< " M: " << 1/MeV*mom.mag() << G4endl;
2607	G4cout <<" trackstatus: "<<kt->GetState()<<G4endl;
2608	}
2609
2610	//----------------------------------------------------------------------------
2611	G4bool G4BinaryCascade::CheckDecay(G4KineticTrackVector * products)
2612	//----------------------------------------------------------------------------
2613	{
2614	G4int A = the3DNucleus->GetMassNumber();
2615	G4int Z = the3DNucleus->GetCharge();
2616
2617	G4FermiMomentum fermiMom;
2618	fermiMom.Init(A, Z);
2619
2620	const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
2621
2622	G4KineticTrackVector::iterator i;
2623	G4ParticleDefinition * definition;
2624
2625	// ------ debug
2626	G4bool myflag = true;
2627	// ------ end debug
2628	for(i = products->begin(); i != products->end(); ++i)
2629	{
2630	definition = (*i)->GetDefinition();
2631	if((definition->GetParticleName() != "delta++" )&&
2632	(definition->GetParticleName() != "delta+" )&&
2633	(definition->GetParticleName() != "delta0" )&&
2634	(definition->GetParticleName() != "delta-" ))
2635	continue;
2636	G4ThreeVector pos = (*i)->GetPosition();
2637	G4double d = density->GetDensity(pos);
2638	G4double pFermi= fermiMom.GetFermiMomentum(d);
2639	G4LorentzVector mom = (*i)->Get4Momentum();
2640	G4LorentzRotation boost(mom.boostVector());
2641	G4ThreeVector pion3(227MeV mom.vect().unit()); // 227 is decay product in rest frame
2642	G4LorentzVector pion(pion3, std::sqrt(sqr(140*MeV) +pion3.mag()));
2643	// G4cout << "pi rest " << pion << G4endl;
2644	pion = boost * pion;
2645	// G4cout << "pi lab " << pion << G4endl;
2646	// ------ debug
2647	// G4cout << "p:[" << (1/MeV)pion.x() << " " << (1/MeV)pion.y() << " "
2648	// << (1/MeV)*pion.z() << "] \|p3\|:"
2649	// << (1/MeV)pion.vect().mag() << " E:" << (1/MeV)pion.t() << " ] m: "
2650	// << (1/MeV)*pion.mag() << " pos["
2651	// << (1/fermi)pos.x() << " "<< (1/fermi)pos.y() << " "
2652	// << (1/fermi)*pos.z() << "] \|Dpos\|: "
2653	// << " Pfermi:"
2654	// << (1/MeV)*pFermi << G4endl;
2655	// ------ end debug
2656
2657	if( pion.vect().mag() < pFermi )
2658	{
2659	// ------ debug
2660	myflag = false;
2661	// ------ end debug
2662	}
2663	}
2664	return myflag;
2665	// return true;
2666	}
2667
2668	//----------------------------------------------------------------------------
2669	G4double G4BinaryCascade::GetIonMass(G4int Z, G4int A)
2670	//----------------------------------------------------------------------------
2671	{
2672	G4double mass(0);
2673	if ( Z > 0 && A >= Z )
2674	{
2675	mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
2676
2677	} else if ( A > 0 && Z>0 )
2678	{
2679	// charge Z > A; will happen for light nuclei with pions involved.
2680	mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(A,A);
2681
2682	} else if ( A >= 0 && Z<=0 )
2683	{
2684	// all neutral, or empty nucleus
2685	mass = A * G4Neutron::Neutron()->GetPDGMass();
2686
2687	} else if ( A == 0 && std::abs(Z)<2 )
2688	{
2689	// empty nucleus, except maybe pions
2690	mass = 0;
2691
2692	} else
2693	{
2694	G4cerr << "G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
2695	<< A << "," << Z << ")" <<G4endl;
2696	G4Exception("G4BinaryCascade::GetIonMass() - giving up");
2697	}
2698	return mass;
2699	}
2700
2701	void G4BinaryCascade::PrintWelcomeMessage()
2702	{
2703	G4cout <<"Thank you for using G4BinaryCascade. "<<G4endl;
2704	}
2705
2706	//----------------------------------------------------------------------------
2707	void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision,
2708	G4KineticTrackVector * products)
2709	{
2710	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2711
2712	G4KineticTrackVector debug1;
2713	debug1.push_back(collision->GetPrimary());
2714	PrintKTVector(&debug1,std::string(" Primary particle"));
2715	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
2716	PrintKTVector(products,std::string(" Scatterer products"));
2717
2718	G4double thisExcitation(0);
2719	// excitation energy from this collision
2720	// initial state:
2721	G4double initial(0);
2722	G4KineticTrack * kt=collision->GetPrimary();
2723	initial += kt->Get4Momentum().e();
2724
2725	initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2726	initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2727	G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
2728	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2729	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2730	<< " " << initial << G4endl;;
2731
2732	G4KineticTrackVector ktv=collision->GetTargetCollection();
2733	for ( unsigned int it=0; it < ktv.size(); it++)
2734	{
2735	kt=ktv[it];
2736	initial += kt->Get4Momentum().e();
2737	thisExcitation += kt->GetDefinition()->GetPDGMass()
2738	- kt->Get4Momentum().e()
2739	- RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2740	// initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2741	// initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2742	G4cout << "Targ. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
2743	<< " " << kt->Get4Momentum().e()
2744	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2745	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2746	<< " " << initial <<" Excit " << thisExcitation << G4endl;;
2747	}
2748
2749	G4double final(0);
2750	G4double mass_out(0);
2751	G4int product_barions(0);
2752	if ( products )
2753	{
2754	for ( unsigned int it=0; it < products->size(); it++)
2755	{
2756	kt=(*products)[it];
2757	final += kt->Get4Momentum().e();
2758	final += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2759	final += RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2760	if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
2761	mass_out += kt->GetDefinition()->GetPDGMass();
2762	G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
2763	<< " " << kt->Get4Momentum().e()
2764	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2765	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2766	<< " " << final << G4endl;;
2767	}
2768	}
2769
2770
2771	G4int finalA = currentA;
2772	G4int finalZ = currentZ;
2773	if ( products )
2774	{
2775	finalA -= product_barions;
2776	finalZ -= GetTotalCharge(*products);
2777	}
2778	G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2779	- (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA)
2780	+ mass_out);
2781	G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ
2782	<< " delta-mass " << delta<<G4endl;
2783	final+=delta;
2784	mass_out = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA);
2785	G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
2786	<< " " << final << " "
2787	<< mass_out<<" "
2788	<< currentInitialEnergy - final - mass_out
2789	<< G4endl;
2790	currentInitialEnergy-=final;
2791
2792	}
2793
2794	//----------------------------------------------------------------------------
2795	void G4BinaryCascade::DebugEpConservation(const G4HadProjectile & aTrack,
2796	G4ReactionProductVector* products)
2797	{
2798	G4ReactionProductVector::iterator iter;
2799	G4double Efinal(0);
2800	G4ThreeVector pFinal(0);
2801	if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
2802	{
2803	G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
2804	}
2805
2806	for(iter = products->begin(); iter != products->end(); ++iter)
2807	{
2808
2809	// G4cout << " Secondary E - Ekin / p " <<
2810	// (*iter)->GetDefinition()->GetParticleName() << " " <<
2811	// (*iter)->GetTotalEnergy() << " - " <<
2812	// (*iter)->GetKineticEnergy()<< " / " <<
2813	// (*iter)->GetMomentum().x() << " " <<
2814	// (*iter)->GetMomentum().y() << " " <<
2815	// (*iter)->GetMomentum().z() << G4endl;
2816	Efinal += (*iter)->GetTotalEnergy();
2817	pFinal += (*iter)->GetMomentum();
2818	}
2819
2820	// G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
2821	G4cout << "BIC E/p delta " <<
2822	(aTrack.Get4Momentum().e()+ the3DNucleus->GetMass() - Efinal)/MeV <<
2823	" MeV / mom " << (aTrack.Get4Momentum() - pFinal ) /MeV << G4endl;
2824
2825	}
2826
2827	//----------------------------------------------------------------------------

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

Download in other formats: