Context Navigation

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

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

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

Download in other formats: