Context Navigation

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

Last change on this file since 1196 was 1196, checked in by garnier, 15 years ago
update CVS release candidate geant4.9.3.01
File size: 101.2 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.GetN(), aNucleus.GetZ());
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 << " 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	G4double initial_Efermi(0);
1084	G4LorentzVector mom4Primary(0);
1085
1086	if (primary->GetState() == G4KineticTrack::inside)
1087	{
1088	initialBaryon = primary->GetDefinition()->GetBaryonNumber();
1089	initialCharge = G4lrint(primary->GetDefinition()->GetPDGCharge());
1090	}
1091
1092	// for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
1093	G4int PDGcode=std::abs(primary->GetDefinition()->GetPDGEncoding());
1094	mom4Primary=primary->Get4Momentum();
1095	initial_Efermi=RKprop->GetField(primary->GetDefinition()->GetPDGEncoding(),primary->GetPosition());
1096
1097	if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1098	{
1099	initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
1100	primary->GetPosition());
1101	primary->Update4Momentum(mom4Primary.e() - initial_Efermi);
1102	}
1103
1104	std::vector<G4KineticTrack *>::iterator titer;
1105	for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
1106	{
1107	G4ParticleDefinition * aDef=(*titer)->GetDefinition();
1108	G4int aCode=aDef->GetPDGEncoding();
1109	G4ThreeVector aPos=(*titer)->GetPosition();
1110	initial_Efermi+= RKprop->GetField(aCode, aPos);
1111	// initial_Efermi+= RKprop->GetField((titer)->GetDefinition()->GetPDGEncoding(),(titer)->GetPosition());
1112	}
1113	//****************************************
1114
1115	G4KineticTrackVector * products=0;
1116	products = collision->GetFinalState();
1117
1118	#ifdef debug_BIC_ApplyCollision
1119	G4bool havePion=false;
1120	for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1121	{
1122	G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
1123	if (std::abs(PDGcode)==211 \|\| PDGcode==111 ) havePion=true;
1124	}
1125	if ( !products \|\| havePion)
1126	{
1127	G4cout << " Collision " << collision << ", type: "<< typeid(*collision->GetGenerator()).name()
1128	<< ", with NO products! " <<G4endl;
1129	G4cout << G4endl<<"Initial condition are these:"<<G4endl;
1130	G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
1131	PrintKTVector(collision->GetPrimary());
1132	for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
1133	{
1134	G4cout << "targ: "
1135	<<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
1136	}
1137	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
1138	}
1139	#endif
1140	G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
1141	// if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl;
1142	// if ( products ) PrintKTVector(products, " reaction products");
1143	//****************************************
1144
1145	// reset primary to initial state
1146	primary->Set4Momentum(mom4Primary);
1147
1148
1149	G4int lateBaryon(0), lateCharge(0);
1150
1151	if ( lateParticleCollision )
1152	{ // for late particles, reset charges
1153	//G4cout << "lateP, initial B C state " << initialBaryon << " "
1154	// << initialCharge<< " " << primary->GetState() << " "<< primary->GetDefinition()->GetParticleName()<< G4endl;
1155	lateBaryon = initialBaryon;
1156	lateCharge = initialCharge;
1157	initialBaryon=initialCharge=0;
1158	}
1159
1160	initialBaryon += collision->GetTargetBaryonNumber();
1161	initialCharge+=G4lrint(collision->GetTargetCharge());
1162	if(!products \|\| products->size()==0 \|\| !CheckPauliPrinciple(products))
1163	{
1164	#ifdef debug_BIC_ApplyCollision
1165	if (products) G4cout << " ======Failed Pauli =====" << G4endl;
1166	G4cerr << "G4BinaryCascade::ApplyCollision blocked"<<G4endl;
1167	#endif
1168	if (products) ClearAndDestroy(products);
1169	if ( ! haveTarget ) FindDecayCollision(primary); // for decay, sample new decay
1170	delete products;
1171	return false;
1172	}
1173
1174	G4double final_Efermi(0);
1175	G4KineticTrackVector resonances;
1176	for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1177	{
1178	G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
1179	// G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
1180	final_Efermi+=RKprop->GetField(PDGcode,(*i)->GetPosition());
1181	if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1182	{
1183	resonances.push_back(*i);
1184	}
1185	}
1186	if ( resonances.size() > 0 )
1187	{
1188	G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
1189	for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
1190	{
1191	G4LorentzVector mom=(*res)->Get4Momentum();
1192	G4double mass2=mom.mag2();
1193	G4double newEnergy=mom.e() + delta_Fermi;
1194	G4double newEnergy2= newEnergy*newEnergy;
1195	//G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
1196	if ( newEnergy2 < mass2 )
1197	{
1198	ClearAndDestroy(products);
1199	if (target_collection.size() == 0 ) FindDecayCollision(primary); // for decay, sample new decay
1200	delete products;
1201	return false;
1202	}
1203	// G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl;
1204	G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
1205	(*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
1206	}
1207	}
1208
1209	#ifdef debug_BIC_ApplyCollision
1210	DebugApplyCollision(collision, products);
1211	#endif
1212
1213	G4int finalBaryon(0);
1214	G4int finalCharge(0);
1215	G4KineticTrackVector toFinalState;
1216	for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1217	{
1218	if ( ! lateParticleCollision )
1219	{
1220	(*i)->SetState(primary->GetState()); // secondaries are created inside nucleus, except for decay in propagate
1221	finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1222	finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
1223	} else {
1224	G4double tin=0., tout=0.;
1225	if (((G4RKPropagation)thePropagator)->GetSphereIntersectionTimes((i),tin,tout))
1226	{
1227	if ( tin > 0 )
1228	{
1229	(*i)->SetState(G4KineticTrack::outside);
1230	}
1231	else if ( tout > 0 )
1232	{
1233	(*i)->SetState(G4KineticTrack::inside);
1234	finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1235	finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
1236	}
1237	else
1238	{
1239	(*i)->SetState(G4KineticTrack::gone_out);
1240	toFinalState.push_back((*i));
1241	}
1242	} else
1243	{
1244	(*i)->SetState(G4KineticTrack::miss_nucleus);
1245	//G4cout << " G4BC - miss -late Part- no intersection found " << G4endl;
1246	toFinalState.push_back((*i));
1247	}
1248
1249	// G4cout << " PDGcode, state " << (i)->GetDefinition()->GetPDGEncoding() << " " << (i)->GetState()<<G4endl;
1250
1251	}
1252	}
1253	if(!toFinalState.empty())
1254	{
1255	theFinalState.insert(theFinalState.end(),
1256	toFinalState.begin(),toFinalState.end());
1257	std::vector<G4KineticTrack *>::iterator iter1, iter2;
1258	for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
1259	++iter1)
1260	{
1261	iter2 = std::find(products->begin(), products->end(),
1262	*iter1);
1263	if ( iter2 != products->end() ) products->erase(iter2);
1264	}
1265	theCollisionMgr->RemoveTracksCollisions(&toFinalState);
1266	}
1267
1268	//G4cout << " currentA, Z be4: " << currentA << " " << currentZ << G4endl;
1269	currentA += finalBaryon-initialBaryon;
1270	currentZ += finalCharge-initialCharge;
1271	//G4cout << " currentA, Z aft: " << currentA << " " << currentZ << G4endl;
1272
1273	G4KineticTrackVector oldSecondaries;
1274	if (primary)
1275	{
1276	oldSecondaries.push_back(primary);
1277	primary->Hit();
1278	}
1279
1280	#ifdef debug_G4BinaryCascade
1281	if ( (finalBaryon-initialBaryon-lateBaryon) != 0 \|\| (finalCharge-initialCharge-lateCharge) != 0 )
1282	{
1283	G4cout << "G4BinaryCascade: Error in Balancing: " << G4endl;
1284	G4cout << "initial/final baryon number, initial/final Charge "
1285	<< initialBaryon <<" "<< finalBaryon <<" "
1286	<< initialCharge <<" "<< finalCharge <<" "
1287	<< " in Collision type: "<< typeid(*collision->GetGenerator()).name()
1288	<< ", with number of products: "<< products->size() <<G4endl;
1289	G4cout << G4endl<<"Initial condition are these:"<<G4endl;
1290	G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
1291	for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
1292	{
1293	G4cout << "targ: "
1294	<<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
1295	}
1296	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
1297	G4cout << G4endl<<G4endl;
1298	}
1299	#endif
1300
1301	G4KineticTrackVector oldTarget = collision->GetTargetCollection();
1302	for(size_t ii=0; ii< oldTarget.size(); ii++)
1303	{
1304	oldTarget[ii]->Hit();
1305	}
1306
1307	debug.push_back("=======> we have hit nucleons <=======");
1308
1309	UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products);
1310	std::for_each(oldSecondaries.begin(), oldSecondaries.end(), Delete<G4KineticTrack>());
1311	std::for_each(oldTarget.begin(), oldTarget.end(), Delete<G4KineticTrack>());
1312
1313	delete products;
1314	return true;
1315	}
1316
1317
1318	//----------------------------------------------------------------------------
1319	G4bool G4BinaryCascade::Absorb()
1320	//----------------------------------------------------------------------------
1321	{
1322	// Do it in two step: cannot change theSecondaryList inside the first loop.
1323	G4Absorber absorber(theCutOnPAbsorb);
1324
1325	// Build the vector of G4KineticTracks that must be absorbed
1326	G4KineticTrackVector absorbList;
1327	std::vector<G4KineticTrack *>::iterator iter;
1328	// PrintKTVector(&theSecondaryList, " testing for Absorb" );
1329	for(iter = theSecondaryList.begin();
1330	iter != theSecondaryList.end(); ++iter)
1331	{
1332	G4KineticTrack * kt = *iter;
1333	if(kt->GetState() == G4KineticTrack::inside)// absorption happens only inside the nucleus
1334	{
1335	if(absorber.WillBeAbsorbed(*kt))
1336	{
1337	absorbList.push_back(kt);
1338	}
1339	}
1340	}
1341
1342	if(absorbList.empty())
1343	return false;
1344
1345	G4KineticTrackVector toDelete;
1346	for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
1347	{
1348	G4KineticTrack * kt = *iter;
1349	if(!absorber.FindAbsorbers(*kt, theTargetList))
1350	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1351
1352	if(!absorber.FindProducts(*kt))
1353	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1354
1355	G4KineticTrackVector * products = absorber.GetProducts();
1356	// ------ debug
1357	G4int count = 0;
1358	// ------ end debug
1359	while(!CheckPauliPrinciple(products))
1360	{
1361	// ------ debug
1362	count++;
1363	// ------ end debug
1364	ClearAndDestroy(products);
1365	if(!absorber.FindProducts(*kt))
1366	throw G4HadronicException(__FILE__, __LINE__,
1367	"G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1368	}
1369	// ------ debug
1370	// G4cerr << "Absorb CheckPauliPrinciple count= " << count << G4endl;
1371	// ------ end debug
1372	G4KineticTrackVector toRemove; // build a vector for UpdateTrack...
1373	toRemove.push_back(kt);
1374	toDelete.push_back(kt); // delete the track later
1375	G4KineticTrackVector * absorbers = absorber.GetAbsorbers();
1376	UpdateTracksAndCollisions(&toRemove, absorbers, products);
1377	ClearAndDestroy(absorbers);
1378	}
1379	ClearAndDestroy(&toDelete);
1380	return true;
1381	}
1382
1383
1384
1385	// Capture all p and n with Energy < theCutOnP
1386	//----------------------------------------------------------------------------
1387	G4bool G4BinaryCascade::Capture(G4bool verbose)
1388	//----------------------------------------------------------------------------
1389	{
1390	G4KineticTrackVector captured;
1391	G4bool capture = false;
1392	std::vector<G4KineticTrack *>::iterator i;
1393
1394	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
1395
1396	G4double capturedEnergy = 0;
1397	G4int particlesAboveCut=0;
1398	G4int particlesBelowCut=0;
1399	if ( verbose ) G4cout << " Capture: secondaries " << theSecondaryList.size() << G4endl;
1400	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1401	{
1402	G4KineticTrack * kt = *i;
1403	if (verbose) G4cout << "Capture position, radius, state " <<kt->GetPosition().mag()<<" "<<theOuterRadius<<" "<<kt->GetState()<<G4endl;
1404	if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1405	{
1406	if((kt->GetDefinition() == G4Proton::Proton()) \|\|
1407	(kt->GetDefinition() == G4Neutron::Neutron()))
1408	{
1409	//GF cut on kinetic energy if(kt->Get4Momentum().vect().mag() >= theCutOnP)
1410	G4double field=RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
1411	-RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
1412	G4double energy= kt->Get4Momentum().e() - kt->GetActualMass() + field;
1413	if (verbose ) G4cout << "Capture: .e(), mass, field, energy" << kt->Get4Momentum().e() <<" "<<kt->GetActualMass()<<" "<<field<<" "<<energy<< G4endl;
1414	// if( energy < theCutOnP )
1415	// {
1416	capturedEnergy+=energy;
1417	++particlesBelowCut;
1418	// } else
1419	// {
1420	// ++particlesAboveCut;
1421	// }
1422	}
1423	}
1424	}
1425	if (verbose) G4cout << "Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
1426	<< particlesAboveCut << " " << particlesBelowCut << " " << capturedEnergy
1427	<< " " << capturedEnergy/particlesBelowCut << " " << 0.2*theCutOnP << G4endl;
1428	// if(particlesAboveCut==0 && particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1429	if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1430	{
1431	capture=true;
1432	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1433	{
1434	G4KineticTrack * kt = *i;
1435	if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1436	{
1437	if((kt->GetDefinition() == G4Proton::Proton()) \|\|
1438	(kt->GetDefinition() == G4Neutron::Neutron()))
1439	{
1440	captured.push_back(kt);
1441	kt->Hit(); //
1442	theCapturedList.push_back(kt);
1443	}
1444	}
1445	}
1446	UpdateTracksAndCollisions(&captured, NULL, NULL);
1447	}
1448
1449	return capture;
1450	}
1451
1452
1453	//----------------------------------------------------------------------------
1454	G4bool G4BinaryCascade::CheckPauliPrinciple(G4KineticTrackVector * products)
1455	//----------------------------------------------------------------------------
1456	{
1457	G4int A = the3DNucleus->GetMassNumber();
1458	G4int Z = the3DNucleus->GetCharge();
1459
1460	G4FermiMomentum fermiMom;
1461	fermiMom.Init(A, Z);
1462
1463	const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
1464
1465	G4KineticTrackVector::iterator i;
1466	G4ParticleDefinition * definition;
1467
1468	// ------ debug
1469	G4bool myflag = true;
1470	// ------ end debug
1471	// G4ThreeVector xpos(0);
1472	for(i = products->begin(); i != products->end(); ++i)
1473	{
1474	definition = (*i)->GetDefinition();
1475	if((definition == G4Proton::Proton()) \|\|
1476	(definition == G4Neutron::Neutron()))
1477	{
1478	G4ThreeVector pos = (*i)->GetPosition();
1479	G4double d = density->GetDensity(pos);
1480	// energy correspondiing to fermi momentum
1481	G4double eFermi = std::sqrt( sqr(fermiMom.GetFermiMomentum(d)) + (*i)->Get4Momentum().mag2() );
1482	if( definition == G4Proton::Proton() )
1483	{
1484	eFermi -= the3DNucleus->CoulombBarrier();
1485	}
1486	G4LorentzVector mom = (*i)->Get4Momentum();
1487	// ------ debug
1488	/*
1489	* G4cout << "p:[" << (1/MeV)mom.x() << " " << (1/MeV)mom.y() << " "
1490	* << (1/MeV)*mom.z() << "] \|p3\|:"
1491	* << (1/MeV)mom.vect().mag() << " E:" << (1/MeV)mom.t() << " ] m: "
1492	* << (1/MeV)*mom.mag() << " pos["
1493	* << (1/fermi)pos.x() << " "<< (1/fermi)pos.y() << " "
1494	* << (1/fermi)*pos.z() << "] \|Dpos\|: "
1495	* << (1/fermi)*(pos-xpos).mag() << " Pfermi:"
1496	* << (1/MeV)*p << G4endl;
1497	* xpos=pos;
1498	*/
1499	// ------ end debug
1500	if(mom.e() < eFermi )
1501	{
1502	// ------ debug
1503	myflag = false;
1504	// ------ end debug
1505	// return false;
1506	}
1507	}
1508	}
1509	#ifdef debug_BIC_CheckPauli
1510	if ( myflag )
1511	{
1512	for(i = products->begin(); i != products->end(); ++i)
1513	{
1514	definition = (*i)->GetDefinition();
1515	if((definition == G4Proton::Proton()) \|\|
1516	(definition == G4Neutron::Neutron()))
1517	{
1518	G4ThreeVector pos = (*i)->GetPosition();
1519	G4double d = density->GetDensity(pos);
1520	G4double pFermi = fermiMom.GetFermiMomentum(d);
1521	G4LorentzVector mom = (*i)->Get4Momentum();
1522	G4double field =((G4RKPropagation*)thePropagator)->GetField(definition->GetPDGEncoding(),pos);
1523	if ( mom.e()-mom.mag()+field > 160*MeV )
1524	{
1525	G4cout << "momentum problem pFermi=" << pFermi
1526	<< " mom, mom.m " << mom << " " << mom.mag()
1527	<< " field " << field << G4endl;
1528	}
1529	}
1530	}
1531	}
1532	#endif
1533
1534	return myflag;
1535	}
1536
1537	//----------------------------------------------------------------------------
1538	void G4BinaryCascade::StepParticlesOut()
1539	//----------------------------------------------------------------------------
1540	{
1541	G4int counter=0;
1542	G4int countreset=0;
1543	//G4cout << " nucl. Radius " << radius << G4endl;
1544	// G4cerr <<"pre-while- theSecondaryList "<<G4endl;
1545	while( theSecondaryList.size() > 0 )
1546	{
1547	G4int nsec=0;
1548	G4double minTimeStep = 1.e-12ns; // about 30fermi/(0.1c_light);1.e-12ns
1549	// i.e. a big step
1550	std::vector<G4KineticTrack *>::iterator i;
1551	for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1552	{
1553	G4KineticTrack * kt = *i;
1554	if( kt->GetState() == G4KineticTrack::inside )
1555	{
1556	nsec++;
1557	G4double tStep(0), tdummy(0);
1558	((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep);
1559	#ifdef debug_BIC_StepParticlesOut
1560	G4cout << " minTimeStep, tStep Particle " <<minTimeStep << " " <<tStep
1561	<< " " <<kt->GetDefinition()->GetParticleName()
1562	<< " 4mom " << kt->GetTrackingMomentum()<<G4endl;
1563	#endif
1564	if(tStep<minTimeStep && tStep> 0 )
1565	{
1566	minTimeStep = tStep;
1567	}
1568	} else if ( kt->GetState() != G4KineticTrack::outside ){
1569	PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
1570	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1571	}
1572	}
1573	minTimeStep *= 1.2;
1574	// G4cerr << "CaptureCount = "<<counter<<" "<<nsec<<" "<<minTimeStep<<" "<<1*ns<<G4endl;
1575	G4double timeToCollision=DBL_MAX;
1576	G4CollisionInitialState * nextCollision=0;
1577	if(theCollisionMgr->Entries() > 0)
1578	{
1579	nextCollision = theCollisionMgr->GetNextCollision();
1580	timeToCollision = nextCollision->GetCollisionTime()-theCurrentTime;
1581	G4cout << " NextCollision * , Time= " << nextCollision << " "
1582	<<timeToCollision<< G4endl;
1583	}
1584	if ( timeToCollision > minTimeStep )
1585	{
1586	DoTimeStep(minTimeStep);
1587	++counter;
1588	} else
1589	{
1590	if (!DoTimeStep(timeToCollision) )
1591	{
1592	// Check if nextCollision is still valid, ie. partcile did not leave nucleus
1593	if (theCollisionMgr->GetNextCollision() != nextCollision )
1594	{
1595	nextCollision = 0;
1596	}
1597	}
1598	// G4cerr <<"post- DoTimeStep 3"<<G4endl;
1599
1600	if(nextCollision)
1601	{
1602	if ( ApplyCollision(nextCollision))
1603	{
1604	// G4cout << "ApplyCollision sucess " << G4endl;
1605	} else
1606	{
1607	theCollisionMgr->RemoveCollision(nextCollision);
1608	}
1609	}
1610	}
1611
1612	if(countreset>100)
1613	{
1614	#ifdef debug_G4BinaryCascade
1615	G4cerr << "G4BinaryCascade.cc: Warning - aborting looping particle(s)" << G4endl;
1616	#endif
1617
1618	// add left secondaries to FinalSate
1619	std::vector<G4KineticTrack *>::iterator iter;
1620	for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
1621	{
1622	theFinalState.push_back(*iter);
1623	}
1624	theSecondaryList.clear();
1625
1626	break;
1627	}
1628
1629	if(Absorb())
1630	{
1631	// haveProducts = true;
1632	// G4cout << "Absorb sucess " << G4endl;
1633	}
1634
1635	if(Capture(false))
1636	{
1637	// haveProducts = true;
1638	#ifdef debug_BIC_StepParticlesOut
1639	G4cout << "Capture sucess " << G4endl;
1640	#endif
1641	}
1642	if ( counter > 100 && theCollisionMgr->Entries() == 0) // no collision, and stepping a while....
1643	{
1644	#ifdef debug_1_BinaryCascade
1645	PrintKTVector(&theSecondaryList,std::string("stepping 100 steps"));
1646	#endif
1647	FindCollisions(&theSecondaryList);
1648	counter=0;
1649	++countreset;
1650	}
1651	}
1652	// G4cerr <<"Finished capture loop "<<G4endl;
1653
1654	//G4cerr <<"pre- DoTimeStep 4"<<G4endl;
1655	DoTimeStep(DBL_MAX);
1656	//G4cerr <<"post- DoTimeStep 4"<<G4endl;
1657
1658
1659	}
1660
1661	//----------------------------------------------------------------------------
1662	void G4BinaryCascade::CorrectFinalPandE()
1663	//----------------------------------------------------------------------------
1664	//
1665	// Modify momenta of outgoing particles.
1666	// Assume two body decay, nucleus(@nominal mass) + sum of final state particles(SFSP).
1667	// momentum of SFSP shall be less than momentum for two body decay.
1668	//
1669	{
1670	#ifdef debug_BIC_CorrectFinalPandE
1671	G4cerr << "BIC: -CorrectFinalPandE called" << G4endl;
1672	#endif
1673
1674	if ( theFinalState.size() == 0 ) return;
1675
1676	G4KineticTrackVector::iterator i;
1677	G4LorentzVector pNucleus=GetFinal4Momentum();
1678	if ( pNucleus.e() == 0 ) return; // check against explicit 0 from GetNucleus4Momentum()
1679	#ifdef debug_BIC_CorrectFinalPandE
1680	G4cerr << " -CorrectFinalPandE 3" << G4endl;
1681	#endif
1682	G4LorentzVector pFinals(0);
1683	G4int nFinals(0);
1684	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1685	{
1686	pFinals += (*i)->Get4Momentum();
1687	++nFinals;
1688	#ifdef debug_BIC_CorrectFinalPandE
1689	G4cout <<"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
1690	<< " 4mom " << (*i)->Get4Momentum()<< G4endl;
1691	#endif
1692	}
1693	#ifdef debug_BIC_CorrectFinalPandE
1694	G4cout << "CorrectFinalPandE pN pF: " <<pNucleus << " " <<pFinals << G4endl;
1695	#endif
1696	G4LorentzVector pCM=pNucleus + pFinals;
1697
1698	G4LorentzRotation toCMS(-pCM.boostVector());
1699	pFinals *=toCMS;
1700
1701	#ifdef debug_BIC_CorrectFinalPandE
1702	G4cout << "CorrectFinalPandE pCM, CMS pCM " << pCM << " " <<toCMS*pCM<< G4endl;
1703	G4cout << "CorrectFinal CMS pN pF " <<toCMS*pNucleus << " "
1704	<<pFinals << G4endl
1705	<< " nucleus initial mass : " <<GetFinal4Momentum().mag()
1706	<<" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus << " " <<pNucleus.mag()<< " "
1707	<< pFinals.mag() << " " << pCM.mag() << G4endl;
1708	#endif
1709
1710	G4LorentzRotation toLab = toCMS.inverse();
1711
1712	G4double s = pCM.mag2();
1713	// G4double m10 = massInNucleus; //pNucleus.mag();
1714	G4double m10 = GetIonMass(currentZ,currentA);
1715	G4double m20 = pFinals.mag();
1716	if( s-(m10+m20)*(m10+m20) < 0 )
1717	{
1718	#ifdef debug_BIC_CorrectFinalPandE
1719	G4cout << "G4BinaryCascade::CorrectFinalPandE() : error! " << G4endl;
1720
1721	G4cout << "not enough mass to correct: mass, A,Z, mass(nucl), mass(finals) "
1722	<< std::sqrt(-s+(m10+m20)*(m10+m20)) << " "
1723	<< currentA << " " << currentZ << " "
1724	<< m10 << " " << m20
1725	<< G4endl;
1726	G4cerr << " -CorrectFinalPandE 4" << G4endl;
1727
1728	PrintKTVector(&theFinalState," mass problem");
1729	#endif
1730	return;
1731	}
1732
1733	// Three momentum in cm system
1734	G4double pInCM = std::sqrt((s-(m10+m20)(m10+m20))(s-(m10-m20)(m10-m20))/(4.s));
1735	#ifdef debug_BIC_CorrectFinalPandE
1736	G4cout <<" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM
1737	<< " " << (pFinals).vect().mag()<< " " << pInCM/(pFinals).vect().mag() << G4endl;
1738	#endif
1739	if ( pFinals.vect().mag() > pInCM )
1740	{
1741	G4ThreeVector p3finals=pInCM*pFinals.vect().unit();
1742
1743	// G4ThreeVector deltap=(p3finals - pFinals.vect() ) / nFinals;
1744	G4double factor=std::max(0.98,pInCM/pFinals.vect().mag()); // small correction
1745	G4LorentzVector qFinals(0);
1746	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1747	{
1748	// G4ThreeVector p3((toCMS(i)->Get4Momentum()).vect() + deltap);
1749	G4ThreeVector p3(factor(toCMS(*i)->Get4Momentum()).vect());
1750	G4LorentzVector p(p3,std::sqrt((*i)->Get4Momentum().mag2() + p3.mag2()));
1751	qFinals += p;
1752	p *= toLab;
1753	#ifdef debug_BIC_CorrectFinalPandE
1754	G4cout << " final p corrected: " << p << G4endl;
1755	#endif
1756	(*i)->Set4Momentum(p);
1757	}
1758	#ifdef debug_BIC_CorrectFinalPandE
1759	G4cout << "CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() << " "
1760	<<GetFinal4Momentum().mag() << G4endl
1761	<< " CMS pFinals , mag, 3.mag : " << qFinals << " " << qFinals.mag() << " " << qFinals.vect().mag()<< G4endl;
1762	G4cerr << " -CorrectFinalPandE 5 " << factor << G4endl;
1763	#endif
1764	}
1765	#ifdef debug_BIC_CorrectFinalPandE
1766	else { G4cerr << " -CorrectFinalPandE 6 - no correction done" << G4endl; }
1767	#endif
1768
1769	}
1770
1771	//----------------------------------------------------------------------------
1772	void G4BinaryCascade::UpdateTracksAndCollisions(
1773	//----------------------------------------------------------------------------
1774	G4KineticTrackVector * oldSecondaries,
1775	G4KineticTrackVector * oldTarget,
1776	G4KineticTrackVector * newSecondaries)
1777	{
1778	// G4cout << "Entering ... "<<oldTarget<<G4endl;
1779	std::vector<G4KineticTrack *>::iterator iter1, iter2;
1780
1781	// remove old secondaries from the secondary list
1782	if(oldSecondaries)
1783	{
1784	if(!oldSecondaries->empty())
1785	{
1786	for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
1787	++iter1)
1788	{
1789	iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(),
1790	*iter1);
1791	if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2);
1792	}
1793	theCollisionMgr->RemoveTracksCollisions(oldSecondaries);
1794	}
1795	}
1796
1797	// remove old target from the target list
1798	if(oldTarget)
1799	{
1800	// G4cout << "################## Debugging 0 "<<G4endl;
1801	if(oldTarget->size()!=0)
1802	{
1803
1804	// G4cout << "################## Debugging 1 "<<oldTarget->size()<<G4endl;
1805	for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
1806	{
1807	iter2 = std::find(theTargetList.begin(), theTargetList.end(),
1808	*iter1);
1809	theTargetList.erase(iter2);
1810	}
1811	theCollisionMgr->RemoveTracksCollisions(oldTarget);
1812	}
1813	}
1814
1815	if(newSecondaries)
1816	{
1817	if(!newSecondaries->empty())
1818	{
1819	// insert new secondaries in the secondary list
1820	for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
1821	++iter1)
1822	{
1823	theSecondaryList.push_back(*iter1);
1824	}
1825	// look for collisions of new secondaries
1826	FindCollisions(newSecondaries);
1827	}
1828	}
1829	// G4cout << "Exiting ... "<<oldTarget<<G4endl;
1830	}
1831
1832
1833	class SelectFromKTV
1834	{
1835	private:
1836	G4KineticTrackVector * ktv;
1837	G4KineticTrack::CascadeState wanted_state;
1838	public:
1839	SelectFromKTV(G4KineticTrackVector * out, G4KineticTrack::CascadeState astate)
1840	:
1841	ktv(out), wanted_state(astate)
1842	{};
1843	void operator() (G4KineticTrack *& kt) const
1844	{
1845	if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt);
1846	};
1847	};
1848
1849
1850
1851	//----------------------------------------------------------------------------
1852	G4bool G4BinaryCascade::DoTimeStep(G4double theTimeStep)
1853	//----------------------------------------------------------------------------
1854	{
1855
1856	#ifdef debug_BIC_DoTimeStep
1857	G4ping debug("debug_G4BinaryCascade");
1858	debug.push_back("======> DoTimeStep 1"); debug.dump();
1859	G4cerr <<"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
1860	<< " , time="<<theCurrentTime << G4endl;
1861	PrintKTVector(&theSecondaryList, std::string("DoTimeStep - theSecondaryList"));
1862	//PrintKTVector(&theTargetList, std::string("DoTimeStep - theTargetList"));
1863	#endif
1864
1865	G4bool success=true;
1866	std::vector<G4KineticTrack *>::iterator iter;
1867
1868	G4KineticTrackVector * kt_outside = new G4KineticTrackVector;
1869	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1870	SelectFromKTV(kt_outside,G4KineticTrack::outside));
1871	//PrintKTVector(kt_outside, std::string("DoTimeStep - found outside"));
1872
1873	G4KineticTrackVector * kt_inside = new G4KineticTrackVector;
1874	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1875	SelectFromKTV(kt_inside, G4KineticTrack::inside));
1876	// PrintKTVector(kt_inside, std::string("DoTimeStep - found inside"));
1877	//-----
1878	G4KineticTrackVector dummy; // needed for re-usability
1879	#ifdef debug_BIC_DoTimeStep
1880	G4cout << "NOW WE ARE ENTERING THE TRANSPORT"<<G4endl;
1881	#endif
1882
1883	// =================== Here we move the particles ===================
1884
1885	thePropagator->Transport(theSecondaryList, dummy, theTimeStep);
1886
1887	// =================== Here we move the particles ===================
1888
1889	//------
1890
1891	theMomentumTransfer += thePropagator->GetMomentumTransfer();
1892	#ifdef debug_BIC_DoTimeStep
1893	G4cout << "DoTimeStep : theMomentumTransfer = " << theMomentumTransfer << G4endl;
1894	PrintKTVector(&theSecondaryList, std::string("DoTimeStep - secondaries aft trsprt"));
1895	#endif
1896
1897	// Partclies which went INTO nucleus
1898
1899	G4KineticTrackVector * kt_gone_in = new G4KineticTrackVector;
1900	std::for_each( kt_outside->begin(),kt_outside->end(),
1901	SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
1902	// PrintKTVector(kt_gone_in, std::string("DoTimeStep - gone in"));
1903
1904
1905	// Partclies which went OUT OF nucleus
1906	G4KineticTrackVector * kt_gone_out = new G4KineticTrackVector;
1907	std::for_each( kt_inside->begin(),kt_inside->end(),
1908	SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
1909
1910	// PrintKTVector(kt_gone_out, std::string("DoTimeStep - gone out"));
1911
1912	G4KineticTrackVector *fail=CorrectBarionsOnBoundary(kt_gone_in,kt_gone_out);
1913
1914	if ( fail )
1915	{
1916	// some particle(s) supposed to enter/leave were miss_nucleus/captured by the correction
1917	kt_gone_in->clear();
1918	std::for_each( kt_outside->begin(),kt_outside->end(),
1919	SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
1920
1921	kt_gone_out->clear();
1922	std::for_each( kt_inside->begin(),kt_inside->end(),
1923	SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
1924
1925	#ifdef debug_BIC_DoTimeStep
1926	PrintKTVector(fail,std::string(" Failed to go in/out -> miss_nucleus/captured"));
1927	PrintKTVector(kt_gone_in, std::string("recreated kt_gone_in"));
1928	PrintKTVector(kt_gone_out, std::string("recreated kt_gone_out"));
1929	#endif
1930	delete fail;
1931	}
1932
1933	// Add tracks missing nucleus and tracks going straight though to addFinals
1934	std::for_each( kt_outside->begin(),kt_outside->end(),
1935	SelectFromKTV(kt_gone_out,G4KineticTrack::miss_nucleus));
1936	//PrintKTVector(kt_gone_out, std::string("miss to append to final state.."));
1937	std::for_each( kt_outside->begin(),kt_outside->end(),
1938	SelectFromKTV(kt_gone_out,G4KineticTrack::gone_out));
1939
1940	#ifdef debug_BIC_DoTimeStep
1941	PrintKTVector(kt_gone_out, std::string("append to final state.."));
1942	#endif
1943
1944	theFinalState.insert(theFinalState.end(),
1945	kt_gone_out->begin(),kt_gone_out->end());
1946
1947	// Partclies which could not leave nucleus, captured...
1948	G4KineticTrackVector * kt_captured = new G4KineticTrackVector;
1949	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
1950	SelectFromKTV(kt_captured, G4KineticTrack::captured));
1951
1952	// Check no track is part in next collision, ie.
1953	// this step was to far, and collisions should not occur any more
1954
1955	if ( theCollisionMgr->Entries()> 0 )
1956	{
1957	if (kt_gone_out->size() )
1958	{
1959	G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
1960	iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
1961	if ( iter != kt_gone_out->end() )
1962	{
1963	success=false;
1964	#ifdef debug_BIC_DoTimeStep
1965	G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
1966	#endif
1967	}
1968	}
1969	if ( kt_captured->size() )
1970	{
1971	G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
1972	iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
1973	if ( iter != kt_captured->end() )
1974	{
1975	success=false;
1976	#ifdef debug_BIC_DoTimeStep
1977	G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
1978	#endif
1979	}
1980	}
1981
1982	}
1983	// PrintKTVector(kt_gone_out," kt_gone_out be4 updatetrack...");
1984	UpdateTracksAndCollisions(kt_gone_out,0 ,0);
1985
1986
1987	if ( kt_captured->size() )
1988	{
1989	theCapturedList.insert(theCapturedList.end(),
1990	kt_captured->begin(),kt_captured->end());
1991	//should be std::for_each(kt_captured->begin(),kt_captured->end(),
1992	// std::mem_fun(&G4KineticTrack::Hit));
1993	// but VC 6 requires:
1994	std::vector<G4KineticTrack *>::iterator i_captured;
1995	for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
1996	{
1997	(*i_captured)->Hit();
1998	}
1999	// PrintKTVector(kt_captured," kt_captured be4 updatetrack...");
2000	UpdateTracksAndCollisions(kt_captured, NULL, NULL);
2001	}
2002
2003	#ifdef debug_G4BinaryCascade
2004	delete kt_inside;
2005	kt_inside = new G4KineticTrackVector;
2006	std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2007	SelectFromKTV(kt_inside, G4KineticTrack::inside));
2008	if ( currentZ != (GetTotalCharge(theTargetList)
2009	+ GetTotalCharge(theCapturedList)
2010	+ GetTotalCharge(*kt_inside)) )
2011	{
2012	G4cout << " error-DoTimeStep aft, A, Z: " << currentA << " " << currentZ
2013	<< " sum(tgt,capt,active) "
2014	<< GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside)
2015	<< " targets: " << GetTotalCharge(theTargetList)
2016	<< " captured: " << GetTotalCharge(theCapturedList)
2017	<< " active: " << GetTotalCharge(*kt_inside)
2018	<< G4endl;
2019	}
2020	#endif
2021
2022	delete kt_inside;
2023	delete kt_outside;
2024	delete kt_captured;
2025	delete kt_gone_in;
2026	delete kt_gone_out;
2027
2028	// G4cerr <<"G4BinaryCascade::DoTimeStep: exit "<<G4endl;
2029	theCurrentTime += theTimeStep;
2030
2031	//debug.push_back("======> DoTimeStep 2"); debug.dump();
2032	return success;
2033
2034	}
2035
2036	//----------------------------------------------------------------------------
2037	G4KineticTrackVector* G4BinaryCascade::CorrectBarionsOnBoundary(
2038	G4KineticTrackVector *in,
2039	G4KineticTrackVector *out)
2040	//----------------------------------------------------------------------------
2041	{
2042	G4KineticTrackVector * kt_fail(0);
2043	std::vector<G4KineticTrack *>::iterator iter;
2044	// G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2045	// << currentZ << " "<< currentA << G4endl;
2046	if (in->size())
2047	{
2048	G4int secondaries_in(0);
2049	G4int secondaryBarions_in(0);
2050	G4int secondaryCharge_in(0);
2051	G4double secondaryMass_in(0);
2052
2053	for ( iter =in->begin(); iter != in->end(); ++iter)
2054	{
2055	++secondaries_in;
2056	secondaryCharge_in += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2057	if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
2058	{
2059	secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
2060	if((*iter)->GetDefinition() == G4Neutron::Neutron() \|\|
2061	(*iter)->GetDefinition() == G4Proton::Proton() )
2062	{
2063	secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
2064	} else {
2065	secondaryMass_in += G4Proton::Proton()->GetPDGMass();
2066	}
2067	}
2068	}
2069	G4double mass_initial= GetIonMass(currentZ,currentA);
2070
2071	currentZ += secondaryCharge_in;
2072	currentA += secondaryBarions_in;
2073
2074	// G4cout << "CorrectBarionsOnBoundary,secondaryCharge_in, secondaryBarions_in "
2075	// << secondaryCharge_in << " "<< secondaryBarions_in << G4endl;
2076
2077	G4double mass_final= GetIonMass(currentZ,currentA);
2078
2079	G4double correction= secondaryMass_in + mass_initial - mass_final;
2080	if (secondaries_in>1)
2081	{correction /= secondaries_in;}
2082
2083	#ifdef debug_BIC_CorrectBarionsOnBoundary
2084	G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2085	<< "secondaryCharge_in,secondaryBarions_in,"
2086	<< "energy correction,m_secondry,m_nucl_init,m_nucl_final "
2087	<< currentZ << " "<< currentA <<" "
2088	<< secondaryCharge_in<<" "<<secondaryBarions_in<<" "
2089	<< correction << " "
2090	<< secondaryMass_in << " "
2091	<< mass_initial << " "
2092	<< mass_final << " "
2093	<< G4endl;
2094	PrintKTVector(in,std::string("in be4 correction"));
2095	#endif
2096
2097	for ( iter = in->begin(); iter != in->end(); ++iter)
2098	{
2099	if ((iter)->GetTrackingMomentum().e()+correction > (iter)->GetActualMass())
2100	{
2101	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + correction);
2102	} else {
2103	//particle cannot go in, put to miss_nucleus
2104	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2105	(*iter)->SetState(G4KineticTrack::miss_nucleus);
2106	// Undo correction for Colomb Barrier
2107	G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2108	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + barrier);
2109	if ( ! kt_fail ) kt_fail=new G4KineticTrackVector;
2110	kt_fail->push_back(*iter);
2111	currentZ -= G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2112	currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
2113
2114	}
2115
2116	}
2117	#ifdef debug_BIC_CorrectBarionsOnBoundary
2118	G4cout << " CorrectBarionsOnBoundary, aft, A, Z, sec-Z,A,m,m_in_nucleus "
2119	<< currentA << " " << currentZ << " "
2120	<< secondaryCharge_in << " " << secondaryBarions_in << " "
2121	<< secondaryMass_in << " "
2122	<< G4endl;
2123	PrintKTVector(in,std::string("in AFT correction"));
2124	#endif
2125
2126	}
2127	//----------------------------------------------
2128	if (out->size())
2129	{
2130	G4int secondaries_out(0);
2131	G4int secondaryBarions_out(0);
2132	G4int secondaryCharge_out(0);
2133	G4double secondaryMass_out(0);
2134
2135	for ( iter =out->begin(); iter != out->end(); ++iter)
2136	{
2137	++secondaries_out;
2138	secondaryCharge_out += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2139	if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
2140	{
2141	secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
2142	if((*iter)->GetDefinition() == G4Neutron::Neutron() \|\|
2143	(*iter)->GetDefinition() == G4Proton::Proton() )
2144	{
2145	secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
2146	} else {
2147	secondaryMass_out += G4Neutron::Neutron()->GetPDGMass();
2148	}
2149	}
2150	}
2151
2152	G4double mass_initial= GetIonMass(currentZ,currentA);
2153	currentA -=secondaryBarions_out;
2154	currentZ -=secondaryCharge_out;
2155
2156	// G4cout << "CorrectBarionsOnBoundary,secondaryCharge_out, secondaryBarions_out"
2157	// << secondaryCharge_out << " "<< secondaryBarions_out << G4endl;
2158
2159	// a delta minus will do currentZ < 0 in light nuclei
2160	// if (currentA < 0 \|\| currentZ < 0 )
2161	if (currentA < 0 )
2162	{
2163	G4cerr << "G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
2164	secondaryBarions_out << " " << secondaryCharge_out << G4endl;
2165	PrintKTVector(&theTargetList,"CorrectBarionsOnBoundary Target");
2166	PrintKTVector(&theCapturedList,"CorrectBarionsOnBoundary Captured");
2167	PrintKTVector(&theSecondaryList,"CorrectBarionsOnBoundary Secondaries");
2168	G4cerr << "G4BinaryCascade - currentA, currentZ " << currentA << " " << currentZ << G4endl;
2169	throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
2170	}
2171	G4double mass_final=GetIonMass(currentZ,currentA);
2172	G4double correction= mass_initial - mass_final - secondaryMass_out;
2173
2174	if (secondaries_out>1) correction /= secondaries_out;
2175	#ifdef debug_BIC_CorrectBarionsOnBoundary
2176	G4cout << "DoTimeStep,currentZ,currentA,"
2177	<< "secondaries_out,"
2178	<<"secondaryCharge_out,secondaryBarions_out,"
2179	<<"energy correction,m_secondry,m_nucl_init,m_nucl_final "
2180	<< " "<< currentZ << " "<< currentA <<" "
2181	<< secondaries_out << " "
2182	<< secondaryCharge_out<<" "<<secondaryBarions_out<<" "
2183	<< correction << " "
2184	<< secondaryMass_out << " "
2185	<< mass_initial << " "
2186	<< mass_final << " "
2187	<< G4endl;
2188	PrintKTVector(out,std::string("out be4 correction"));
2189	#endif
2190
2191	for ( iter = out->begin(); iter != out->end(); ++iter)
2192	{
2193	if ((iter)->GetTrackingMomentum().e()+correction > (iter)->GetActualMass())
2194	{
2195	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() + correction);
2196	} else
2197	{
2198	// particle cannot go out due to change of nuclear potential!
2199	// capture protons and neutrons;
2200	if(((*iter)->GetDefinition() == G4Proton::Proton()) \|\|
2201	((*iter)->GetDefinition() == G4Neutron::Neutron()))
2202	{
2203	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2204	(*iter)->SetState(G4KineticTrack::captured);
2205	// Undo correction for Colomb Barrier
2206	G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2207	(iter)->UpdateTrackingMomentum((iter)->GetTrackingMomentum().e() - barrier);
2208	if ( kt_fail == 0 ) kt_fail=new G4KineticTrackVector;
2209	kt_fail->push_back(*iter);
2210	currentZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge());
2211	currentA += (*iter)->GetDefinition()->GetBaryonNumber();
2212	}
2213	#ifdef debug_BIC_CorrectBarionsOnBoundary
2214	else
2215	{
2216	G4cout << "Not correcting outgoing " << *iter << " "
2217	<< (*iter)->GetDefinition()->GetPDGEncoding() << " "
2218	<< (*iter)->GetDefinition()->GetParticleName() << G4endl;
2219	PrintKTVector(out,std::string("outgoing, one not corrected"));
2220	}
2221	#endif
2222	}
2223	}
2224
2225	#ifdef debug_BIC_CorrectBarionsOnBoundary
2226	PrintKTVector(out,std::string("out AFTER correction"));
2227	G4cout << " DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
2228	<< currentA << " "<< currentZ << " "
2229	<< secondaryCharge_out << " "<< secondaryBarions_out << " "<<
2230	secondaryMass_out << " "
2231	<< massInNucleus << " "
2232	<< G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2233	<< " " << massInNucleus -G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2234	<< G4endl;
2235	#endif
2236	}
2237
2238	return kt_fail;
2239	}
2240
2241
2242	//----------------------------------------------------------------------------
2243
2244	G4Fragment * G4BinaryCascade::FindFragments()
2245	//----------------------------------------------------------------------------
2246	{
2247
2248	G4int a = theTargetList.size()+theCapturedList.size();
2249	#ifdef debug_BIC_FindFragments
2250	G4cout << "target, captured, secondary: "
2251	<< theTargetList.size() << " "
2252	<< theCapturedList.size()<< " "
2253	<< theSecondaryList.size()
2254	<< G4endl;
2255	#endif
2256	G4int zTarget = 0;
2257	G4KineticTrackVector::iterator i;
2258	for(i = theTargetList.begin(); i != theTargetList.end(); ++i)
2259	{
2260	if((*i)->GetDefinition()->GetPDGCharge() == eplus)
2261	{
2262	zTarget++;
2263	}
2264	}
2265
2266	G4int zCaptured = 0;
2267	G4LorentzVector CapturedMomentum=0;
2268	for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2269	{
2270	CapturedMomentum += (*i)->Get4Momentum();
2271	if((*i)->GetDefinition()->GetPDGCharge() == eplus)
2272	{
2273	zCaptured++;
2274	}
2275	}
2276
2277	G4int z = zTarget+zCaptured;
2278
2279	#ifdef debug_G4BinaryCascade
2280	if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) )
2281	{
2282	G4cout << " FindFragment Counting error z a " << z << " " <<a << " "
2283	<< GetTotalCharge(theTargetList) << " " << GetTotalCharge(theCapturedList)<<
2284	G4endl;
2285	PrintKTVector(&theTargetList, std::string("theTargetList"));
2286	PrintKTVector(&theCapturedList, std::string("theCapturedList"));
2287	}
2288	#endif
2289	//debug
2290	/*
2291	* G4cout << " Fragment mass table / real "
2292	* << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
2293	* << " / " << GetFinal4Momentum().mag()
2294	* << " difference "
2295	* << GetFinal4Momentum().mag() - G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
2296	* << G4endl;
2297	*/
2298	//
2299	// if(getenv("BCDEBUG") ) G4cerr << "Fragment A, Z "<< a <<" "<< z<<G4endl;
2300	if ( z < 1 ) return 0;
2301
2302	G4int holes = the3DNucleus->GetMassNumber() - theTargetList.size();
2303	G4int excitons = theCapturedList.size();
2304	#ifdef debug_BIC_FindFragments
2305	G4cout << "Fragment: a= " << a
2306	<< " z= " << z
2307	<< " particles= " << excitons
2308	<< " Charged= " << zCaptured
2309	<< " holes= " << holes
2310	<< " excitE= " <<GetExcitationEnergy()
2311	<< " Final4Momentum= " << GetFinalNucleusMomentum()
2312	<< " capturMomentum= " << CapturedMomentum
2313	<< G4endl;
2314	#endif
2315
2316	G4Fragment * fragment = new G4Fragment(a,z,GetFinalNucleusMomentum());
2317	fragment->SetNumberOfHoles(holes);
2318
2319	//GF fragment->SetNumberOfParticles(excitons-holes);
2320	fragment->SetNumberOfParticles(excitons);
2321	fragment->SetNumberOfCharged(zCaptured);
2322	G4ParticleDefinition * aIonDefinition =
2323	G4ParticleTable::GetParticleTable()->FindIon(a,z,0,z);
2324	fragment->SetParticleDefinition(aIonDefinition);
2325
2326	return fragment;
2327	}
2328
2329	//----------------------------------------------------------------------------
2330
2331	G4LorentzVector G4BinaryCascade::GetFinal4Momentum()
2332	//----------------------------------------------------------------------------
2333	{
2334	// the initial 3-momentum will differ from 0, if nucleus created by string model.
2335	G4LorentzVector final4Momentum = theInitial4Mom;
2336	G4KineticTrackVector::iterator i;
2337	for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
2338	{
2339	final4Momentum += (*i)->GetTrackingMomentum();
2340	//G4cout << "Initial state: "<<(*i)->Get4Momentum()<<G4endl;
2341	}
2342	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2343	{
2344	final4Momentum -= (*i)->Get4Momentum();
2345	}
2346
2347	if((final4Momentum.vect()/final4Momentum.e()).mag()>1.0 && currentA > 0)
2348	{
2349	# ifdef debug_BIC_Final4Momentum
2350	G4cerr << G4endl;
2351	G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl;
2352	G4KineticTrackVector::iterator i;
2353	G4cerr <<" Initial nucleus "<<theInitial4Mom<<G4endl;
2354	for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
2355	{
2356	G4cerr << " Initial state (get4M), (trackingM): "
2357	<<(i)->Get4Momentum()<<", " << (i)->GetTrackingMomentum() <<","
2358	<<(*i)->GetDefinition()->GetParticleName()<<G4endl;
2359	}
2360	for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2361	{
2362	G4cerr <<" Final state: "<<(i)->Get4Momentum()<<(i)->GetDefinition()->GetParticleName()<<G4endl;
2363	}
2364	G4cerr<< " Final4Momentum = "<<final4Momentum <<" "<<final4Momentum.m()<<G4endl;
2365	G4cerr <<" current A, Z = "<< currentA<<", "<<currentZ<<G4endl;
2366	G4cerr << G4endl;
2367	# endif
2368
2369	final4Momentum=G4LorentzVector(0);
2370	}
2371	return final4Momentum;
2372	}
2373
2374	//----------------------------------------------------------------------------
2375	G4LorentzVector G4BinaryCascade::GetFinalNucleusMomentum()
2376	//----------------------------------------------------------------------------
2377	{
2378	// return momentum of nucleus for use with precompound model; also keep transformation to
2379	// apply to precompoud products.
2380
2381	G4LorentzVector CapturedMomentum=0;
2382	G4KineticTrackVector::iterator i;
2383	// G4cout << "GetFinalNucleusMomentum Captured size: " <<theCapturedList.size() << G4endl;
2384	for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2385	{
2386	CapturedMomentum += (*i)->Get4Momentum();
2387	}
2388	//G4cout << "GetFinalNucleusMomentum CapturedMomentum= " <<CapturedMomentum << G4endl;
2389	// G4cerr << "it 9"<<G4endl;
2390
2391	G4LorentzVector NucleusMomentum = GetFinal4Momentum();
2392	if ( NucleusMomentum.e() > 0 )
2393	{
2394	// G4cout << "GetFinalNucleusMomentum GetFinal4Momentum= " <<NucleusMomentum <<" "<<NucleusMomentum.mag()<<G4endl;
2395	// boost nucleus to a frame such that the momentum of nucleus == momentum of Captured
2396	G4ThreeVector boost= (NucleusMomentum.vect() -CapturedMomentum.vect())/NucleusMomentum.e();
2397	if(boost.mag2()>1.0)
2398	{
2399	# ifdef debug_BIC_FinalNucleusMomentum
2400	G4cerr << "G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<G4endl;
2401	G4cerr << "it 0"<<boost <<G4endl;
2402	G4cerr << "it 01"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2403	G4cout <<" testing boost "<<boost<<" "<<boost.mag()<<G4endl;
2404	# endif
2405	boost=G4ThreeVector(0);
2406	NucleusMomentum=G4LorentzVector(0);
2407	}
2408	G4LorentzRotation nucleusBoost( -boost );
2409	precompoundLorentzboost.set( boost );
2410	#ifdef debug_debug_BIC_FinalNucleusMomentum
2411	G4cout << "GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2412	#endif
2413	NucleusMomentum *= nucleusBoost;
2414	#ifdef debug_BIC_FinalNucleusMomentum
2415	G4cout << "GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<G4endl;
2416	#endif
2417	}
2418	return NucleusMomentum;
2419	}
2420
2421	//----------------------------------------------------------------------------
2422	G4ReactionProductVector * G4BinaryCascade::Propagate1H1(
2423	//----------------------------------------------------------------------------
2424	G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus)
2425	{
2426	G4ReactionProductVector * products = new G4ReactionProductVector;
2427	G4ParticleDefinition * aHTarg = G4Proton::ProtonDefinition();
2428	G4double mass = aHTarg->GetPDGMass();
2429	if (nucleus->GetCharge() == 0) aHTarg = G4Neutron::NeutronDefinition();
2430	mass = aHTarg->GetPDGMass();
2431	G4KineticTrackVector * secs = 0;
2432	G4ThreeVector pos(0,0,0);
2433	G4LorentzVector mom(mass);
2434	G4KineticTrack aTarget(aHTarg, 0., pos, mom);
2435	G4bool done(false);
2436	std::vector<G4KineticTrack *>::iterator iter, jter;
2437	// data member static G4Scatterer theH1Scatterer;
2438	//G4cout << " start 1H1 for " << (*secondaries).front()->GetDefinition()->GetParticleName()
2439	// << " on " << aHTarg->GetParticleName() << G4endl;
2440	G4int tryCount(0);
2441	while(!done && tryCount++ <200)
2442	{
2443	if(secs)
2444	{
2445	std::for_each(secs->begin(), secs->end(), DeleteKineticTrack());
2446	delete secs;
2447	}
2448	secs = theH1Scatterer->Scatter((secondaries).front(), aTarget);
2449	for(size_t ss=0; secs && ss<secs->size(); ss++)
2450	{
2451	// G4cout << "1H1 " << (*secs)[ss]->GetDefinition()->GetParticleName()
2452	// << ", shortlived? "<< (*secs)[ss]->GetDefinition()->IsShortLived()<< G4endl;
2453	if((*secs)[ss]->GetDefinition()->IsShortLived()) done = true;
2454	}
2455	// G4cout << G4endl;
2456	}
2457	size_t current(0);
2458	for(current=0; current<secs->size(); current++)
2459	{
2460	if((*secs)[current]->GetDefinition()->IsShortLived())
2461	{
2462	G4KineticTrackVector * dec = (*secs)[current]->Decay();
2463	for(jter=dec->begin(); jter != dec->end(); jter++)
2464	{
2465	//G4cout << "Decay"<<G4endl;
2466	secs->push_back(*jter);
2467	//G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<<G4endl;
2468	}
2469	delete (*secs)[current];
2470	delete dec;
2471	}
2472	else
2473	{
2474	//G4cout << "FS "<<G4endl;
2475	//G4cout << "FS "<<(*secs)[current]->GetDefinition()->GetParticleName()<<G4endl;
2476	theFinalState.push_back((*secs)[current]);
2477	}
2478	}
2479	//G4cout << "Through loop"<<G4endl;
2480	delete secs;
2481	for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2482	{
2483	G4KineticTrack * kt = *iter;
2484	G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
2485	aNew->SetMomentum(kt->Get4Momentum().vect());
2486	aNew->SetTotalEnergy(kt->Get4Momentum().e());
2487	products->push_back(aNew);
2488	#ifdef debug_1_BinaryCascade
2489	if (! kt->GetDefinition()->GetPDGStable() )
2490	{
2491	if (kt->GetDefinition()->IsShortLived())
2492	{
2493	G4cout << "final shortlived : ";
2494	} else
2495	{
2496	G4cout << "final un stable : ";
2497	}
2498	G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
2499	}
2500	#endif
2501	delete kt;
2502	}
2503	theFinalState.clear();
2504	return products;
2505
2506	}
2507
2508	//----------------------------------------------------------------------------
2509	G4ThreeVector G4BinaryCascade::GetSpherePoint(
2510	G4double r, const G4LorentzVector & mom4)
2511	//----------------------------------------------------------------------------
2512	{
2513	// Get a point outside radius.
2514	// point is random in plane (circle of radius r) orthogonal to mom,
2515	// plus -1rmom->vect()->unit();
2516	G4ThreeVector o1, o2;
2517	G4ThreeVector mom = mom4.vect();
2518
2519	o1= mom.orthogonal(); // we simply need any vector non parallel
2520	o2= mom.cross(o1); // o2 is now orthogonal to mom and o1, ie. o1 and o2 define plane.
2521
2522	G4double x2, x1;
2523
2524	do
2525	{
2526	x1=(G4UniformRand()-.5)*2;
2527	x2=(G4UniformRand()-.5)*2;
2528	} while (sqr(x1) +sqr(x2) > 1.);
2529
2530	return G4ThreeVector(r(x1o1.unit() + x2o2.unit() - 1.5 mom.unit()));
2531
2532
2533
2534	/*
2535	* // Get a point uniformly distributed on the surface of a sphere,
2536	* // with z < 0.
2537	* G4double b = r*G4UniformRand(); // impact parameter
2538	* G4double phi = G4UniformRand()2pi;
2539	* G4double x = b*std::cos(phi);
2540	* G4double y = b*std::sin(phi);
2541	* G4double z = -std::sqrt(rr-bb);
2542	* z *= 1.001; // Get position a little bit out of the sphere...
2543	* point.setX(x);
2544	* point.setY(y);
2545	* point.setZ(z);
2546	*/
2547	}
2548
2549	//----------------------------------------------------------------------------
2550
2551	void G4BinaryCascade::ClearAndDestroy(G4KineticTrackVector * ktv)
2552	//----------------------------------------------------------------------------
2553	{
2554	std::vector<G4KineticTrack *>::iterator i;
2555	for(i = ktv->begin(); i != ktv->end(); ++i)
2556	delete (*i);
2557	ktv->clear();
2558	}
2559
2560	//----------------------------------------------------------------------------
2561	void G4BinaryCascade::ClearAndDestroy(G4ReactionProductVector * rpv)
2562	//----------------------------------------------------------------------------
2563	{
2564	std::vector<G4ReactionProduct *>::iterator i;
2565	for(i = rpv->begin(); i != rpv->end(); ++i)
2566	delete (*i);
2567	rpv->clear();
2568	}
2569
2570	//----------------------------------------------------------------------------
2571	void G4BinaryCascade::PrintKTVector(G4KineticTrackVector * ktv, std::string comment)
2572	//----------------------------------------------------------------------------
2573	{
2574	if (comment.size() > 0 ) G4cout << comment << G4endl;
2575	G4cout << " vector: " << ktv << ", number of tracks: " << ktv->size()
2576	<< G4endl;
2577	std::vector<G4KineticTrack *>::iterator i;
2578	G4int count;
2579
2580	for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
2581	{
2582	G4KineticTrack * kt = *i;
2583	G4cout << " track n. " << count;
2584	PrintKTVector(kt);
2585	}
2586	}
2587	//----------------------------------------------------------------------------
2588	void G4BinaryCascade::PrintKTVector(G4KineticTrack * kt, std::string comment)
2589	//----------------------------------------------------------------------------
2590	{
2591	if (comment.size() > 0 ) G4cout << comment << G4endl;
2592	G4cout << ", id: " << kt << G4endl;
2593	G4ThreeVector pos = kt->GetPosition();
2594	G4LorentzVector mom = kt->Get4Momentum();
2595	G4LorentzVector tmom = kt->GetTrackingMomentum();
2596	G4ParticleDefinition * definition = kt->GetDefinition();
2597	G4cout << " definition: " << definition->GetPDGEncoding() << " pos: "
2598	<< 1/fermipos << " R: " << 1/fermipos.mag() << " 4mom: "
2599	<< 1/MeVmom <<"Tr_mom" << 1/MeVtmom << " P: " << 1/MeV*mom.vect().mag()
2600	<< " M: " << 1/MeV*mom.mag() << G4endl;
2601	G4cout <<" trackstatus: "<<kt->GetState()<<G4endl;
2602	}
2603
2604	//----------------------------------------------------------------------------
2605	G4bool G4BinaryCascade::CheckDecay(G4KineticTrackVector * products)
2606	//----------------------------------------------------------------------------
2607	{
2608	G4int A = the3DNucleus->GetMassNumber();
2609	G4int Z = the3DNucleus->GetCharge();
2610
2611	G4FermiMomentum fermiMom;
2612	fermiMom.Init(A, Z);
2613
2614	const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
2615
2616	G4KineticTrackVector::iterator i;
2617	G4ParticleDefinition * definition;
2618
2619	// ------ debug
2620	G4bool myflag = true;
2621	// ------ end debug
2622	for(i = products->begin(); i != products->end(); ++i)
2623	{
2624	definition = (*i)->GetDefinition();
2625	if((definition->GetParticleName() != "delta++" )&&
2626	(definition->GetParticleName() != "delta+" )&&
2627	(definition->GetParticleName() != "delta0" )&&
2628	(definition->GetParticleName() != "delta-" ))
2629	continue;
2630	G4ThreeVector pos = (*i)->GetPosition();
2631	G4double d = density->GetDensity(pos);
2632	G4double pFermi= fermiMom.GetFermiMomentum(d);
2633	G4LorentzVector mom = (*i)->Get4Momentum();
2634	G4LorentzRotation boost(mom.boostVector());
2635	G4ThreeVector pion3(227MeV mom.vect().unit()); // 227 is decay product in rest frame
2636	G4LorentzVector pion(pion3, std::sqrt(sqr(140*MeV) +pion3.mag()));
2637	// G4cout << "pi rest " << pion << G4endl;
2638	pion = boost * pion;
2639	// G4cout << "pi lab " << pion << G4endl;
2640	// ------ debug
2641	// G4cout << "p:[" << (1/MeV)pion.x() << " " << (1/MeV)pion.y() << " "
2642	// << (1/MeV)*pion.z() << "] \|p3\|:"
2643	// << (1/MeV)pion.vect().mag() << " E:" << (1/MeV)pion.t() << " ] m: "
2644	// << (1/MeV)*pion.mag() << " pos["
2645	// << (1/fermi)pos.x() << " "<< (1/fermi)pos.y() << " "
2646	// << (1/fermi)*pos.z() << "] \|Dpos\|: "
2647	// << " Pfermi:"
2648	// << (1/MeV)*pFermi << G4endl;
2649	// ------ end debug
2650
2651	if( pion.vect().mag() < pFermi )
2652	{
2653	// ------ debug
2654	myflag = false;
2655	// ------ end debug
2656	}
2657	}
2658	return myflag;
2659	// return true;
2660	}
2661
2662	//----------------------------------------------------------------------------
2663	G4double G4BinaryCascade::GetIonMass(G4int Z, G4int A)
2664	//----------------------------------------------------------------------------
2665	{
2666	G4double mass(0);
2667	if ( Z > 0 && A >= Z )
2668	{
2669	mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
2670
2671	} else if ( A > 0 && Z>0 )
2672	{
2673	// charge Z > A; will happen for light nuclei with pions involved.
2674	mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(A,A);
2675
2676	} else if ( A >= 0 && Z<=0 )
2677	{
2678	// all neutral, or empty nucleus
2679	mass = A * G4Neutron::Neutron()->GetPDGMass();
2680
2681	} else if ( A == 0 && std::abs(Z)<2 )
2682	{
2683	// empty nucleus, except maybe pions
2684	mass = 0;
2685
2686	} else
2687	{
2688	G4cerr << "G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
2689	<< A << "," << Z << ")" <<G4endl;
2690	G4Exception("G4BinaryCascade::GetIonMass() - giving up");
2691	}
2692	return mass;
2693	}
2694
2695	void G4BinaryCascade::PrintWelcomeMessage()
2696	{
2697	G4cout <<"Thank you for using G4BinaryCascade. "<<G4endl;
2698	}
2699
2700	//----------------------------------------------------------------------------
2701	void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision,
2702	G4KineticTrackVector * products)
2703	{
2704	G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2705
2706	G4KineticTrackVector debug1;
2707	debug1.push_back(collision->GetPrimary());
2708	PrintKTVector(&debug1,std::string(" Primary particle"));
2709	PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
2710	PrintKTVector(products,std::string(" Scatterer products"));
2711
2712	G4double thisExcitation(0);
2713	// excitation energy from this collision
2714	// initial state:
2715	G4double initial(0);
2716	G4KineticTrack * kt=collision->GetPrimary();
2717	initial += kt->Get4Momentum().e();
2718
2719	initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2720	initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2721	G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
2722	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2723	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2724	<< " " << initial << G4endl;;
2725
2726	G4KineticTrackVector ktv=collision->GetTargetCollection();
2727	for ( unsigned int it=0; it < ktv.size(); it++)
2728	{
2729	kt=ktv[it];
2730	initial += kt->Get4Momentum().e();
2731	thisExcitation += kt->GetDefinition()->GetPDGMass()
2732	- kt->Get4Momentum().e()
2733	- RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2734	// initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2735	// initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2736	G4cout << "Targ. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
2737	<< " " << kt->Get4Momentum().e()
2738	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2739	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2740	<< " " << initial <<" Excit " << thisExcitation << G4endl;;
2741	}
2742
2743	G4double final(0);
2744	G4double mass_out(0);
2745	G4int product_barions(0);
2746	if ( products )
2747	{
2748	for ( unsigned int it=0; it < products->size(); it++)
2749	{
2750	kt=(*products)[it];
2751	final += kt->Get4Momentum().e();
2752	final += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
2753	final += RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
2754	if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
2755	mass_out += kt->GetDefinition()->GetPDGMass();
2756	G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
2757	<< " " << kt->Get4Momentum().e()
2758	<< " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
2759	<< " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
2760	<< " " << final << G4endl;;
2761	}
2762	}
2763
2764
2765	G4int finalA = currentA;
2766	G4int finalZ = currentZ;
2767	if ( products )
2768	{
2769	finalA -= product_barions;
2770	finalZ -= GetTotalCharge(*products);
2771	}
2772	G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
2773	- (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA)
2774	+ mass_out);
2775	G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ
2776	<< " delta-mass " << delta<<G4endl;
2777	final+=delta;
2778	mass_out = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA);
2779	G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
2780	<< " " << final << " "
2781	<< mass_out<<" "
2782	<< currentInitialEnergy - final - mass_out
2783	<< G4endl;
2784	currentInitialEnergy-=final;
2785
2786	}
2787
2788	//----------------------------------------------------------------------------
2789	void G4BinaryCascade::DebugEpConservation(const G4HadProjectile & aTrack,
2790	G4ReactionProductVector* products)
2791	{
2792	G4ReactionProductVector::iterator iter;
2793	G4double Efinal(0);
2794	G4ThreeVector pFinal(0);
2795	if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
2796	{
2797	G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
2798	}
2799
2800	for(iter = products->begin(); iter != products->end(); ++iter)
2801	{
2802
2803	// G4cout << " Secondary E - Ekin / p " <<
2804	// (*iter)->GetDefinition()->GetParticleName() << " " <<
2805	// (*iter)->GetTotalEnergy() << " - " <<
2806	// (*iter)->GetKineticEnergy()<< " / " <<
2807	// (*iter)->GetMomentum().x() << " " <<
2808	// (*iter)->GetMomentum().y() << " " <<
2809	// (*iter)->GetMomentum().z() << G4endl;
2810	Efinal += (*iter)->GetTotalEnergy();
2811	pFinal += (*iter)->GetMomentum();
2812	}
2813
2814	// G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
2815	G4cout << "BIC E/p delta " <<
2816	(aTrack.Get4Momentum().e()+ the3DNucleus->GetMass() - Efinal)/MeV <<
2817	" MeV / mom " << (aTrack.Get4Momentum() - pFinal ) /MeV << G4endl;
2818
2819	}
2820
2821	//----------------------------------------------------------------------------

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

Download in other formats: