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

Last change on this file since 1353 was 1347, checked in by garnier, 15 years ago

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