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

Last change on this file since 960 was 819, checked in by garnier, 17 years ago

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