Context Navigation

G4BinaryCascade.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago
geant4 tag 9.4
File size: 102.0 KB

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

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

Context Navigation

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

Download in other formats: