Context Navigation

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

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

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

Download in other formats: