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

Last change on this file since 1246 was 1228, checked in by garnier, 16 years ago

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