| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 #include "globals.hh"
28 #include "G4PhysicalConstants.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4Proton.hh"
31 #include "G4Neutron.hh"
32 #include "G4LorentzRotation.hh"
33 #include "G4BinaryCascade.hh"
34 #include "G4KineticTrackVector.hh"
35 #include "G4DecayKineticTracks.hh"
36 #include "G4ReactionProductVector.hh"
37 #include "G4Track.hh"
38 #include "G4V3DNucleus.hh"
39 #include "G4Fancy3DNucleus.hh"
40 #include "G4Scatterer.hh"
41 #include "G4MesonAbsorption.hh"
42 #include "G4ping.hh"
43 #include "G4Delete.hh"
44
45 #include "G4CollisionManager.hh"
46 #include "G4Absorber.hh"
47
48 #include "G4CollisionInitialState.hh"
49 #include "G4ListOfCollisions.hh"
50 #include "G4Fragment.hh"
51 #include "G4RKPropagation.hh"
52
53 #include "G4NuclearShellModelDensity.hh"
54 #include "G4NuclearFermiDensity.hh"
55 #include "G4FermiMomentum.hh"
56
57 #include "G4PreCompoundModel.hh"
58 #include "G4ExcitationHandler.hh"
59 #include "G4HadronicInteractionRegistry.hh"
60
61 #include "G4FermiPhaseSpaceDecay.hh"
62
63 #include "G4PreCompoundModel.hh"
64 #include "G4HadronicParameters.hh"
65
66 #include <algorithm>
67 #include "G4ShortLivedConstructor.hh"
68 #include <typeinfo>
69
70 #include "G4PhysicsModelCatalog.hh"
71
72 // turn on general debugging info, and consi
73
74 //#define debug_G4BinaryCascade 1
75
76 // more detailed debugging -- deprecated
77 //#define debug_H1_BinaryCascade 1
78
79 // specific debugging info per method or func
80 //#define debug_BIC_ApplyCollision 1
81 //#define debug_BIC_CheckPauli 1
82 //#define debug_BIC_CorrectFinalPandE 1
83 //#define debug_BIC_Propagate 1
84 //#define debug_BIC_Propagate_Excitation 1
85 //#define debug_BIC_Propagate_Collisions 1
86 //#define debug_BIC_Propagate_finals 1
87 //#define debug_BIC_DoTimeStep 1
88 //#define debug_BIC_CorrectBarionsOnBoundary 1
89 //#define debug_BIC_GetExcitationEnergy 1
90 //#define debug_BIC_DeexcitationProducts 1
91 //#define debug_BIC_FinalNucleusMomentum 1
92 //#define debug_BIC_Final4Momentum 1
93 //#define debug_BIC_FillVoidnucleus 1
94 //#define debug_BIC_FindFragments 1
95 //#define debug_BIC_BuildTargetList 1
96 //#define debug_BIC_FindCollision 1
97 //#define debug_BIC_return 1
98
99 //-------
100 //#if defined(debug_G4BinaryCascade)
101 #if 0
102 #define _CheckChargeAndBaryonNumber_(val) Ch
103 //#define debugCheckChargeAndBaryonNumberver
104 #else
105 #define _CheckChargeAndBaryonNumber_(val)
106 #endif
107 //#if defined(debug_G4BinaryCascade)
108 #if 0
109 #define _DebugEpConservation(val) DebugEpCo
110 //#define debugCheckChargeAndBaryonNumberver
111 #else
112 #define _DebugEpConservation(val)
113 #endif
114
115 G4int G4BinaryCascade::theBIC_ID = -1;
116
117 //
118 // C O N S T R U C T O R S A N D D E S T
119 //
120 G4BinaryCascade::G4BinaryCascade(G4VPreCompoun
121 G4VIntraNuclearTransportModel("Binary Cascade"
122 {
123 // initialise the resonance sector
124 G4ShortLivedConstructor ShortLived;
125 ShortLived.ConstructParticle();
126
127 theCollisionMgr = new G4CollisionManager;
128 theDecay=new G4BCDecay;
129 theImR.push_back(theDecay);
130 theLateParticle= new G4BCLateParticle;
131 G4MesonAbsorption * aAb=new G4MesonAbsorpt
132 theImR.push_back(aAb);
133 G4Scatterer * aSc=new G4Scatterer;
134 theH1Scatterer = new G4Scatterer;
135 theImR.push_back(aSc);
136
137 thePropagator = new G4RKPropagation;
138 theCurrentTime = 0.;
139 theBCminP = 45*MeV;
140 theCutOnP = 90*MeV;
141 theCutOnPAbsorb= 0*MeV; // No Absorption
142
143 // reuse existing pre-compound model
144 if(!ptr) {
145 G4HadronicInteraction* p =
146 G4HadronicInteractionRegistry::Instance()->F
147 G4VPreCompoundModel* pre = static_cast<G
148 if(!pre) { pre = new G4PreCompoundModel(
149 SetDeExcitation(pre);
150 }
151 theExcitationHandler = GetDeExcitation()->
152 SetMinEnergy(0.0*GeV);
153 SetMaxEnergy(10.1*GeV);
154 //PrintWelcomeMessage();
155 thePrimaryEscape = true;
156 thePrimaryType = 0;
157
158 SetEnergyMomentumCheckLevels(1.0*perCent,
159
160 // init data members
161 currentA=currentZ=0;
162 lateA=lateZ=0;
163 initialA=initialZ=0;
164 projectileA=projectileZ=0;
165 currentInitialEnergy=initial_nuclear_mass=
166 massInNucleus=0.;
167 theOuterRadius=0.;
168 theBIC_ID = G4PhysicsModelCatalog::GetMode
169 fBCDEBUG = G4HadronicParameters::Instance(
170 }
171
172 G4BinaryCascade::~G4BinaryCascade()
173 {
174 ClearAndDestroy(&theTargetList);
175 ClearAndDestroy(&theSecondaryList);
176 ClearAndDestroy(&theCapturedList);
177 delete thePropagator;
178 delete theCollisionMgr;
179 for(auto & ptr : theImR) { delete ptr; }
180 theImR.clear();
181 delete theLateParticle;
182 delete theH1Scatterer;
183 }
184
185 void G4BinaryCascade::ModelDescription(std::os
186 {
187 outFile << "G4BinaryCascade is an intra-nu
188 << "an incident hadron collides wi
189 << "final-state particles, one or
190 << "The resonances then decay hadr
191 << "are then propagated through th
192 << "trajectories until they re-int
193 << "This model is valid for incide
194 << "nucleons up to 10 GeV.\n"
195 << "The remaining excited nucleus
196 if (theDeExcitation)
197 {
198 outFile << theDeExcitation->GetM
199 theDeExcitation->DeExciteModelDe
200 }
201 else if (theExcitationHandler)
202 {
203 outFile << "G4ExcitationHandler
204 theExcitationHandler->ModelDesc
205 }
206 else
207 {
208 outFile << "void.\n";
209 }
210 outFile<< " \n";
211 }
212 void G4BinaryCascade::PropagateModelDescriptio
213 {
214 outFile << "G4BinaryCascade propagtes seco
215 << "energy model through the wound
216 << "Secondaries are followed after
217 << "within the nucleus are propaga
218 << "potential along curved traject
219 << "with a nucleon, decay, or leav
220 << "An interaction of a secondary
221 << "final-state particles, one or
222 << "Resonances decay hadronically
223 << "are in turn propagated through
224 << "trajectories until they re-int
225 << "This model is valid for pions
226 << "nucleons up to about 3.5 GeV.\
227 << "The remaining excited nucleus
228 if (theDeExcitation) // pre
229 {
230 outFile << theDeExcitation->GetModelName
231 theDeExcitation->DeExciteModelDescriptio
232 }
233 else if (theExcitationHandler) // de-ex
234 {
235 outFile << "G4ExcitationHandler"; //
236 theExcitationHandler->ModelDescription(
237 }
238 else
239 {
240 outFile << "void.\n";
241 }
242 outFile<< " \n";
243 }
244
245 //--------------------------------------------
246
247 //
248 // I M P L E M E N T A T I O N
249 //
250
251
252 //--------------------------------------------
253 G4HadFinalState * G4BinaryCascade::ApplyYourse
254 G4Nucleus & aNucleus)
255 //--------------------------------------------
256 {
257 if(fBCDEBUG) G4cerr << " ######### Binary
258
259 G4LorentzVector initial4Momentum = aTrack.
260 const G4ParticleDefinition * definition =
261
262 if(initial4Momentum.e()-initial4Momentum.m
263 ( definition==G4Neutron::NeutronDe
264 {
265 return theDeExcitation->ApplyYourself(
266 }
267
268 theParticleChange.Clear();
269 // initialize the G4V3DNucleus from G4Nucl
270 the3DNucleus = new G4Fancy3DNucleus;
271
272 // Build a KineticTrackVector with the G4T
273 G4KineticTrackVector * secondaries;// = ne
274 G4ThreeVector initialPosition(0., 0., 0.);
275
276 if(!fBCDEBUG)
277 {
278 if(definition!=G4Neutron::NeutronDefin
279 definition!=G4Proton::ProtonDe
280 definition!=G4PionPlus::PionPl
281 definition!=G4PionMinus::PionM
282 {
283 G4cerr << "You are trying to use G
284 G4cerr << "G4BinaryCascade should
285 G4cerr << "If you want to continue
286 G4cerr << "setenv I_Am_G4BinaryCas
287 throw G4HadronicException(__FILE__
288 }
289 }
290
291 // keep primary
292 thePrimaryType = definition;
293 thePrimaryEscape = false;
294
295 G4double timePrimary=aTrack.GetGlobalTime(
296
297 // try until an interaction will happen
298 G4ReactionProductVector * products = nullp
299 G4int interactionCounter = 0,collisionLoop
300 do
301 {
302 // reset status that could be changed
303 theCollisionMgr->ClearAndDestroy();
304
305 if(products != nullptr)
306 { // free memory from previous loop e
307 ClearAndDestroy(products);
308 delete products;
309 }
310
311 G4int massNumber=aNucleus.GetA_asInt()
312 the3DNucleus->Init(massNumber, aNucleu
313 thePropagator->Init(the3DNucleus);
314 G4KineticTrack * kt;
315 collisionLoopMaxCount = 200;
316 do // sample impact p
317 {
318 theCurrentTime=0;
319 G4double radius = the3DNucleus->Ge
320 initialPosition=GetSpherePoint(1.1
321 kt = new G4KineticTrack(definition
322 kt->SetState(G4KineticTrack::outsi
323 // secondaries has been cleared by
324 secondaries= new G4KineticTrackVec
325 secondaries->push_back(kt);
326 if(massNumber > 1) // 1H1 is speci
327 {
328 products = Propagate(secondari
329 } else {
330 products = Propagate1H1(second
331 }
332 // until we FIND a collision ... o
333 } while(! products && --collisionLoopM
334
335 if(++interactionCounter>99) break;
336 // ...until we find an ALLOWED collisi
337 } while(products && products->size() == 0)
338
339 if(products && products->size()>0)
340 {
341 // G4cout << "BIC Applyyourself: numb
342
343 // Fill the G4ParticleChange * with pr
344 theParticleChange.SetStatusChange(stop
345 G4ReactionProductVector::iterator iter
346
347 for(iter = products->begin(); iter !=
348 {
349 G4DynamicParticle * aNewDP =
350 new G4DynamicParticle((*it
351 (*iter)->GetTotalE
352 (*iter)->GetMoment
353 G4HadSecondary aNew = G4HadSecondary
354 G4double time=(*iter)->GetFormatio
355 if(time < 0.0) { time = 0.0; }
356 aNew.SetTime(timePrimary + time);
357 aNew.SetCreatorModelID((*iter)->Ge
358 aNew.SetParentResonanceDef((*iter)
359 aNew.SetParentResonanceID((*iter)-
360 theParticleChange.AddSecondary(aNe
361 }
362
363 //DebugFinalEpConservation(aTrack, pr
364
365
366 } else { // no interaction, return primar
367 if(fBCDEBUG) G4cerr << " ######### Bin
368 theParticleChange.SetStatusChange(isAl
369 theParticleChange.SetEnergyChange(aTra
370 theParticleChange.SetMomentumChange(aT
371 }
372
373 if ( products )
374 {
375 ClearAndDestroy(products);
376 delete products;
377 }
378
379 delete the3DNucleus;
380 the3DNucleus = nullptr;
381
382 if(fBCDEBUG) G4cerr << " ######### Binary
383
384 return &theParticleChange;
385 }
386 //--------------------------------------------
387 G4ReactionProductVector * G4BinaryCascade::Pro
388 G4KineticTrackVector * secondaries, G4
389 //--------------------------------------------
390 {
391 G4ping debug("debug_G4BinaryCascade");
392 #ifdef debug_BIC_Propagate
393 G4cout << "G4BinaryCascade Propagate start
394 #endif
395
396 the3DNucleus=aNucleus;
397 G4ReactionProductVector * products = new G
398 theOuterRadius = the3DNucleus->GetOuterRad
399 theCurrentTime=0;
400 theProjectile4Momentum=G4LorentzVector(0,0
401 theMomentumTransfer=G4ThreeVector(0,0,0);
402 // build theSecondaryList, theProjectileLi
403 ClearAndDestroy(&theCapturedList);
404 ClearAndDestroy(&theSecondaryList);
405 theSecondaryList.clear();
406 ClearAndDestroy(&theFinalState);
407 std::vector<G4KineticTrack *>::iterator it
408 theCollisionMgr->ClearAndDestroy();
409
410 theCutOnP=90*MeV;
411 if(the3DNucleus->GetMass()>30) theCutOnP =
412 if(the3DNucleus->GetMass()>60) theCutOnP =
413 if(the3DNucleus->GetMass()>120) theCutOnP
414
415
416 BuildTargetList();
417
418 #ifdef debug_BIC_GetExcitationEnergy
419 G4cout << "ExcitationEnergy0 " << GetExcit
420 #endif
421
422 thePropagator->Init(the3DNucleus);
423
424 G4bool success = BuildLateParticleCollisio
425 if (! success ) // fails if no excitatio
426 {
427 products=HighEnergyModelFSProducts(prod
428 ClearAndDestroy(secondaries);
429 delete secondaries;
430
431 #ifdef debug_G4BinaryCascade
432 G4cout << "G4BinaryCascade::Propagate:
433 #endif
434
435 return products;
436 }
437 // check baryon and charge ...
438
439 _CheckChargeAndBaryonNumber_("lateparticle
440 _DebugEpConservation(" be4 findcollisions"
441
442 // if called stand alone find first collis
443 FindCollisions(&theSecondaryList);
444
445
446 if(theCollisionMgr->Entries() == 0 )
447 {
448 //G4cout << " no collsions -> return 0
449 delete products;
450 #ifdef debug_BIC_return
451 G4cout << "return @ begin2, no collis
452 #endif
453 return nullptr;
454 }
455
456 // end of initialization: do the job now
457 // loop until there are no more collisions
458
459
460 G4bool haveProducts = false;
461 #ifdef debug_BIC_Propagate_Collisions
462 G4int collisionCount=0;
463 #endif
464 G4int collisionLoopMaxCount=1000000;
465 while(theCollisionMgr->Entries() > 0 && cu
466 {
467 if(Absorb()) { // absorb secondaries, p
468 haveProducts = true;
469 }
470 if(Capture()) { // capture secondaries,
471 haveProducts = true;
472 }
473
474 // propagate to the next collision if
475 // by previous absorption or capture)
476 if(theCollisionMgr->Entries() > 0)
477 {
478 G4CollisionInitialState *
479 nextCollision = theCollisionMgr->G
480 #ifdef debug_BIC_Propagate_Collisions
481 G4cout << " NextCollision * , Tim
482 <<nextCollision->GetCollis
483 theCurrentTime<< G4endl;
484 #endif
485 if (!DoTimeStep(nextCollision->Get
486 {
487 // Check if nextCollision is s
488 if (theCollisionMgr->GetNextCo
489 {
490 nextCollision = nullptr;
491 }
492 }
493 //_DebugEpConservation("Stepped");
494
495 if( nextCollision )
496 {
497 if (ApplyCollision(nextCollisi
498 {
499 //G4cerr << "ApplyCollisio
500 haveProducts = true;
501 #ifdef debug_BIC_Propagate_Collisions
502 collisionCount++;
503 #endif
504
505 } else {
506 //G4cerr << "ApplyCollisio
507 theCollisionMgr->RemoveCol
508 }
509 }
510 }
511 }
512
513 //--------- end of on Collisions
514 //G4cout << "currentZ @ end loop " << curr
515 G4int nProtons(0);
516 for(iter = theTargetList.begin(); iter !=
517 {
518 if ( (*iter)->GetDefinition() == G4Proto
519 }
520 if ( ! theTargetList.size() || ! nProtons
521 // nucleus completely destroyed, fill
522 products = FillVoidNucleusProducts(prod
523 #ifdef debug_BIC_return
524 G4cout << "return @ Z=0 after collisio
525 PrintKTVector(&theSecondaryList,std::s
526 G4cout << "theTargetList size: " << th
527 PrintKTVector(&theTargetList,std::stri
528 PrintKTVector(&theCapturedList,std::st
529
530 G4cout << " ExcitE be4 Correct : " <<G
531 G4cout << " Mom Transfered to nucleus
532 PrintKTVector(&theFinalState,std::stri
533 G4cout << "returned products: " << pro
534 _CheckChargeAndBaryonNumber_("destroye
535 _DebugEpConservation("destroyed Nucleu
536 #endif
537
538 return products;
539 }
540
541 // No more collisions: absorb, capture and
542 if(Absorb()) {
543 haveProducts = true;
544 // G4cout << "Absorb sucess " << G4end
545 }
546
547 if(Capture()) {
548 haveProducts = true;
549 // G4cout << "Capture sucess " << G4en
550 }
551
552 if(!haveProducts) // no collisions happen
553 {
554 #ifdef debug_BIC_return
555 G4cout << "return 3, no products "<< G
556 #endif
557 return products;
558 }
559
560
561 #ifdef debug_BIC_Propagate
562 G4cout << " Momentum transfer to Nucleus "
563 G4cout << " Stepping particles out......
564 #endif
565
566 StepParticlesOut();
567 _DebugEpConservation("stepped out");
568
569
570 if ( theSecondaryList.size() > 0 )
571 {
572 #ifdef debug_G4BinaryCascade
573 G4cerr << "G4BinaryCascade: Warning, h
574 PrintKTVector(&theSecondaryList, "acti
575 #endif
576 // add left secondaries to FinalSate
577 for ( iter =theSecondaryList.begin();
578 {
579 theFinalState.push_back(*iter);
580 }
581 theSecondaryList.clear();
582
583 }
584 while ( theCollisionMgr->Entries() > 0 )
585 {
586 #ifdef debug_G4BinaryCascade
587 G4cerr << " Warning: remove left over
588 #endif
589 theCollisionMgr->RemoveCollision(theCo
590 }
591
592 #ifdef debug_BIC_Propagate_Excitation
593
594 PrintKTVector(&theSecondaryList,std::strin
595 G4cout << "theTargetList size: " << theTar
596 // PrintKTVector(&theTargetList,std::stri
597 PrintKTVector(&theCapturedList,std::string
598
599 G4cout << " ExcitE be4 Correct : " <<GetEx
600 G4cout << " Mom Transfered to nucleus : "
601 PrintKTVector(&theFinalState,std::string("
602 #endif
603
604 //
605
606
607 G4double ExcitationEnergy=GetExcitationEne
608
609 #ifdef debug_BIC_Propagate_finals
610 PrintKTVector(&theFinalState,std::string("
611 G4cout << " Excitation Energy prefinal, #
612 << ExcitationEnergy << " "
613 << collisionCount << " "
614 << theFinalState.size() << " "
615 << theCapturedList.size()<<G4endl;
616 #endif
617
618 if (ExcitationEnergy < 0 )
619 {
620 G4int maxtry=5, ntry=0;
621 do {
622 CorrectFinalPandE();
623 ExcitationEnergy=GetExcitationEner
624 } while ( ++ntry < maxtry && Excitatio
625 }
626 _DebugEpConservation("corrected");
627
628 #ifdef debug_BIC_Propagate_finals
629 PrintKTVector(&theFinalState,std::string("
630 G4cout << " Excitation Energy final, #col
631 << ExcitationEnergy << " "
632 << collisionCount << " "
633 << theFinalState.size() << " "
634 << theCapturedList.size()<<G4endl;
635 #endif
636
637
638 if ( ExcitationEnergy < 0. )
639 {
640 #ifdef debug_G4BinaryCascade
641 G4cerr << "G4BinaryCascade-W
642 G4cerr <<ExcitationEnergy<<G
643 PrintKTVector(&theFinalState,
644 PrintKTVector(&theCapturedLis
645 G4cout << "negative ExE:Final
646 << " "<< GetFinal4Mome
647 << "negative ExE:Final
648 << " "<< GetFinalNucleusMome
649 #endif
650 #ifdef debug_BIC_return
651 G4cout << " negative Exci
652 G4cout << "return 4, excit
653 #endif
654
655 ClearAndDestroy(products);
656 return products; // return empty pro
657 }
658
659 G4ReactionProductVector * precompoundProdu
660
661
662 G4DecayKineticTracks decay(&theFinalState)
663 _DebugEpConservation("decayed");
664
665 products= ProductsAddFinalState(products,
666
667 products= ProductsAddPrecompound(products,
668
669 // products=ProductsAddFakeGamma(products);
670
671
672 thePrimaryEscape = true;
673
674 #ifdef debug_BIC_return
675 G4cout << "BIC: return @end, all ok "<< G4
676 //G4cout << " return products " << produc
677 #endif
678
679 return products;
680 }
681
682 //--------------------------------------------
683 G4double G4BinaryCascade::GetExcitationEnergy(
684 //--------------------------------------------
685 {
686
687 // get A and Z for the residual nucleus
688 #if defined(debug_G4BinaryCascade) || defined(
689 G4int finalA = theTargetList.size()+theCap
690 G4int finalZ = GetTotalCharge(theTargetLis
691 if ( (currentA - finalA) != 0 || (currentZ
692 {
693 G4cerr << "G4BIC:GetExcitationEnergy()
694 << "("<< currentA << "," << fi
695 }
696
697 #endif
698
699 G4double excitationE(0);
700 G4double nucleusMass(0);
701 if(currentZ>.5)
702 {
703 nucleusMass = GetIonMass(currentZ,curr
704 }
705 else if (currentZ==0 )
706 {
707 if(currentA == 1) {nucleusMass = G4Neu
708 else {nucleusMass = GetFi
709 }
710 else
711 {
712 #ifdef debug_G4BinaryCascade
713 G4cout << "G4BinaryCascade::GetExcitat
714 << currentA << "," << currentZ
715 #endif
716 return 0;
717 }
718
719 #ifdef debug_BIC_GetExcitationEnergy
720 G4ping debug("debug_ExcitationEnergy");
721 debug.push_back("====> current A, Z");
722 debug.push_back(currentZ);
723 debug.push_back(currentA);
724 debug.push_back("====> final A, Z");
725 debug.push_back(finalZ);
726 debug.push_back(finalA);
727 debug.push_back(nucleusMass);
728 debug.push_back(GetFinalNucleusMomentum().
729 debug.dump();
730 // PrintKTVector(&theTargetList, std::str
731 //PrintKTVector(&theCapturedList, std::str
732 #endif
733
734 excitationE = GetFinalNucleusMomentum().ma
735
736 //G4double exE2 = GetFinal4Momentum().mag(
737
738 //G4cout << "old/new excitE " << excitatio
739
740 #ifdef debug_BIC_GetExcitationEnergy
741 // ------ debug
742 if ( excitationE < 0 )
743 {
744 G4cout << "negative ExE final Ion mass
745 G4LorentzVector Nucl_mom=GetFinalNucle
746 if(finalZ>.5) G4cout << " Final nuclmo
747 << " (A,Z)=("<< finalA
748 << " mass " << nucleusM
749 << " excitE " << excita
750
751
752 G4int A = the3DNucleus->GetMassNumber(
753 G4int Z = the3DNucleus->GetCharge();
754 G4double initialExc(0);
755 if(Z>.5)
756 {
757 initialExc = theInitial4Mom.mag()-
758 G4cout << "GetExcitationEnergy: In
759 }
760 }
761
762 #endif
763
764 return excitationE;
765 }
766
767
768 //--------------------------------------------
769 //
770 // P R I V A T E M E M B E R F U N C
771 //
772 //--------------------------------------------
773
774 //--------------------------------------------
775 void G4BinaryCascade::BuildTargetList()
776 //--------------------------------------------
777 {
778
779 if(!the3DNucleus->StartLoop())
780 {
781 // G4cerr << "G4BinaryCascade::Buil
782 // << G4endl;
783 return;
784 }
785
786 ClearAndDestroy(&theTargetList); // clear
787
788 G4Nucleon * nucleon;
789 const G4ParticleDefinition * definition;
790 G4ThreeVector pos;
791 G4LorentzVector mom;
792 // if there are nucleon hit by higher ener
793 initialZ=the3DNucleus->GetCharge();
794 initialA=the3DNucleus->GetMassNumber();
795 initial_nuclear_mass=GetIonMass(initialZ,i
796 theInitial4Mom = G4LorentzVector(0,0,0,ini
797 currentA=0;
798 currentZ=0;
799 while((nucleon = the3DNucleus->GetNextNucl
800 {
801 // check if nucleon is hit by higher e
802 if ( ! nucleon->AreYouHit() )
803 {
804 definition = nucleon->GetDefinitio
805 pos = nucleon->GetPosition();
806 mom = nucleon->GetMomentum();
807 // G4cout << "Nucleus " << pos.
808 //theInitial4Mom += mom;
809 // the potential inside the
810 mom.setE( std::sqrt( mom.vect().ma
811 G4KineticTrack * kt = new G4Kineti
812 kt->SetState(G4KineticTrack::insid
813 kt->SetNucleon(nucleon);
814 theTargetList.push_back(kt);
815 ++currentA;
816 if (definition->GetPDGCharge() > .
817 }
818
819 #ifdef debug_BIC_BuildTargetList
820 else { G4cout << "nucleon is hit" << n
821 #endif
822 }
823 massInNucleus = 0;
824 if(currentZ>.5)
825 {
826 massInNucleus = GetIonMass(currentZ,cu
827 } else if (currentZ==0 && currentA>=1 )
828 {
829 massInNucleus = currentA * G4Neutron::
830 } else
831 {
832 G4cerr << "G4BinaryCascade::BuildTarge
833 << currentA << "," << currentZ
834 throw G4HadronicException(__FILE__, __
835 }
836 currentInitialEnergy= theInitial4Mom.e() +
837
838 #ifdef debug_BIC_BuildTargetList
839 G4cout << "G4BinaryCascade::BuildTargetLis
840 << currentA << "," << currentZ <<
841 ", theInitial4Mom " << theInitial4
842 ", currentInitialEnergy " << curre
843 #endif
844
845 }
846
847 //--------------------------------------------
848 G4bool G4BinaryCascade::BuildLateParticleColl
849 //--------------------------------------------
850 {
851 G4bool success(false);
852 std::vector<G4KineticTrack *>::iterator ite
853
854 lateA=lateZ=0;
855 projectileA=projectileZ=0;
856
857 G4double StartingTime=DBL_MAX; // Se
858 for(iter = secondaries->begin(); iter != se
859 {
860 if((*iter)->GetFormationTime() < Startin
861 StartingTime = (*iter)->GetFormationT
862 }
863
864 //PrintKTVector(secondaries, "initial late
865 G4LorentzVector lateParticles4Momentum(0,0,
866 for(iter = secondaries->begin(); iter != se
867 {
868 // G4cout << " Formation time : " << (*
869 // << (*iter)->GetFormationTime() << G
870 G4double FormTime = (*iter)->GetFormatio
871 (*iter)->SetFormationTime(FormTime);
872 if( (*iter)->GetState() == G4KineticTrac
873 {
874 FindLateParticleCollision(*iter);
875 lateParticles4Momentum += (*iter)->Ge
876 lateA += (*iter)->GetDefinition()->Ge
877 lateZ += G4lrint((*iter)->GetDefiniti
878 //PrintKTVector(*iter, "late particle
879 } else
880 {
881 theSecondaryList.push_back(*iter);
882 //PrintKTVector(*iter, "incoming part
883 theProjectile4Momentum += (*iter)->Ge
884 projectileA += (*iter)->GetDefinition
885 projectileZ += G4lrint((*iter)->GetDe
886 #ifdef debug_BIC_Propagate
887 G4cout << " Adding initial secondary
888 << " time" << (*iter)->GetForma
889 << ", state " << (*iter)->GetSt
890 #endif
891 }
892 }
893 //theCollisionMgr->Print();
894 const G4HadProjectile * primary = GetPrimar
895
896 if (primary){
897 G4LorentzVector mom=primary->Get4Momentu
898 theProjectile4Momentum += mom;
899 projectileA = primary->GetDefinition()->
900 projectileZ = G4lrint(primary->GetDefini
901 // now check if "excitation" energy left
902 G4double excitation= theProjectile4Momen
903 #ifdef debug_BIC_GetExcitationEnergy
904 G4cout << "BIC: Proj.e, nucl initial, nu
905 << theProjectile4Momentum << ", "
906 << initial_nuclear_mass<< ", " <<
907 << lateParticles4Momentum << G4end
908 G4cout << "BIC: Proj.e / initial excitat
909 #endif
910 success = excitation > 0;
911 #ifdef debug_G4BinaryCascade
912 if ( ! success ) {
913 G4cout << "G4BinaryCascade::BuildLate
914 //PrintKTVector(secondaries);
915 }
916 #endif
917 } else {
918 // no primary from HE model -> cascade
919 success=true;
920 }
921
922 if (success) {
923 secondaries->clear(); // Don't leave "G4
924 delete secondaries;
925 }
926 return success;
927 }
928
929 //--------------------------------------------
930 G4ReactionProductVector * G4BinaryCascade::De
931 //--------------------------------------------
932 {
933 // find a fragment and call the precompound
934 G4Fragment * fragment = nullptr;
935 G4ReactionProductVector * precompoundProduc
936
937 G4LorentzVector pFragment(0);
938 // G4cout << " final4mon " << GetFinal4Mome
939
940 fragment = FindFragments();
941
942 if(fragment)
943 {
944 if(fragment->GetA_asInt() >1)
945 {
946 pFragment=fragment->GetMomentum();
947 // G4cout << " going to preco with fr
948 if (theDeExcitation) /
949 {
950 precompoundProducts= theDeExcitati
951 }
952 else if (theExcitationHandler) //
953 {
954 precompoundProducts=theExcitationH
955 }
956
957 } else
958 { // f
959 if (theTargetList.size() + theCapture
960 throw G4HadronicException(__FILE__
961 }
962
963 std::vector<G4KineticTrack *>::iterat
964 if ( theTargetList.size() == 1 ) {i=
965 if ( theCapturedList.size() == 1 ) {i
966 G4ReactionProduct * aNew = new G4Reac
967 aNew->SetTotalEnergy((*i)->GetDefinit
968 aNew->SetCreatorModelID(theBIC_ID);
969 aNew->SetParentResonanceDef((*i)->GetParent
970 aNew->SetParentResonanceID((*i)->GetParentR
971 aNew->SetMomentum(G4ThreeVector(0));/
972 precompoundProducts = new G4ReactionP
973 precompoundProducts->push_back(aNew);
974 } // End of f
975 delete fragment;
976 fragment=nullptr;
977
978 } else // End of
979 { // No fra
980
981 precompoundProducts = DecayVoidNucleus()
982 }
983 #ifdef debug_BIC_DeexcitationProducts
984
985 G4LorentzVector fragment_momentum=GetFi
986 G4LorentzVector Preco_momentum;
987 if ( precompoundProducts )
988 {
989 std::vector<G4ReactionProduct *>::it
990 for(j = precompoundProducts->begin()
991 {
992 G4LorentzVector pProduct((*j)->Ge
993 Preco_momentum += pProduct;
994 }
995 }
996 G4cout << "finalNuclMom / sum preco pro
997 << " delta E "<< fragment_momentum.
998
999 #endif
1000
1001 return precompoundProducts;
1002 }
1003
1004 //-------------------------------------------
1005 G4ReactionProductVector * G4BinaryCascade::D
1006 //-------------------------------------------
1007 {
1008 G4ReactionProductVector * result = nullptr
1009 if ( (theTargetList.size()+theCapturedList
1010 {
1011 result = new G4ReactionProductVector;
1012 std::vector<G4KineticTrack *>::iterator
1013 G4LorentzVector aVec;
1014 std::vector<G4double> masses;
1015 G4double sumMass(0);
1016
1017 if ( theTargetList.size() != 0)
1018 {
1019 for ( aNuc=theTargetList.begin(); aN
1020 {
1021 G4double mass=(*aNuc)->GetDefinit
1022 masses.push_back(mass);
1023 sumMass += mass;
1024 }
1025 }
1026
1027 if ( theCapturedList.size() != 0)
1028 {
1029 for(aNuc = theCapturedList.begin();
1030 {
1031 G4double mass=(*aNuc)->GetDefinit
1032 masses.push_back(mass);
1033 sumMass += mass;
1034 }
1035 }
1036
1037 G4LorentzVector finalP=GetFinal4Momentu
1038 G4FermiPhaseSpaceDecay decay;
1039 // G4cout << " some neutrons? " << mass
1040 // G4cout<< theTargetList.size()<<" "<<
1041
1042 G4double eCMS=finalP.mag();
1043 if ( eCMS < sumMass )
1044 {
1045 eCMS=sumMass + 2*MeV*masses.size();
1046 finalP.setE(std::sqrt(finalP.vect().
1047 }
1048
1049 precompoundLorentzboost.set(finalP.boos
1050 std::vector<G4LorentzVector*> * momenta
1051 std::vector<G4LorentzVector*>::iterator
1052
1053
1054 if ( theTargetList.size() != 0)
1055 {
1056 for ( aNuc=theTargetList.begin();
1057 (aNuc != theTargetList.end())
1058 aNuc++, aMom++ )
1059 {
1060 G4ReactionProduct * aNew = new G4
1061 aNew->SetTotalEnergy((*aMom)->e()
1062 aNew->SetMomentum((*aMom)->vect()
1063 aNew->SetCreatorModelID(theBIC_ID
1064 aNew->SetParentResonanceDef((*aNuc)->Ge
1065 aNew->SetParentResonanceID((*aNuc)->Get
1066 result->push_back(aNew);
1067 delete *aMom;
1068 }
1069 }
1070
1071 if ( theCapturedList.size() != 0)
1072 {
1073 for ( aNuc=theCapturedList.begin();
1074 (aNuc != theCapturedList.end()
1075 aNuc++, aMom++ )
1076 {
1077 G4ReactionProduct * aNew = new G4
1078 aNew->SetTotalEnergy((*aMom)->e()
1079 aNew->SetMomentum((*aMom)->vect()
1080 aNew->SetCreatorModelID(theBIC_ID
1081 aNew->SetParentResonanceDef((*aNu
1082 aNew->SetParentResonanceID((*aNuc
1083 result->push_back(aNew);
1084 delete *aMom;
1085 }
1086 }
1087
1088 delete momenta;
1089 }
1090 return result;
1091 } // End if(!fragment)
1092
1093 //-------------------------------------------
1094 G4ReactionProductVector * G4BinaryCascade::Pr
1095 //-------------------------------------------
1096 {
1097 // fill in products the outgoing particles
1098 std::size_t i(0);
1099 #ifdef debug_BIC_Propagate_finals
1100 G4LorentzVector mom_fs;
1101 #endif
1102 for(i = 0; i< fs.size(); i++)
1103 {
1104 G4KineticTrack * kt = fs[i];
1105 G4ReactionProduct * aNew = new G4Reac
1106 aNew->SetMomentum(kt->Get4Momentum().
1107 aNew->SetTotalEnergy(kt->Get4Momentum
1108 aNew->SetNewlyAdded(kt->IsParticipant
1109 aNew->SetCreatorModelID(theBIC_ID);
1110 aNew->SetParentResonanceDef(kt->GetPa
1111 aNew->SetParentResonanceID(kt->GetPar
1112 products->push_back(aNew);
1113
1114 #ifdef debug_BIC_Propagate_finals
1115 mom_fs += kt->Get4Momentum();
1116 G4cout <<kt->GetDefinition()->GetPart
1117 G4cout << " Particle Ekin " << aNew->
1118 G4cout << ", is " << (kt->GetDefiniti
1119 (kt->GetDefinition()->IsShort
1120 G4cout << G4endl;
1121 #endif
1122
1123 }
1124 #ifdef debug_BIC_Propagate_finals
1125 G4cout << " Final state momentum " << mom
1126 #endif
1127
1128 return products;
1129 }
1130 //-------------------------------------------
1131 G4ReactionProductVector * G4BinaryCascade::Pr
1132 //-------------------------------------------
1133 {
1134 G4LorentzVector pSumPreco(0), pPreco(0);
1135
1136 if ( precompoundProducts )
1137 {
1138 for(auto j = precompoundProducts->cbegi
1139 {
1140 // boost back to system of moving nu
1141 G4LorentzVector pProduct((*j)->GetMo
1142 pPreco+= pProduct;
1143 #ifdef debug_BIC_Propagate_finals
1144 G4cout << "BIC: pProduct be4 boost "
1145 #endif
1146 pProduct *= precompoundLorentzboost;
1147 #ifdef debug_BIC_Propagate_finals
1148 G4cout << "BIC: pProduct aft boost "
1149 #endif
1150 pSumPreco += pProduct;
1151 (*j)->SetTotalEnergy(pProduct.e());
1152 (*j)->SetMomentum(pProduct.vect());
1153 (*j)->SetNewlyAdded(true);
1154 products->push_back(*j);
1155 }
1156 // G4cout << " unboosted preco result m
1157 // G4cout << " preco result mom " << pS
1158 precompoundProducts->clear();
1159 delete precompoundProducts;
1160 }
1161 return products;
1162 }
1163 //-------------------------------------------
1164 void G4BinaryCascade::FindCollisions(G4Kinet
1165 //-------------------------------------------
1166 {
1167 for(auto i=secondaries->cbegin(); i!=seco
1168 {
1169 for(auto j=theImR.cbegin(); j!=theImR
1170 {
1171 // G4cout << "G4BinaryCascad
1172 const std::vector<G4CollisionInit
1173 = (*j)->GetCollisions(*i, theTarg
1174 for(std::size_t count=0; count<aC
1175 {
1176 theCollisionMgr->AddCollision
1177 //4cout << "=================
1178 //theCollisionMgr->Print();
1179 }
1180 }
1181 }
1182 }
1183
1184
1185 //-------------------------------------------
1186 void G4BinaryCascade::FindDecayCollision(G4K
1187 //-------------------------------------------
1188 {
1189 const auto& aCandList = theDecay->GetColl
1190 for(std::size_t count=0; count<aCandList.
1191 {
1192 theCollisionMgr->AddCollision(aCandLi
1193 }
1194 }
1195
1196 //-------------------------------------------
1197 void G4BinaryCascade::FindLateParticleCollis
1198 //-------------------------------------------
1199 {
1200
1201 G4double tin=0., tout=0.;
1202 if (((G4RKPropagation*)thePropagator)->Ge
1203 {
1204 if ( tin > 0 )
1205 {
1206 secondary->SetState(G4KineticTrac
1207 } else if ( tout > 0 )
1208 {
1209 secondary->SetState(G4KineticTrac
1210 } else {
1211 //G4cout << "G4BC set miss , tin,
1212 secondary->SetState(G4KineticTrac
1213 }
1214 } else {
1215 secondary->SetState(G4KineticTrack::m
1216 //G4cout << "G4BC set miss ,no inters
1217 }
1218
1219
1220 #ifdef debug_BIC_FindCollision
1221 G4cout << "FindLateP Particle, 4-mom, tim
1222 << secondary->GetDefinition()->Ge
1223 << secondary->Get4Momentum()
1224 << " times " << tin << " " << to
1225 << secondary->GetState() << G4end
1226 #endif
1227
1228 const auto& aCandList = theLateParticle->
1229 for(std::size_t count=0; count<aCandList.
1230 {
1231 #ifdef debug_BIC_FindCollision
1232 G4cout << " Adding a late Col : " <<
1233 #endif
1234 theCollisionMgr->AddCollision(aCandLi
1235 }
1236 }
1237
1238
1239 //-------------------------------------------
1240 G4bool G4BinaryCascade::ApplyCollision(G4Coll
1241 //-------------------------------------------
1242 {
1243 G4KineticTrack * primary = collision->Get
1244
1245 #ifdef debug_BIC_ApplyCollision
1246 G4cerr << "G4BinaryCascade::ApplyCollisio
1247 theCollisionMgr->Print();
1248 G4cout << "ApplyCollisions : projte 4mom
1249 #endif
1250
1251 G4KineticTrackVector target_collection=co
1252 G4bool haveTarget=target_collection.size(
1253 if( haveTarget && (primary->GetState() !=
1254 {
1255 #ifdef debug_G4BinaryCascade
1256 G4cout << "G4BinaryCasacde::ApplyColl
1257 PrintKTVector(primary,std::string("pr
1258 PrintKTVector(&target_collection,std:
1259 collision->Print();
1260 G4cout << G4endl;
1261 theCollisionMgr->Print();
1262 //*GF* throw G4HadronicException(
1263 #endif
1264 return false;
1265 // } else {
1266 // G4cout << "G4BinaryCasacde::ApplyCol
1267 // PrintKTVector(primary,std::string("p
1268 // G4double tin=0., tout=0.;
1269 // if (((G4RKPropagation*)thePropagator
1270 // {
1271 // G4cout << "tin tout: " << tin <<
1272 // }
1273
1274 }
1275
1276 G4LorentzVector mom4Primary=primary->Get4
1277
1278 G4int initialBaryon(0);
1279 G4int initialCharge(0);
1280 if (primary->GetState() == G4KineticTrack
1281 {
1282 initialBaryon = primary->GetDefinitio
1283 initialCharge = G4lrint(primary->GetD
1284 }
1285
1286 // for primary resonances, subtract neutr
1287 G4double initial_Efermi=CorrectShortlived
1288 //***************************************
1289
1290
1291 G4KineticTrackVector * products = collisi
1292
1293 #ifdef debug_BIC_ApplyCollision
1294 DebugApplyCollisionFail(collision, produc
1295 #endif
1296
1297 // reset primary to initial state, in cas
1298 primary->Set4Momentum(mom4Primary);
1299
1300 G4bool lateParticleCollision= (!haveTarge
1301 G4bool decayCollision= (!haveTarget) && p
1302 G4bool Success(true);
1303
1304
1305 #ifdef debug_G4BinaryCascade
1306 G4int lateBaryon(0), lateCharge(0);
1307 #endif
1308
1309 if ( lateParticleCollision )
1310 { // for late particles, reset charges
1311 //G4cout << "lateP, initial B C state
1312 // << initialCharge<< " " << p
1313 #ifdef debug_G4BinaryCascade
1314 lateBaryon = initialBaryon;
1315 lateCharge = initialCharge;
1316 #endif
1317 initialBaryon=initialCharge=0;
1318 lateA -= primary->GetDefinition()->Ge
1319 lateZ -= G4lrint(primary->GetDefiniti
1320 }
1321
1322 initialBaryon += collision->GetTargetBary
1323 initialCharge += G4lrint(collision->GetTa
1324 if (!lateParticleCollision)
1325 {
1326 if( !products || products->size()==0 |
1327 {
1328 #ifdef debug_BIC_ApplyCollision
1329 if (products) G4cout << " ======Fai
1330 G4cerr << "G4BinaryCascade::ApplyCo
1331 #endif
1332 Success=false;
1333 }
1334
1335
1336
1337 if (Success && primary->GetState() ==
1338 if (! CorrectShortlivedFinalsForFer
1339 Success=false;
1340 }
1341 }
1342 }
1343
1344 #ifdef debug_BIC_ApplyCollision
1345 DebugApplyCollision(collision, products);
1346 #endif
1347
1348 if ( ! Success ){
1349 if (products) ClearAndDestroy(product
1350 if ( decayCollision ) FindDecayCollis
1351 delete products;
1352 products=nullptr;
1353 return false;
1354 }
1355
1356 G4int finalBaryon(0);
1357 G4int finalCharge(0);
1358 G4KineticTrackVector toFinalState;
1359 for(auto i=products->cbegin(); i!=product
1360 {
1361 if ( ! lateParticleCollision )
1362 {
1363 (*i)->SetState(primary->GetState(
1364 if ( (*i)->GetState() == G4Kineti
1365 finalBaryon+=(*i)->GetDefinit
1366 finalCharge+=G4lrint((*i)->Ge
1367 } else {
1368 G4double tin=0., tout=0.;
1369 if (((G4RKPropagation*)theProp
1370 tin < 0 && tout > 0 )
1371 {
1372 PrintKTVector((*i),"particl
1373 G4cout << "tin tout: " <<
1374 }
1375 }
1376 } else {
1377 G4double tin=0., tout=0.;
1378 if (((G4RKPropagation*)thePropaga
1379 {
1380 //G4cout << "tin tout: " << t
1381 if ( tin > 0 )
1382 {
1383 (*i)->SetState(G4KineticT
1384 }
1385 else if ( tout > 0 )
1386 {
1387 (*i)->SetState(G4KineticT
1388 finalBaryon+=(*i)->GetDef
1389 finalCharge+=G4lrint((*i)
1390 }
1391 else
1392 {
1393 (*i)->SetState(G4KineticT
1394 toFinalState.push_back((*
1395 }
1396 } else
1397 {
1398 (*i)->SetState(G4KineticTrack
1399 //G4cout << " G4BC - miss -la
1400 toFinalState.push_back((*i));
1401 }
1402
1403 }
1404 }
1405 if(!toFinalState.empty())
1406 {
1407 theFinalState.insert(theFinalState.ce
1408 toFinalState.cbegin(),toFinal
1409 std::vector<G4KineticTrack *>::iterat
1410 for(auto iter1=toFinalState.cbegin();
1411 {
1412 iter2 = std::find(products->begin
1413 if ( iter2 != products->cend() )
1414 }
1415 theCollisionMgr->RemoveTracksCollisio
1416 }
1417
1418 //G4cout << " currentA, Z be4: " << curre
1419 currentA += finalBaryon-initialBaryon;
1420 currentZ += finalCharge-initialCharge;
1421 //G4cout << " ApplyCollision currentA, Z
1422
1423 G4KineticTrackVector oldSecondaries;
1424 oldSecondaries.push_back(primary);
1425 primary->Hit();
1426
1427 #ifdef debug_G4BinaryCascade
1428 if ( (finalBaryon-initialBaryon-lateBaryo
1429 {
1430 G4cout << "G4BinaryCascade: Error in
1431 G4cout << "initial/final baryon numbe
1432 << initialBaryon <<" "<< fina
1433 << initialCharge <<" "<< fina
1434 << " in Collision type: "<< t
1435 << ", with number of products
1436 G4cout << G4endl<<"Initial condition
1437 G4cout << "proj: "<<collision->GetPri
1438 for(std::size_t it=0; it<collision->G
1439 {
1440 G4cout << "targ: "
1441 <<collision->GetTargetCol
1442 }
1443 PrintKTVector(&collision->GetTargetCo
1444 G4cout << G4endl<<G4endl;
1445 }
1446 #endif
1447
1448 G4KineticTrackVector oldTarget = collisio
1449 for(std::size_t ii=0; ii< oldTarget.size(
1450 {
1451 oldTarget[ii]->Hit();
1452 }
1453
1454 UpdateTracksAndCollisions(&oldSecondaries
1455 std::for_each(oldSecondaries.cbegin(), ol
1456 std::for_each(oldTarget.cbegin(), oldTarg
1457
1458 delete products;
1459 return true;
1460 }
1461
1462
1463 //-------------------------------------------
1464 G4bool G4BinaryCascade::Absorb()
1465 //-------------------------------------------
1466 {
1467 // Do it in two step: cannot change theSe
1468 G4Absorber absorber(theCutOnPAbsorb);
1469
1470 // Build the vector of G4KineticTracks th
1471 G4KineticTrackVector absorbList;
1472 std::vector<G4KineticTrack *>::const_iter
1473 // PrintKTVector(&theSecondaryList, " te
1474 for(iter = theSecondaryList.cbegin();
1475 iter != theSecondaryList.cend();
1476 {
1477 G4KineticTrack * kt = *iter;
1478 if(kt->GetState() == G4KineticTrack::
1479 {
1480 if(absorber.WillBeAbsorbed(*kt))
1481 {
1482 absorbList.push_back(kt);
1483 }
1484 }
1485 }
1486
1487 if(absorbList.empty())
1488 return false;
1489
1490 G4KineticTrackVector toDelete;
1491 for(iter = absorbList.cbegin(); iter != a
1492 {
1493 G4KineticTrack * kt = *iter;
1494 if(!absorber.FindAbsorbers(*kt, theTa
1495 throw G4HadronicException(__FILE_
1496
1497 if(!absorber.FindProducts(*kt))
1498 throw G4HadronicException(__FILE_
1499
1500 G4KineticTrackVector * products = abs
1501 G4int maxLoopCount = 1000;
1502 while(!CheckPauliPrinciple(products)
1503 {
1504 ClearAndDestroy(products);
1505 if(!absorber.FindProducts(*kt))
1506 throw G4HadronicException(__F
1507 "G4BinaryCascade::Abs
1508 }
1509 if ( --maxLoopCount < 0 ) throw G4Hadro
1510 // ------ debug
1511 // G4cerr << "Absorb CheckPauliPri
1512 // ------ end debug
1513 G4KineticTrackVector toRemove; // bu
1514 toRemove.push_back(kt);
1515 toDelete.push_back(kt); // delete th
1516 G4KineticTrackVector * absorbers = ab
1517 UpdateTracksAndCollisions(&toRemove,
1518 ClearAndDestroy(absorbers);
1519 }
1520 ClearAndDestroy(&toDelete);
1521 return true;
1522 }
1523
1524
1525
1526 // Capture all p and n with Energy < theCutOn
1527 //-------------------------------------------
1528 G4bool G4BinaryCascade::Capture(G4bool verbos
1529 //-------------------------------------------
1530 {
1531 G4KineticTrackVector captured;
1532 G4bool capture = false;
1533 std::vector<G4KineticTrack *>::const_iter
1534
1535 G4RKPropagation * RKprop=(G4RKPropagation
1536
1537 G4double capturedEnergy = 0;
1538 G4int particlesAboveCut=0;
1539 G4int particlesBelowCut=0;
1540 if ( verbose ) G4cout << " Capture: secon
1541 for(i = theSecondaryList.cbegin(); i != t
1542 {
1543 G4KineticTrack * kt = *i;
1544 if (verbose) G4cout << "Capture posit
1545 if(kt->GetState() == G4KineticTrack::
1546 {
1547 if((kt->GetDefinition() == G4Prot
1548 (kt->GetDefinition() == G
1549 {
1550 //GF cut on kinetic energy
1551 G4double field=RKprop->GetFie
1552 -RKprop->GetBarr
1553 G4double energy= kt->Get4Mome
1554 if (verbose ) G4cout << "Capt
1555 // if( energy < theCutOnP )
1556 // {
1557 capturedEnergy+=energy;
1558 ++particlesBelowCut;
1559 // } else
1560 // {
1561 // ++particlesAboveCut;
1562 // }
1563 }
1564 }
1565 }
1566 if (verbose) G4cout << "Capture particles
1567 << particlesAboveCut << " " << pa
1568 << " " << G4endl;
1569 // << " " << (particlesBelowCut>0) ? (capt
1570 // if(particlesAboveCut==0 && particlesB
1571 if(particlesBelowCut>0 && capturedEnergy/
1572 {
1573 capture=true;
1574 for(i = theSecondaryList.cbegin(); i
1575 {
1576 G4KineticTrack * kt = *i;
1577 if(kt->GetState() == G4KineticTra
1578 {
1579 if((kt->GetDefinition() == G4
1580 (kt->GetDefinition()
1581 {
1582 captured.push_back(kt);
1583 kt->Hit(); //
1584 theCapturedList.push_back
1585 }
1586 }
1587 }
1588 UpdateTracksAndCollisions(&captured,
1589 }
1590
1591 return capture;
1592 }
1593
1594
1595 //-------------------------------------------
1596 G4bool G4BinaryCascade::CheckPauliPrinciple(G
1597 //-------------------------------------------
1598 {
1599 G4int A = the3DNucleus->GetMassNumber();
1600 G4int Z = the3DNucleus->GetCharge();
1601
1602 G4FermiMomentum fermiMom;
1603 fermiMom.Init(A, Z);
1604
1605 const G4VNuclearDensity * density=the3DNu
1606
1607 G4KineticTrackVector::const_iterator i;
1608 const G4ParticleDefinition * definition;
1609
1610 // ------ debug
1611 G4bool myflag = true;
1612 // ------ end debug
1613 // G4ThreeVector xpos(0);
1614 for(i = products->cbegin(); i != products
1615 {
1616 definition = (*i)->GetDefinition();
1617 if((definition == G4Proton::Proton())
1618 (definition == G4Neutron::Neu
1619 {
1620 G4ThreeVector pos = (*i)->GetPosi
1621 G4double d = density->GetDensity(
1622 // energy correspondiing to fermi
1623 G4double eFermi = std::sqrt( sqr(
1624 if( definition == G4Proton::Proto
1625 {
1626 eFermi -= the3DNucleus->Coulo
1627 }
1628 G4LorentzVector mom = (*i)->Get4M
1629 // ------ debug
1630 /*
1631 * G4cout << "p:[" << (1/M
1632 * << (1/MeV)*mom.z()
1633 * << (1/MeV)*mom.vect
1634 * << (1/MeV)*mom.mag(
1635 * << (1/fermi)*pos.x(
1636 * << (1/fermi)*pos.z(
1637 * << (1/fermi)*(pos-x
1638 * << (1/MeV)*p << G4e
1639 * xpos=pos;
1640 */
1641 // ------ end debug
1642 if(mom.e() < eFermi )
1643 {
1644 // ------ debug
1645 myflag = false;
1646 // ------ end debug
1647 // return false;
1648 }
1649 }
1650 }
1651 #ifdef debug_BIC_CheckPauli
1652 if ( myflag )
1653 {
1654 for(i = products->cbegin(); i != prod
1655 {
1656 definition = (*i)->GetDefinition(
1657 if((definition == G4Proton::Proto
1658 (definition == G4Neutron:
1659 {
1660 G4ThreeVector pos = (*i)->Get
1661 G4double d = density->GetDens
1662 G4double pFermi = fermiMom.Ge
1663 G4LorentzVector mom = (*i)->G
1664 G4double field =((G4RKPropaga
1665 if ( mom.e()-mom.mag()+field
1666 {
1667 G4cout << "momentum probl
1668 << " mom, mom.m "
1669 << " field " << f
1670 }
1671 }
1672 }
1673 }
1674 #endif
1675
1676 return myflag;
1677 }
1678
1679 //-------------------------------------------
1680 void G4BinaryCascade::StepParticlesOut()
1681 //-------------------------------------------
1682 {
1683 G4int counter=0;
1684 G4int countreset=0;
1685 //G4cout << " nucl. Radius " << radius <<
1686 // G4cerr <<"pre-while- theSecondaryList
1687 while( theSecondaryList.size() > 0 )
1688
1689 {
1690 G4double minTimeStep = 1.e-12*ns; /
1691 /
1692 for(auto i = theSecondaryList.cbegin(
1693 {
1694 G4KineticTrack * kt = *i;
1695 if( kt->GetState() == G4KineticTr
1696 {
1697 G4double tStep(0), tdummy(0);
1698 G4bool intersect =
1699 ((G4RKPropagation*)th
1700 #ifdef debug_BIC_StepParticlesOut
1701 G4cout << " minTimeStep, tSte
1702 << " " <<kt->GetDefin
1703 << " 4mom " << kt->Ge
1704 if ( ! intersect );
1705 {
1706 PrintKTVector(&theSeconda
1707 throw G4HadronicException
1708 }
1709 #endif
1710 if(intersect && tStep<minTime
1711 {
1712 minTimeStep = tStep;
1713 }
1714 } else if ( kt->GetState() != G4K
1715 PrintKTVector(&theSecondaryLi
1716 throw G4HadronicException(__F
1717 }
1718 }
1719 minTimeStep *= 1.2;
1720 G4double timeToCollision=DBL_MAX;
1721 G4CollisionInitialState * nextCollisi
1722 if(theCollisionMgr->Entries() > 0)
1723 {
1724 nextCollision = theCollisionMgr->
1725 timeToCollision = nextCollision->
1726 // G4cout << " NextCollision *
1727 }
1728 if ( timeToCollision > minTimeStep )
1729 {
1730 DoTimeStep(minTimeStep);
1731 ++counter;
1732 } else
1733 {
1734 if (!DoTimeStep(timeToCollision)
1735 {
1736 // Check if nextCollision is
1737 if (theCollisionMgr->GetNextC
1738 {
1739 nextCollision = nullptr;
1740 }
1741 }
1742 // G4cerr <<"post- DoTimeStep 3"<
1743
1744 if(nextCollision)
1745 {
1746 if ( ApplyCollision(nextColl
1747 {
1748 // G4cout << "ApplyCollis
1749 } else
1750 {
1751 theCollisionMgr->RemoveCo
1752 }
1753 }
1754 }
1755
1756 if(countreset>100)
1757 {
1758 #ifdef debug_G4BinaryCascade
1759 G4cerr << "G4BinaryCascade.cc: Wa
1760 PrintKTVector(&theSecondaryList,"
1761 #endif
1762
1763 // add left secondaries to Final
1764 for (auto iter=theSecondaryList.c
1765 {
1766 theFinalState.push_back(*iter
1767 }
1768 theSecondaryList.clear();
1769
1770 break;
1771 }
1772
1773 if(Absorb())
1774 {
1775 // haveProducts = true;
1776 // G4cout << "Absorb sucess " <<
1777 }
1778
1779 if(Capture(false))
1780 {
1781 // haveProducts = true;
1782 #ifdef debug_BIC_StepParticlesOut
1783 G4cout << "Capture sucess " << G4
1784 #endif
1785 }
1786 if ( counter > 100 && theCollisionMgr
1787 {
1788 #ifdef debug_BIC_StepParticlesOut
1789 PrintKTVector(&theSecondaryList,s
1790 #endif
1791 FindCollisions(&theSecondaryList)
1792 counter=0;
1793 ++countreset;
1794 }
1795 //G4cout << "currentZ @ end loop " <<
1796 if ( ! currentZ ){
1797 // nucleus completely destroyed,
1798 // products = FillVoidNucleusProdu
1799 #ifdef debug_BIC_return
1800 G4cout << "return @ Z=0 after col
1801 PrintKTVector(&theSecondaryList,s
1802 G4cout << "theTargetList size: "
1803 PrintKTVector(&theTargetList,std:
1804 PrintKTVector(&theCapturedList,st
1805
1806 G4cout << " ExcitE be4 Correct :
1807 G4cout << " Mom Transfered to nuc
1808 PrintKTVector(&theFinalState,std:
1809 // G4cout << "returned products: "
1810 #endif
1811 }
1812
1813 }
1814 // G4cerr <<"Finished capture loop "<<G4
1815
1816 //G4cerr <<"pre- DoTimeStep 4"<<G4endl;
1817 DoTimeStep(DBL_MAX);
1818 //G4cerr <<"post- DoTimeStep 4"<<G4endl;
1819 }
1820
1821 //-------------------------------------------
1822 G4double G4BinaryCascade::CorrectShortlivedPr
1823 G4KineticTrack* primary,G4KineticTrac
1824 //-------------------------------------------
1825 {
1826 G4double Efermi(0);
1827 if (primary->GetState() == G4KineticTrack
1828 G4int PDGcode=primary->GetDefinition(
1829 Efermi=((G4RKPropagation *)thePropaga
1830
1831 if ( std::abs(PDGcode) > 1000 && PDGc
1832 {
1833 Efermi = ((G4RKPropagation *)theP
1834 G4LorentzVector mom4Primary=prima
1835 primary->Update4Momentum(mom4Prim
1836 }
1837
1838 for (auto titer=target_collection.cbe
1839 {
1840 const G4ParticleDefinition * aDef
1841 G4int aCode=aDef->GetPDGEncoding(
1842 G4ThreeVector aPos=(*titer)->GetP
1843 Efermi+= ((G4RKPropagation *)theP
1844 }
1845 }
1846 return Efermi;
1847 }
1848
1849 //-------------------------------------------
1850 G4bool G4BinaryCascade::CorrectShortlivedFina
1851 G4double initial_Efermi)
1852 //-------------------------------------------
1853 {
1854 G4double final_Efermi(0);
1855 G4KineticTrackVector resonances;
1856 for (auto i =products->cbegin(); i != pro
1857 {
1858 G4int PDGcode=(*i)->GetDefinition()->
1859 // G4cout << " PDGcode, state "
1860 final_Efermi+=((G4RKPropagation *)the
1861 if ( std::abs(PDGcode) > 1000 && PDGc
1862 {
1863 resonances.push_back(*i);
1864 }
1865 }
1866 if ( resonances.size() > 0 )
1867 {
1868 G4double delta_Fermi= (initial_Efermi
1869 for (auto res=resonances.cbegin(); re
1870 {
1871 G4LorentzVector mom=(*res)->Get4M
1872 G4double mass2=mom.mag2();
1873 G4double newEnergy=mom.e() + delt
1874 G4double newEnergy2= newEnergy*ne
1875 //G4cout << "mom = " << mom <<" n
1876 if ( newEnergy2 < mass2 )
1877 {
1878 return false;
1879 }
1880 G4ThreeVector mom3=std::sqrt(newE
1881 (*res)->Set4Momentum(G4LorentzVec
1882 //G4cout << " correct resonan
1883 // " 3mom from/to " <<
1884 }
1885 }
1886 return true;
1887 }
1888
1889 //-------------------------------------------
1890 void G4BinaryCascade::CorrectFinalPandE()
1891 //-------------------------------------------
1892 //
1893 // Modify momenta of outgoing particles.
1894 // Assume two body decay, nucleus(@nominal
1895 // momentum of SFSP shall be less than mome
1896 //
1897 {
1898 #ifdef debug_BIC_CorrectFinalPandE
1899 G4cerr << "BIC: -CorrectFinalPandE called
1900 #endif
1901
1902 if ( theFinalState.size() == 0 ) return;
1903
1904 G4KineticTrackVector::const_iterator i;
1905 G4LorentzVector pNucleus=GetFinal4Momentu
1906 if ( pNucleus.e() == 0 ) return; // ch
1907 #ifdef debug_BIC_CorrectFinalPandE
1908 G4cerr << " -CorrectFinalPandE 3" << G4en
1909 #endif
1910 G4LorentzVector pFinals(0);
1911 for(i = theFinalState.cbegin(); i != theF
1912 {
1913 pFinals += (*i)->Get4Momentum();
1914 #ifdef debug_BIC_CorrectFinalPandE
1915 G4cout <<"CorrectFinalPandE a final "
1916 << " 4mom " << (*i)->G
1917 #endif
1918 }
1919 #ifdef debug_BIC_CorrectFinalPandE
1920 G4cout << "CorrectFinalPandE pN pF: " <<p
1921 #endif
1922 G4LorentzVector pCM=pNucleus + pFinals;
1923
1924 G4LorentzRotation toCMS(-pCM.boostVector(
1925 pFinals *=toCMS;
1926 #ifdef debug_BIC_CorrectFinalPandE
1927 G4cout << "CorrectFinalPandE pCM, CMS pCM
1928 G4cout << "CorrectFinal CMS pN pF " <<toC
1929 <<pFinals << G4endl
1930 << " nucleus initial mass : " <<G
1931 <<" massInNucleus m(nucleus) m(fi
1932 << pFinals.mag() << " " << pCM.ma
1933 #endif
1934
1935 G4LorentzRotation toLab = toCMS.inverse()
1936
1937 G4double s0 = pCM.mag2();
1938 G4double m10 = GetIonMass(currentZ,curren
1939 G4double m20 = pFinals.mag();
1940 if( s0-(m10+m20)*(m10+m20) < 0 )
1941 {
1942 #ifdef debug_BIC_CorrectFinalPandE
1943 G4cout << "G4BinaryCascade::CorrectFi
1944
1945 G4cout << "not enough mass to correct
1946 << (s0-(m10+m20)*(m10+m20)) <
1947 << currentA << " " << current
1948 << m10 << " " << m20
1949 << G4endl;
1950 G4cerr << " -CorrectFinalPandE 4" <<
1951
1952 PrintKTVector(&theFinalState," mass p
1953 #endif
1954 return;
1955 }
1956
1957 // Three momentum in cm system
1958 G4double pInCM = std::sqrt((s0-(m10+m20)*
1959 #ifdef debug_BIC_CorrectFinalPandE
1960 G4cout <<" CorrectFinalPandE pInCM new,
1961 << " " << (pFinals).vect().mag()<
1962 #endif
1963 if ( pFinals.vect().mag() > pInCM )
1964 {
1965 G4ThreeVector p3finals=pInCM*pFinals.
1966
1967 G4double factor=std::max(0.98,pInCM/p
1968 G4LorentzVector qFinals(0);
1969 for(i = theFinalState.cbegin(); i !=
1970 {
1971 // G4ThreeVector p3((toCMS*(
1972 G4ThreeVector p3(factor*(toCMS*(*
1973 G4LorentzVector p(p3,std::sqrt((*
1974 qFinals += p;
1975 p *= toLab;
1976 #ifdef debug_BIC_CorrectFinalPandE
1977 G4cout << " final p corrected: "
1978 #endif
1979 (*i)->Set4Momentum(p);
1980 }
1981 #ifdef debug_BIC_CorrectFinalPandE
1982 G4cout << "CorrectFinalPandE nucleus
1983 <<GetFinal4Momentum().mag() <
1984 << " CMS pFinals , mag, 3.mag
1985 G4cerr << " -CorrectFinalPandE 5 " <<
1986 #endif
1987 }
1988 #ifdef debug_BIC_CorrectFinalPandE
1989 else { G4cerr << " -CorrectFinalPandE 6 -
1990 #endif
1991
1992 }
1993
1994 //-------------------------------------------
1995 void G4BinaryCascade::UpdateTracksAndCollisio
1996 //-----------------------------------
1997 G4KineticTrackVector * oldSecondaries
1998 G4KineticTrackVector * oldTarget,
1999 G4KineticTrackVector * newSecondaries
2000 {
2001 std::vector<G4KineticTrack *>::const_iter
2002
2003 // remove old secondaries from the second
2004 if(oldSecondaries)
2005 {
2006 if(!oldSecondaries->empty())
2007 {
2008 for(auto iter1=oldSecondaries->cb
2009 {
2010 iter2 = std::find(theSecondar
2011 if ( iter2 != theSecondaryLis
2012 }
2013 theCollisionMgr->RemoveTracksColl
2014 }
2015 }
2016
2017 // remove old target from the target list
2018 if(oldTarget)
2019 {
2020 // G4cout << "################## Debu
2021 if(oldTarget->size()!=0)
2022 {
2023
2024 // G4cout << "##################
2025 for(auto iter1 = oldTarget->cbegi
2026 {
2027 iter2 = std::find(theTargetLi
2028 theTargetList.erase(iter2);
2029 }
2030 theCollisionMgr->RemoveTracksColl
2031 }
2032 }
2033
2034 if(newSecondaries)
2035 {
2036 if(!newSecondaries->empty())
2037 {
2038 // insert new secondaries in the
2039 for(auto iter1 = newSecondaries->
2040 {
2041 theSecondaryList.push_back(*i
2042 if ((*iter1)->GetState() == G
2043 {
2044 PrintKTVector(*iter1, "und
2045 }
2046
2047
2048 }
2049 // look for collisions of new sec
2050 FindCollisions(newSecondaries);
2051 }
2052 }
2053 // G4cout << "Exiting ... "<<oldTarget<<G
2054 }
2055
2056
2057 class SelectFromKTV
2058 {
2059 private:
2060 G4KineticTrackVector * ktv;
2061 G4KineticTrack::CascadeState wanted_state
2062 public:
2063 SelectFromKTV(G4KineticTrackVector * out,
2064 :
2065 ktv(out), wanted_state(astate)
2066 {};
2067 void operator() (G4KineticTrack *& kt) co
2068 {
2069 if ( (kt)->GetState() == wanted_state
2070 };
2071 };
2072
2073
2074
2075 //-------------------------------------------
2076 G4bool G4BinaryCascade::DoTimeStep(G4double t
2077 //-------------------------------------------
2078 {
2079
2080 #ifdef debug_BIC_DoTimeStep
2081 G4ping debug("debug_G4BinaryCascade");
2082 debug.push_back("======> DoTimeStep 1");
2083 G4cerr <<"G4BinaryCascade::DoTimeStep: en
2084 << " , time="<<theCurrentTime <<
2085 PrintKTVector(&theSecondaryList, std::str
2086 //PrintKTVector(&theTargetList, std::stri
2087 #endif
2088
2089 G4bool success=true;
2090 std::vector<G4KineticTrack *>::const_iter
2091
2092 G4KineticTrackVector * kt_outside = new G
2093 std::for_each( theSecondaryList.begin(),t
2094 SelectFromKTV(kt_outside,G4Kineti
2095 //PrintKTVector(kt_outside, std::string("
2096
2097 G4KineticTrackVector * kt_inside = new G4
2098 std::for_each( theSecondaryList.begin(),t
2099 SelectFromKTV(kt_inside, G4Kineti
2100 // PrintKTVector(kt_inside, std::string(
2101 //-----
2102 G4KineticTrackVector dummy; // needed f
2103 #ifdef debug_BIC_DoTimeStep
2104 G4cout << "NOW WE ARE ENTERING THE TRANSP
2105 #endif
2106
2107 // =================== Here we move the p
2108
2109 thePropagator->Transport(theSecondaryList
2110
2111 // =================== Here we move the p
2112
2113 //------
2114
2115 theMomentumTransfer += thePropagator->Get
2116 #ifdef debug_BIC_DoTimeStep
2117 G4cout << "DoTimeStep : theMomentumTransf
2118 PrintKTVector(&theSecondaryList, std::str
2119 #endif
2120
2121 //_DebugEpConservation(" after stepping")
2122
2123 // Partclies which went INTO nucleus
2124
2125 G4KineticTrackVector * kt_gone_in = new G
2126 std::for_each( kt_outside->begin(),kt_out
2127 SelectFromKTV(kt_gone_in,G4Kineti
2128 // PrintKTVector(kt_gone_in, std::string
2129
2130
2131 // Partclies which went OUT OF nucleus
2132 G4KineticTrackVector * kt_gone_out = new
2133 std::for_each( kt_inside->begin(),kt_insi
2134 SelectFromKTV(kt_gone_out, G4Kine
2135
2136 // PrintKTVector(kt_gone_out, std::strin
2137
2138 G4KineticTrackVector *fail=CorrectBarions
2139
2140 if ( fail )
2141 {
2142 // some particle(s) supposed to enter
2143 kt_gone_in->clear();
2144 std::for_each( kt_outside->begin(),kt
2145 SelectFromKTV(kt_gone_in,G4Ki
2146
2147 kt_gone_out->clear();
2148 std::for_each( kt_inside->begin(),kt_
2149 SelectFromKTV(kt_gone_out, G4
2150
2151 #ifdef debug_BIC_DoTimeStep
2152 PrintKTVector(fail,std::string(" Fail
2153 PrintKTVector(kt_gone_in, std::string
2154 PrintKTVector(kt_gone_out, std::strin
2155 #endif
2156 delete fail;
2157 }
2158
2159 // Add tracks missing nucleus and tracks
2160 std::for_each( kt_outside->begin(),kt_out
2161 SelectFromKTV(kt_gone_out,G4Kinet
2162 //PrintKTVector(kt_gone_out, std::string(
2163 std::for_each( kt_outside->begin(),kt_out
2164 SelectFromKTV(kt_gone_out,G4Kinet
2165
2166 #ifdef debug_BIC_DoTimeStep
2167 PrintKTVector(kt_gone_out, std::string("a
2168 #endif
2169
2170 theFinalState.insert(theFinalState.end(),
2171 kt_gone_out->begin(),kt_gone_out-
2172
2173 // Partclies which could not leave nucleu
2174 G4KineticTrackVector * kt_captured = new
2175 std::for_each( theSecondaryList.begin(),t
2176 SelectFromKTV(kt_captured, G4Kine
2177
2178 // Check no track is part in next collisi
2179 // this step was to far, and collisions
2180
2181 if ( theCollisionMgr->Entries()> 0 )
2182 {
2183 if (kt_gone_out->size() )
2184 {
2185 G4KineticTrack * nextPrimary = th
2186 iter = std::find(kt_gone_out->beg
2187 if ( iter != kt_gone_out->cend()
2188 {
2189 success=false;
2190 #ifdef debug_BIC_DoTimeStep
2191 G4cout << " DoTimeStep - WARN
2192 #endif
2193 }
2194 }
2195 if ( kt_captured->size() )
2196 {
2197 G4KineticTrack * nextPrimary = th
2198 iter = std::find(kt_captured->beg
2199 if ( iter != kt_captured->cend()
2200 {
2201 success=false;
2202 #ifdef debug_BIC_DoTimeStep
2203 G4cout << " DoTimeStep - WARN
2204 #endif
2205 }
2206 }
2207
2208 }
2209 // PrintKTVector(kt_gone_out," kt_gone_ou
2210 UpdateTracksAndCollisions(kt_gone_out,0 ,
2211
2212
2213 if ( kt_captured->size() )
2214 {
2215 theCapturedList.insert(theCapturedLis
2216 kt_captured->begin(),kt_captu
2217 //should be std::for_each(kt_cap
2218 // std::mem_fun(&G4Kinet
2219 // but VC 6 requires:
2220 for(auto i_captured=kt_captured->cbeg
2221 {
2222 (*i_captured)->Hit();
2223 }
2224 // PrintKTVector(kt_captured," kt
2225 UpdateTracksAndCollisions(kt_captured
2226 }
2227
2228 #ifdef debug_G4BinaryCascade
2229 delete kt_inside;
2230 kt_inside = new G4KineticTrackVector;
2231 std::for_each( theSecondaryList.begin(),t
2232 SelectFromKTV(kt_inside, G4Kineti
2233 if ( currentZ != (GetTotalCharge(theTarge
2234 + GetTotalCharge(theCapturedList)
2235 + GetTotalCharge(*kt_inside)) )
2236 {
2237 G4cout << " error-DoTimeStep aft, A,
2238 << " sum(tgt,capt,active) "
2239 << GetTotalCharge(theTargetLi
2240 << " targets: " << GetTotalC
2241 << " captured: " << GetTotalC
2242 << " active: " << GetTotalC
2243 << G4endl;
2244 }
2245 #endif
2246
2247 delete kt_inside;
2248 delete kt_outside;
2249 delete kt_captured;
2250 delete kt_gone_in;
2251 delete kt_gone_out;
2252
2253 // G4cerr <<"G4BinaryCascade::DoTimeStep
2254 theCurrentTime += theTimeStep;
2255
2256 //debug.push_back("======> DoTimeStep 2")
2257 return success;
2258
2259 }
2260
2261 //-------------------------------------------
2262 G4KineticTrackVector* G4BinaryCascade::Correc
2263 G4KineticTrackVector *in,
2264 G4KineticTrackVector *out)
2265 //-------------------------------------------
2266 {
2267 G4KineticTrackVector * kt_fail(nullptr);
2268 std::vector<G4KineticTrack *>::const_iter
2269 // G4cout << "CorrectBarionsOnBoundary,c
2270 // << currentZ << " "<< currentA
2271 if (in->size())
2272 {
2273 G4int secondaries_in(0);
2274 G4int secondaryBarions_in(0);
2275 G4int secondaryCharge_in(0);
2276 G4double secondaryMass_in(0);
2277
2278 for ( iter =in->cbegin(); iter != in-
2279 {
2280 ++secondaries_in;
2281 secondaryCharge_in += G4lrint((*i
2282 if ((*iter)->GetDefinition()->Get
2283 {
2284 secondaryBarions_in += (*iter
2285 if((*iter)->GetDefinition() =
2286 (*iter)->GetDefinitio
2287 {
2288 secondaryMass_in += (*ite
2289 } else {
2290 secondaryMass_in += G4Pro
2291 }
2292 }
2293 }
2294 G4double mass_initial= GetIonMass(cur
2295
2296 currentZ += secondaryCharge_in;
2297 currentA += secondaryBarions_in;
2298
2299 // G4cout << "CorrectBarionsOnBounda
2300 // << secondaryCharge_in <
2301
2302 G4double mass_final= GetIonMass(curre
2303
2304 G4double correction= secondaryMass_in
2305 if (secondaries_in>1)
2306 {correction /= secondaries_in;}
2307
2308 #ifdef debug_BIC_CorrectBarionsOnBoundary
2309 G4cout << "CorrectBarionsOnBoundary,c
2310 << "secondaryCharge_in,second
2311 << "energy correction,m_secon
2312 << currentZ << " "<< currentA
2313 << secondaryCharge_in<<" "<<s
2314 << correction << " "
2315 << secondaryMass_in << " "
2316 << mass_initial << " "
2317 << mass_final << " "
2318 << G4endl;
2319 PrintKTVector(in,std::string("in be4
2320 #endif
2321 for ( iter = in->cbegin(); iter != in
2322 {
2323 if ((*iter)->GetTrackingMomentum(
2324 {
2325 (*iter)->UpdateTrackingMoment
2326 } else {
2327 //particle cannot go in, put
2328 G4RKPropagation * RKprop=(G4R
2329 (*iter)->SetState(G4KineticTr
2330 // Undo correction for Colomb
2331 G4double barrier=RKprop->GetB
2332 (*iter)->UpdateTrackingMoment
2333 if ( ! kt_fail ) kt_fail=new
2334 kt_fail->push_back(*iter);
2335 currentZ -= G4lrint((*iter)->
2336 currentA -= (*iter)->GetDefin
2337
2338 }
2339
2340 }
2341 #ifdef debug_BIC_CorrectBarionsOnBoundary
2342 G4cout << " CorrectBarionsOnBoundary,
2343 << currentZ << " " << current
2344 << secondaryCharge_in << " "
2345 << secondaryMass_in << " "
2346 << G4endl;
2347 PrintKTVector(in,std::string("in AFT
2348 #endif
2349
2350 }
2351 //---------------------------------------
2352 if (out->size())
2353 {
2354 G4int secondaries_out(0);
2355 G4int secondaryBarions_out(0);
2356 G4int secondaryCharge_out(0);
2357 G4double secondaryMass_out(0);
2358
2359 for ( iter = out->cbegin(); iter != o
2360 {
2361 ++secondaries_out;
2362 secondaryCharge_out += G4lrint((*
2363 if ((*iter)->GetDefinition()->Get
2364 {
2365 secondaryBarions_out += (*ite
2366 if((*iter)->GetDefinition() =
2367 (*iter)->GetDefinitio
2368 {
2369 secondaryMass_out += (*it
2370 } else {
2371 secondaryMass_out += G4Ne
2372 }
2373 }
2374 }
2375
2376 G4double mass_initial= GetIonMass(cu
2377 currentA -=secondaryBarions_out;
2378 currentZ -=secondaryCharge_out;
2379
2380 // G4cout << "CorrectBarionsOnBounda
2381 // << secondaryCharge_out
2382
2383 // a delta min
2384 // if (currentA < 0 || currentZ <
2385 if (currentA < 0 )
2386 {
2387 G4cerr << "G4BinaryCascade - seco
2388 secondaryBarions_out << "
2389 PrintKTVector(&theTargetList,"Cor
2390 PrintKTVector(&theCapturedList,"C
2391 PrintKTVector(&theSecondaryList,"
2392 G4cerr << "G4BinaryCascade - curr
2393 throw G4HadronicException(__FILE_
2394 }
2395 G4double mass_final=GetIonMass(curren
2396 G4double correction= mass_initial - m
2397 // G4cout << "G4BinaryCascade::Correc
2398
2399 if (secondaries_out>1) correction /=
2400 #ifdef debug_BIC_CorrectBarionsOnBoundary
2401 G4cout << "DoTimeStep,(current Z,A),"
2402 << "(secondaries out,Charge,B
2403 <<"* energy correction,(m_sec
2404 << "("<< currentZ << ","<< cu
2405 << secondaries_out << ","
2406 << secondaryCharge_out<<","<<
2407 << correction << " ("
2408 << secondaryMass_out << ", "
2409 << mass_initial << ", "
2410 << mass_final << ")"
2411 << G4endl;
2412 PrintKTVector(out,std::string("out be
2413 #endif
2414
2415 for ( iter = out->cbegin(); iter != o
2416 {
2417 if ((*iter)->GetTrackingMomentum(
2418 {
2419 (*iter)->UpdateTrackingMoment
2420 } else
2421 {
2422 // particle cannot go out due
2423 // capture protons and neutr
2424 if(((*iter)->GetDefinition()
2425 ((*iter)->GetDefiniti
2426 {
2427 (*iter)->SetState(G4Kinet
2428 // Undo correction for Co
2429 G4double barrier=((G4RKPr
2430 (*iter)->UpdateTrackingMo
2431 if ( kt_fail == 0 ) kt_fa
2432 kt_fail->push_back(*iter)
2433 currentZ += G4lrint((*ite
2434 currentA += (*iter)->GetD
2435 }
2436 #ifdef debug_BIC_CorrectBarionsOnBoundary
2437 else
2438 {
2439 G4cout << "Not correcting
2440 << (*iter)->GetDe
2441 << (*iter)->GetDe
2442 PrintKTVector(out,std::st
2443 }
2444 #endif
2445 }
2446 }
2447
2448 #ifdef debug_BIC_CorrectBarionsOnBoundary
2449 PrintKTVector(out,std::string("out AF
2450 G4cout << " DoTimeStep, nucl-update,
2451 << currentA << " "<< currentZ
2452 << secondaryCharge_out << " "
2453 secondaryMass_out << " "
2454 << massInNucleus << " "
2455 << GetIonMass(currentZ,curren
2456 << " " << massInNucleus - Get
2457 << G4endl;
2458 #endif
2459 }
2460
2461 return kt_fail;
2462 }
2463
2464
2465 //-------------------------------------------
2466
2467 G4Fragment * G4BinaryCascade::FindFragments()
2468 //-------------------------------------------
2469 {
2470
2471 #ifdef debug_BIC_FindFragments
2472 G4cout << "target, captured, secondary: "
2473 << theTargetList.size() << " "
2474 << theCapturedList.size()<< " "
2475 << theSecondaryList.size()
2476 << G4endl;
2477 #endif
2478
2479 G4int a = G4int(theTargetList.size()+theC
2480 G4int zTarget = 0;
2481 for(auto i = theTargetList.cbegin(); i !=
2482 {
2483 if(G4lrint((*i)->GetDefinition()->Get
2484 {
2485 zTarget++;
2486 }
2487 }
2488
2489 G4int zCaptured = 0;
2490 G4LorentzVector CapturedMomentum(0.,0.,0.
2491 for(auto i = theCapturedList.cbegin(); i
2492 {
2493 CapturedMomentum += (*i)->Get4Momentu
2494 if(G4lrint((*i)->GetDefinition()->Get
2495 {
2496 zCaptured++;
2497 }
2498 }
2499
2500 G4int z = zTarget+zCaptured;
2501
2502 #ifdef debug_G4BinaryCascade
2503 if ( z != (GetTotalCharge(theTargetList)
2504 {
2505 G4cout << " FindFragment Counting err
2506 << GetTotalCharge(theTargetLi
2507 G4endl;
2508 PrintKTVector(&theTargetList, std::st
2509 PrintKTVector(&theCapturedList, std::
2510 }
2511 #endif
2512 //debug
2513 /*
2514 * G4cout << " Fragment mass table / re
2515 * << GetIonMass(z, a)
2516 * << " / " << GetFinal4Momentum().mag(
2517 * << " difference "
2518 * << GetFinal4Momentum().mag() - GetI
2519 * << G4endl;
2520 */
2521 //
2522 // if(fBCDEBUG) G4cerr << "Fragment A, Z
2523 if ( z < 1 ) return 0;
2524
2525 G4int holes = G4int(the3DNucleus->GetMass
2526 G4int excitons = (G4int)theCapturedList.s
2527 #ifdef debug_BIC_FindFragments
2528 G4cout << "Fragment: a= " << a << " z= "
2529 << " Charged= " << zCaptured << "
2530 << " excitE= " <<GetExcitationEne
2531 << " Final4Momentum= " << GetFina
2532 << G4endl;
2533 #endif
2534
2535 G4Fragment * fragment = new G4Fragment(a,
2536 fragment->SetNumberOfHoles(holes);
2537
2538 //GF fragment->SetNumberOfParticles(exci
2539 fragment->SetNumberOfParticles(excitons);
2540 fragment->SetNumberOfCharged(zCaptured);
2541 fragment->SetCreatorModelID(theBIC_ID);
2542
2543 return fragment;
2544 }
2545
2546 //-------------------------------------------
2547
2548 G4LorentzVector G4BinaryCascade::GetFinal4Mom
2549 //-------------------------------------------
2550 // Return momentum of reminder nulceus;
2551 // ie. difference of (initial state(primarie
2552 {
2553 G4LorentzVector final4Momentum = theIniti
2554 G4LorentzVector finals(0,0,0,0);
2555 for(auto i = theFinalState.cbegin(); i !=
2556 {
2557 final4Momentum -= (*i)->Get4Momentum(
2558 finals += (*i)->Get4Momentum();
2559 }
2560
2561 if(final4Momentum.e()> 0 && (final4Moment
2562 {
2563 #ifdef debug_BIC_Final4Momentum
2564 G4cerr << G4endl;
2565 G4cerr << "G4BinaryCascade::GetFinal4
2566 G4KineticTrackVector::iterator i;
2567 G4cerr <<"Total initial 4-momentum "
2568 G4cerr <<" GetFinal4Momentum: Initial
2569 for(i = theFinalState.begin(); i != t
2570 {
2571 G4cerr <<" Final state: "<<(*i)->
2572 }
2573 G4cerr << "Sum( 4-mom ) finals " << f
2574 G4cerr<< " Final4Momentum = "<<final4
2575 G4cerr <<" current A, Z = "<< current
2576 G4cerr << G4endl;
2577 #endif
2578
2579 final4Momentum=G4LorentzVector(0,0,0,
2580 }
2581 return final4Momentum;
2582 }
2583
2584 //-------------------------------------------
2585 G4LorentzVector G4BinaryCascade::GetFinalNucl
2586 //-------------------------------------------
2587 {
2588 // return momentum of nucleus for use wit
2589 // apply to precompoud products.
2590
2591 G4LorentzVector CapturedMomentum(0,0,0,0)
2592 // G4cout << "GetFinalNucleusMomentum Ca
2593 for(auto i = theCapturedList.cbegin(); i
2594 {
2595 CapturedMomentum += (*i)->Get4Momentu
2596 }
2597 //G4cout << "GetFinalNucleusMomentum Capt
2598 // G4cerr << "it 9"<<G4endl;
2599
2600 G4LorentzVector NucleusMomentum = GetFina
2601 if ( NucleusMomentum.e() > 0 )
2602 {
2603 // G4cout << "GetFinalNucleusMomentum
2604 // boost nucleus to a frame such that
2605 G4ThreeVector boost= (NucleusMomentum
2606 if(boost.mag2()>1.0)
2607 {
2608 # ifdef debug_BIC_FinalNucleusMomentum
2609 G4cerr << "G4BinaryCascade::GetFi
2610 G4cerr << "it 0"<<boost <<G4endl;
2611 G4cerr << "it 01"<<NucleusMomentu
2612 G4cout <<" testing boost "<<boost
2613 # endif
2614 boost=G4ThreeVector(0,0,0);
2615 NucleusMomentum=G4LorentzVector(0
2616 }
2617 G4LorentzRotation nucleusBoost( -boo
2618 precompoundLorentzboost.set( boost );
2619 #ifdef debug_debug_BIC_FinalNucleusMomentum
2620 G4cout << "GetFinalNucleusMomentum be
2621 #endif
2622 NucleusMomentum *= nucleusBoost;
2623 #ifdef debug_BIC_FinalNucleusMomentum
2624 G4cout << "GetFinalNucleusMomentum af
2625 #endif
2626 }
2627 return NucleusMomentum;
2628 }
2629
2630 //-------------------------------------------
2631 G4ReactionProductVector * G4BinaryCascade::Pr
2632 //-----------------------------------
2633 G4KineticTrackVector * secondaries, G
2634 {
2635 G4ReactionProductVector * products = new
2636 const G4ParticleDefinition * aHTarg = G4P
2637 if (nucleus->GetCharge() == 0) aHTarg = G
2638 G4double mass = aHTarg->GetPDGMass();
2639 G4KineticTrackVector * secs = nullptr;
2640 G4ThreeVector pos(0,0,0);
2641 G4LorentzVector mom(mass);
2642 G4KineticTrack aTarget(aHTarg, 0., pos, m
2643 G4bool done(false);
2644 // data member static G4Scatterer theH
2645 //G4cout << " start 1H1 for " << (*second
2646 // << " on " << aHTarg->GetParticle
2647 G4int tryCount(0);
2648 while(!done && tryCount++ <200)
2649 {
2650 if(secs)
2651 {
2652 std::for_each(secs->begin(), secs
2653 delete secs;
2654 }
2655 secs = theH1Scatterer->Scatter(*(*sec
2656 #ifdef debug_H1_BinaryCascade
2657 PrintKTVector(secs," From Scatter");
2658 #endif
2659 for(std::size_t ss=0; secs && ss<secs
2660 {
2661 // must have one resonance in fin
2662 if((*secs)[ss]->GetDefinition()->
2663 }
2664 }
2665
2666 ClearAndDestroy(&theFinalState);
2667 ClearAndDestroy(secondaries);
2668 delete secondaries;
2669
2670 for(std::size_t current=0; secs && curren
2671 {
2672 if((*secs)[current]->GetDefinition()-
2673 {
2674 done = true; // must have one re
2675 G4KineticTrackVector * dec = (*se
2676 for(auto jter=dec->cbegin(); jter
2677 {
2678 //G4cout << "Decay"<<G4endl;
2679 secs->push_back(*jter);
2680 //G4cout << "decay "<<(*jter)
2681 }
2682 delete (*secs)[current];
2683 delete dec;
2684 }
2685 else
2686 {
2687 //G4cout << "FS "<<G4endl;
2688 //G4cout << "FS "<<(*secs)[curren
2689 theFinalState.push_back((*secs)[c
2690 }
2691 }
2692
2693 delete secs;
2694 #ifdef debug_H1_BinaryCascade
2695 PrintKTVector(&theFinalState," FinalState
2696 #endif
2697 for(auto iter = theFinalState.cbegin(); i
2698 {
2699 G4KineticTrack * kt = *iter;
2700 G4ReactionProduct * aNew = new G4Reac
2701 aNew->SetMomentum(kt->Get4Momentum().
2702 aNew->SetTotalEnergy(kt->Get4Momentum
2703 aNew->SetCreatorModelID(theBIC_ID);
2704 aNew->SetParentResonanceDef(kt->GetPa
2705 aNew->SetParentResonanceID(kt->GetPar
2706 products->push_back(aNew);
2707 #ifdef debug_H1_BinaryCascade
2708 if (! kt->GetDefinition()->GetPDGStab
2709 {
2710 if (kt->GetDefinition()->IsShortL
2711 {
2712 G4cout << "final shortlived :
2713 } else
2714 {
2715 G4cout << "final un stable :
2716 }
2717 G4cout <<kt->GetDefinition()->Get
2718 }
2719 #endif
2720 delete kt;
2721 }
2722 theFinalState.clear();
2723 return products;
2724
2725 }
2726
2727 //-------------------------------------------
2728 G4ThreeVector G4BinaryCascade::GetSpherePoint
2729 G4double r, const G4LorentzVector & m
2730 //-------------------------------------------
2731 {
2732 // Get a point outside radius.
2733 // point is random in plane (circle o
2734 // plus -1*r*mom->vect()->unit();
2735 G4ThreeVector o1, o2;
2736 G4ThreeVector mom = mom4.vect();
2737
2738 o1= mom.orthogonal(); // we simply need
2739 o2= mom.cross(o1); // o2 is now orth
2740
2741 G4double x2, x1;
2742
2743 do
2744 {
2745 x1=(G4UniformRand()-.5)*2;
2746 x2=(G4UniformRand()-.5)*2;
2747 } while (sqr(x1) +sqr(x2) > 1.);
2748
2749 return G4ThreeVector(r*(x1*o1.unit() + x2
2750
2751
2752
2753 /*
2754 * // Get a point uniformly distributed o
2755 * // with z < 0.
2756 * G4double b = r*G4UniformRand(); //
2757 * G4double phi = G4UniformRand()*2*pi;
2758 * G4double x = b*std::cos(phi);
2759 * G4double y = b*std::sin(phi);
2760 * G4double z = -std::sqrt(r*r-b*b);
2761 * z *= 1.001; // Get position a little
2762 * point.setX(x);
2763 * point.setY(y);
2764 * point.setZ(z);
2765 */
2766 }
2767
2768 //-------------------------------------------
2769
2770 void G4BinaryCascade::ClearAndDestroy(G4Kinet
2771 //-------------------------------------------
2772 {
2773 for(auto i = ktv->cbegin(); i != ktv->cen
2774 delete (*i);
2775 ktv->clear();
2776 }
2777
2778 //-------------------------------------------
2779 void G4BinaryCascade::ClearAndDestroy(G4React
2780 //-------------------------------------------
2781 {
2782 for(auto i = rpv->cbegin(); i != rpv->cen
2783 delete (*i);
2784 rpv->clear();
2785 }
2786
2787 //-------------------------------------------
2788 void G4BinaryCascade::PrintKTVector(G4Kinetic
2789 //-------------------------------------------
2790 {
2791 if (comment.size() > 0 ) G4cout << "G4Bin
2792 if (ktv) {
2793 G4cout << " vector: " << ktv << ", n
2794 << G4endl;
2795 std::vector<G4KineticTrack *>::const_
2796 G4int count;
2797
2798 for(count = 0, i = ktv->cbegin(); i !
2799 {
2800 G4KineticTrack * kt = *i;
2801 G4cout << " track n. " << count;
2802 PrintKTVector(kt);
2803 }
2804 } else {
2805 G4cout << "G4BinaryCascade::PrintKTVe
2806 }
2807 }
2808 //-------------------------------------------
2809 void G4BinaryCascade::PrintKTVector(G4Kinetic
2810 //-------------------------------------------
2811 {
2812 if (comment.size() > 0 ) G4cout << "G4Bin
2813 if ( kt ){
2814 G4cout << ", id: " << kt << G4endl;
2815 G4ThreeVector pos = kt->GetPosition()
2816 G4LorentzVector mom = kt->Get4Momentu
2817 G4LorentzVector tmom = kt->GetTrackin
2818 const G4ParticleDefinition * definiti
2819 G4cout << " definition: " << defin
2820 << 1/fermi*pos << " R: " << 1
2821 << 1/MeV*mom <<"Tr_mom" << 1
2822 << " M: " << 1/MeV*mom.mag()
2823 G4cout <<" trackstatus: "<<kt->Get
2824 } else {
2825 G4cout << "G4BinaryCascade::PrintKTVe
2826 }
2827 }
2828
2829
2830 //-------------------------------------------
2831 G4double G4BinaryCascade::GetIonMass(G4int Z,
2832 //-------------------------------------------
2833 {
2834 G4double mass(0);
2835 if ( Z > 0 && A >= Z )
2836 {
2837 mass = G4ParticleTable::GetParticleTa
2838
2839 } else if ( A > 0 && Z>0 )
2840 {
2841 // charge Z > A; will happen for ligh
2842 mass = G4ParticleTable::GetParticleTa
2843
2844 } else if ( A >= 0 && Z<=0 )
2845 {
2846 // all neutral, or empty nucleus
2847 mass = A * G4Neutron::Neutron()->GetP
2848
2849 } else if ( A == 0 )
2850 {
2851 // empty nucleus, except maybe pions
2852 mass = 0;
2853
2854 } else
2855 {
2856 G4cerr << "G4BinaryCascade::GetIonMas
2857 << A << "," << Z << ")" <<G4e
2858 throw G4HadronicException(__FILE__, _
2859
2860 }
2861 // G4cout << "G4BinaryCascade::GetIonMass
2862 return mass;
2863 }
2864 G4ReactionProductVector * G4BinaryCascade::Fi
2865 {
2866 // return product when nucleus is destroy
2867 G4double Esecondaries(0.);
2868 G4LorentzVector psecondaries;
2869 std::vector<G4KineticTrack *>::const_iter
2870 std::vector<G4ReactionProduct *>::const_i
2871 decayKTV.Decay(&theFinalState);
2872
2873 for(iter = theFinalState.cbegin(); iter !
2874 {
2875 G4ReactionProduct * aNew = new G4Reac
2876 aNew->SetMomentum((*iter)->Get4Moment
2877 aNew->SetTotalEnergy((*iter)->Get4Mom
2878 aNew->SetCreatorModelID(theBIC_ID);
2879 aNew->SetParentResonanceDef((*iter)->GetPar
2880 aNew->SetParentResonanceID((*iter)->GetPare
2881 Esecondaries +=(*iter)->Get4Momentum(
2882 psecondaries +=(*iter)->Get4Momentum(
2883 aNew->SetNewlyAdded(true);
2884 //G4cout << " Particle Ekin " << aNew
2885 products->push_back(aNew);
2886 }
2887
2888 // pull out late particles from collision
2889 //theCollisionMgr->Print();
2890 while(theCollisionMgr->Entries() > 0)
2891 {
2892 G4CollisionInitialState *
2893 collision = theCollisionMgr->GetNextC
2894
2895 if ( ! collision->GetTargetCollection
2896 G4KineticTrackVector * lates = co
2897 if ( lates->size() == 1 ) {
2898 G4KineticTrack * atrack=*(lat
2899 //PrintKTVector(atrack, " lat
2900
2901 G4ReactionProduct * aNew = ne
2902 aNew->SetMomentum(atrack->Get
2903 aNew->SetTotalEnergy(atrack->
2904 aNew->SetCreatorModelID(atrac
2905 aNew->SetParentResonanceDef(atrack->GetPa
2906 aNew->SetParentResonanceID(atrack->GetPar
2907 Esecondaries +=atrack->Get4Mo
2908 psecondaries +=atrack->Get4Mo
2909 aNew->SetNewlyAdded(true);
2910 products->push_back(aNew);
2911
2912 }
2913 }
2914 theCollisionMgr->RemoveCollision(coll
2915
2916 }
2917
2918 // decay must be after loop on Collisions
2919 // to by Collisions.
2920 decayKTV.Decay(&theSecondaryList);
2921
2922 // Correct for momentum transfered to Nuc
2923 G4ThreeVector transferCorrection(0);
2924 if ( (theSecondaryList.size() + theCaptur
2925 {
2926 transferCorrection= theMomentumTransfer
2927 }
2928
2929 for(iter = theSecondaryList.cbegin(); ite
2930 {
2931 G4ReactionProduct * aNew = new G4Reac
2932 (*iter)->Update4Momentum((*iter)->Get
2933 aNew->SetMomentum((*iter)->Get4Moment
2934 aNew->SetTotalEnergy((*iter)->Get4Mom
2935 aNew->SetCreatorModelID(theBIC_ID);
2936 aNew->SetParentResonanceDef((*iter)->GetPar
2937 aNew->SetParentResonanceID((*iter)->GetPare
2938 Esecondaries +=(*iter)->Get4Momentum(
2939 psecondaries +=(*iter)->Get4Momentum(
2940 if ( (*iter)->IsParticipant() ) aNew-
2941 products->push_back(aNew);
2942 }
2943
2944 for(iter = theCapturedList.cbegin(); iter
2945 {
2946 G4ReactionProduct * aNew = new G4Reac
2947 (*iter)->Update4Momentum((*iter)->Get
2948 aNew->SetMomentum((*iter)->Get4Moment
2949 aNew->SetTotalEnergy((*iter)->Get4Mom
2950 aNew->SetCreatorModelID(theBIC_ID);
2951 aNew->SetParentResonanceDef((*iter)->GetPar
2952 aNew->SetParentResonanceID((*iter)->GetPare
2953 Esecondaries +=(*iter)->Get4Momentum(
2954 psecondaries +=(*iter)->Get4Momentum(
2955 aNew->SetNewlyAdded(true);
2956 products->push_back(aNew);
2957 }
2958
2959 G4double SumMassNucleons(0.);
2960 G4LorentzVector pNucleons(0.);
2961 for(iter = theTargetList.cbegin(); iter !
2962 {
2963 SumMassNucleons += (*iter)->GetDefini
2964 pNucleons += (*iter)->Get4Momentum();
2965 }
2966
2967 G4double Ekinetic=theProjectile4Momentum.
2968 #ifdef debug_BIC_FillVoidnucleus
2969 G4LorentzVector deltaP=theProjectile4
2970 psecondaries - pNucle
2971 //G4cout << "BIC::FillVoidNucleus() n
2972 // ", deltaP " << deltaP << " de
2973 #endif
2974 if (Ekinetic > 0. && theTargetList.size()
2975 Ekinetic /= theTargetList.size();
2976 } else {
2977 G4double Ekineticrdm(0);
2978 if (theTargetList.size()) Ekineticrdm
2979 G4double TotalEkin(Ekineticrdm);
2980 for (rpiter=products->cbegin(); rpite
2981 TotalEkin+=(*rpiter)->GetKineticEne
2982 }
2983 G4double correction(1.);
2984 if ( std::abs(Ekinetic) < 20*perCent
2985 correction=1. + (Ekinetic-Ekineticr
2986 }
2987 #ifdef debug_G4BinaryCascade
2988 else {
2989 G4cout << "BLIC::FillVoidNucleus() fa
2990 }
2991 #endif
2992
2993 for (rpiter=products->cbegin(); rpite
2994 (*rpiter)->SetKineticEnergy((*rpiter)
2995 (*rpiter)->SetMomentum((*rpiter)->G
2996
2997 }
2998
2999 Ekinetic=Ekineticrdm*correction;
3000 if (theTargetList.size())Ekinetic /=
3001
3002 }
3003
3004 for(iter = theTargetList.cbegin(); iter !
3005 // set Nucleon it to be hit - as it is
3006 (*iter)->Hit();
3007 G4ReactionProduct * aNew = new G4Reac
3008 aNew->SetKineticEnergy(Ekinetic);
3009 aNew->SetMomentum(aNew->GetTotalMomen
3010 aNew->SetNewlyAdded(true);
3011 aNew->SetCreatorModelID(theBIC_ID);
3012 aNew->SetParentResonanceDef((*iter)->GetPar
3013 aNew->SetParentResonanceID((*iter)->GetPare
3014 products->push_back(aNew);
3015 Esecondaries += aNew->GetTotalEnergy(
3016 psecondaries += G4LorentzVector(aNew-
3017 }
3018 psecondaries=G4LorentzVector(0);
3019 for (rpiter=products->cbegin(); rpiter!=p
3020 psecondaries += G4LorentzVector((*rpite
3021 }
3022
3023 G4LorentzVector initial4Mom=theProjectile
3024
3025 //G4cout << "::FillVoidNucleus()final e/p
3026 // << " final " << psecondaries << " del
3027
3028 G4ThreeVector SumMom=psecondaries.vect();
3029
3030 SumMom=initial4Mom.vect()-SumMom;
3031 G4int loopcount(0);
3032
3033 // reverse_iterator reverse - start to co
3034 while ( SumMom.mag() > 0.1*MeV && loopcou
3035 {
3036 G4int index=(G4int)products->size();
3037 for (auto reverse=products->crbegin();
3038 SumMom=initial4Mom.vect();
3039 for (rpiter=products->cbegin(); rpite
3040 SumMom-=(*rpiter)->GetMomentum();
3041 }
3042 G4double p=((*reverse)->GetMomentum()
3043 (*reverse)->SetMomentum( p*(((*rever
3044 }
3045 }
3046 return products;
3047 }
3048
3049 G4ReactionProductVector * G4BinaryCascade::Hi
3050 G4KineticTrackVector * secondaries)
3051 {
3052 for(auto iter = secondaries->cbegin(); it
3053 {
3054 G4ReactionProduct * aNew = new G4Reac
3055 aNew->SetMomentum((*iter)->Get4Moment
3056 aNew->SetTotalEnergy((*iter)->Get4Mom
3057 aNew->SetNewlyAdded(true);
3058 aNew->SetCreatorModelID((*iter)->GetC
3059 aNew->SetParentResonanceDef((*iter)->
3060 aNew->SetParentResonanceID((*iter)->G
3061 //G4cout << " Particle Ekin " << aNew
3062 products->push_back(aNew);
3063 }
3064 const G4ParticleDefinition* fragment = 0;
3065 if (currentA == 1 && currentZ == 0) {
3066 fragment = G4Neutron::NeutronDefiniti
3067 } else if (currentA == 1 && currentZ == 1
3068 fragment = G4Proton::ProtonDefinition
3069 } else if (currentA == 2 && currentZ == 1
3070 fragment = G4Deuteron::DeuteronDefini
3071 } else if (currentA == 3 && currentZ == 1
3072 fragment = G4Triton::TritonDefinition
3073 } else if (currentA == 3 && currentZ == 2
3074 fragment = G4He3::He3Definition();
3075 } else if (currentA == 4 && currentZ == 2
3076 fragment = G4Alpha::AlphaDefinition()
3077 } else {
3078 fragment =
3079 G4ParticleTable::GetParticleT
3080 }
3081 if (fragment != 0) {
3082 G4ReactionProduct * theNew = new G4Re
3083 theNew->SetMomentum(G4ThreeVector(0,0
3084 theNew->SetTotalEnergy(massInNucleus)
3085 theNew->SetCreatorModelID(theBIC_ID);
3086 //theNew->SetFormationTime(??0.??);
3087 //G4cout << " Nucleus (" << currentZ
3088 products->push_back(theNew);
3089 }
3090 return products;
3091 }
3092
3093 void G4BinaryCascade::PrintWelcomeMessage()
3094 {
3095 G4cout <<"Thank you for using G4BinaryCas
3096 }
3097
3098 //-------------------------------------------
3099 void G4BinaryCascade::DebugApplyCollisionFail
3100 G4KineticTrackVector * products)
3101 {
3102 G4bool havePion=false;
3103 if (products)
3104 {
3105 for ( auto i =products->cbegin(); i !=
3106 {
3107 G4int PDGcode=std::abs((*i)->GetDefi
3108 if (std::abs(PDGcode)==211 || PDGcod
3109 }
3110 }
3111 if ( !products || havePion)
3112 {
3113 const G4BCAction &action= *collision->G
3114 G4cout << " Collision " << collision <<
3115 << ", with
3116 G4cout << G4endl<<"Initial condition ar
3117 G4cout << "proj: "<<collision->GetPrima
3118 PrintKTVector(collision->GetPrimary());
3119 for(std::size_t it=0; it<collision->Get
3120 {
3121 G4cout << "targ: "
3122 <<collision->GetTargetCollecti
3123 }
3124 PrintKTVector(&collision->GetTargetColl
3125 }
3126 // if ( lateParticleCollision ) G4cout <<
3127 // if ( lateParticleCollision && products
3128 }
3129
3130 //-------------------------------------------
3131
3132 G4bool G4BinaryCascade::CheckChargeAndBaryonN
3133 {
3134 static G4int lastdA(0), lastdZ(0);
3135 G4int iStateA = the3DNucleus->GetMassNumbe
3136 G4int iStateZ = the3DNucleus->GetCharge()
3137
3138 G4int fStateA(0);
3139 G4int fStateZ(0);
3140
3141 G4int CapturedA(0), CapturedZ(0);
3142 G4int secsA(0), secsZ(0);
3143 for (auto i=theCapturedList.cbegin(); i!=t
3144 CapturedA += (*i)->GetDefinition()->Get
3145 CapturedZ += G4lrint((*i)->GetDefinitio
3146 }
3147
3148 for (auto i=theSecondaryList.cbegin(); i!=
3149 if ( (*i)->GetState() != G4KineticTrack
3150 secsA += (*i)->GetDefinition()->GetB
3151 secsZ += G4lrint((*i)->GetDefinition
3152 }
3153 }
3154
3155 for (auto i=theFinalState.cbegin(); i!=the
3156 fStateA += (*i)->GetDefinition()->GetBa
3157 fStateZ += G4lrint((*i)->GetDefinition(
3158 }
3159
3160 G4int deltaA= iStateA - secsA - fStateA -
3161 G4int deltaZ= iStateZ - secsZ - fStateZ -
3162
3163 #ifdef debugCheckChargeAndBaryonNumberverbose
3164 G4cout << where <<" A: iState= "<< iStateA<
3165 G4cout << where <<" Z: iState= "<< iStateZ<
3166 #endif
3167
3168 if (deltaA != 0 || deltaZ!=0 ) {
3169 if (deltaA != lastdA || deltaZ != lastd
3170 G4cout << "baryon/charge imbalance -
3171 << "deltaA " <<deltaA<<", iSta
3172 << ", fStateA "<<fStateA << ",
3173 << "deltaZ "<<deltaZ<<", iStat
3174 << ", fStateZ "<<fStateZ << ",
3175 lastdA=deltaA;
3176 lastdZ=deltaZ;
3177 }
3178 } else { lastdA=lastdZ=0;}
3179
3180 return true;
3181 }
3182 //-------------------------------------------
3183 void G4BinaryCascade::DebugApplyCollision(G4C
3184 G4KineticTrackVector * products)
3185 {
3186
3187 PrintKTVector(collision->GetPrimary(),std
3188 PrintKTVector(&collision->GetTargetCollec
3189 PrintKTVector(products,std::string(" Scat
3190
3191 #ifdef dontUse
3192 G4double thisExcitation(0);
3193 // excitation energy from this collision
3194 // initial state:
3195 G4double initial(0);
3196 G4KineticTrack * kt=collision->GetPrimary
3197 initial += kt->Get4Momentum().e();
3198
3199 G4RKPropagation * RKprop=(G4RKPropagation
3200
3201 initial += RKprop->GetField(kt->GetDefin
3202 initial -= RKprop->GetBarrier(kt->GetDef
3203 G4cout << "prim. E/field/Barr/Sum " << kt
3204 << " " << RKprop->GetFi
3205 << " " << RKprop->GetBa
3206 << " " << initial << G4
3207
3208 G4KineticTrackVector ktv=collision->GetTa
3209 for ( unsigned int it=0; it < ktv.size();
3210 {
3211 kt=ktv[it];
3212 initial += kt->Get4Momentum().e();
3213 thisExcitation += kt->GetDefinition()
3214 - kt->Get4Momentum()
3215 - RKprop->GetField(k
3216 // initial += RKprop->GetField(k
3217 // initial -= RKprop->GetBarrier
3218 G4cout << "Targ. def/E/field/Barr/Sum
3219 << " " << kt->Get4Momentum(
3220 << " " << RKprop->GetField(
3221 << " " << RKprop->GetBarrie
3222 << " " << initial <<" Excit
3223 }
3224
3225 G4double fin(0);
3226 G4double mass_out(0);
3227 G4int product_barions(0);
3228 if ( products )
3229 {
3230 for ( unsigned int it=0; it < product
3231 {
3232 kt=(*products)[it];
3233 fin += kt->Get4Momentum().e();
3234 fin += RKprop->GetField(kt->GetD
3235 fin += RKprop->GetBarrier(kt->Ge
3236 if ( kt->GetDefinition()->GetBary
3237 mass_out += kt->GetDefinition()->
3238 G4cout << "sec. def/E/field/Barr/
3239 << " " << kt->Get4Moment
3240 << " " << RKprop->GetFie
3241 << " " << RKprop->GetBar
3242 << " " << fin << G4endl;
3243 }
3244 }
3245
3246
3247 G4int finalA = currentA;
3248 G4int finalZ = currentZ;
3249 if ( products )
3250 {
3251 finalA -= product_barions;
3252 finalZ -= GetTotalCharge(*products);
3253 }
3254 G4double delta = GetIonMass(currentZ,curr
3255 G4cout << " current/final a,z " << curren
3256 << " delta-mass " << delta<<G4en
3257 fin+=delta;
3258 mass_out = GetIonMass(finalZ,finalA);
3259 G4cout << " initE/ E_out/ Mfinal/ Excit "
3260 << " " << fin << " "
3261 << mass_out<<" "
3262 << currentInitialEnergy - fin -
3263 << G4endl;
3264 currentInitialEnergy-=fin;
3265 #endif
3266 }
3267
3268 //-------------------------------------------
3269 G4bool G4BinaryCascade::DebugFinalEpConservat
3270 G4ReactionProductVector* products)
3271 //-------------------------------------------
3272 {
3273 G4double Efinal(0);
3274 G4ThreeVector pFinal(0);
3275 if (std::abs(theParticleChange.GetWeightC
3276 {
3277 G4cout <<" BIC-weight change " << the
3278 }
3279
3280 for(auto iter = products->cbegin(); iter
3281 {
3282
3283 G4cout << " Secondary E - Ekin / p "
3284 (*iter)->GetDefinition()->GetPa
3285 (*iter)->GetTotalEnergy() << "
3286 (*iter)->GetKineticEnergy()<< "
3287 (*iter)->GetMomentum().x() << "
3288 (*iter)->GetMomentum().y() << "
3289 (*iter)->GetMomentum().z() << G
3290 Efinal += (*iter)->GetTotalEnergy();
3291 pFinal += (*iter)->GetMomentum();
3292 }
3293
3294 G4cout << "e outgoing/ total : " << Efi
3295 G4cout << "BIC E/p delta " <<
3296 (aTrack.Get4Momentum().e()+theIni
3297 " MeV / mom " << (aTrack.Get4Mome
3298
3299 return (aTrack.Get4Momentum().e() + theIn
3300
3301 }
3302 //-------------------------------------------
3303 G4bool G4BinaryCascade::DebugEpConservation(c
3304 //-------------------------------------------
3305 {
3306 G4cout << where << G4endl;
3307 G4LorentzVector psecs, ptgts, pcpts
3308 if (std::abs(theParticleChange.GetWeightC
3309 {
3310 G4cout <<" BIC-weight change " << the
3311 }
3312
3313 std::vector<G4KineticTrack *>::const_iter
3314 for(ktiter = theSecondaryList.cbegin(); k
3315 {
3316
3317 G4cout << " Secondary E - Ekin
3318 (*ktiter)->GetDefinition()->
3319 (*ktiter)->Get4Momentum().e(
3320 (*ktiter)->Get4Momentum().e(
3321 (*ktiter)->Get4Momentum().ve
3322 psecs += (*ktiter)->Get4Momentum()
3323 }
3324
3325 for(ktiter = theTargetList.cbegin(); ktit
3326 {
3327
3328 G4cout << " Target E - Ekin / p
3329 (*ktiter)->GetDefinition()->
3330 (*ktiter)->Get4Momentum().e(
3331 (*ktiter)->Get4Momentum().e(
3332 (*ktiter)->Get4Momentum().ve
3333 ptgts += (*ktiter)->Get4Momentum()
3334 }
3335
3336 for(ktiter = theCapturedList.cbegin(); kt
3337 {
3338
3339 G4cout << " Captured E - Ekin
3340 (*ktiter)->GetDefinition()-
3341 (*ktiter)->Get4Momentum().e
3342 (*ktiter)->Get4Momentum().e
3343 (*ktiter)->Get4Momentum().v
3344 pcpts += (*ktiter)->Get4Momentum(
3345 }
3346
3347 for(ktiter = theFinalState.cbegin(); ktit
3348 {
3349
3350 G4cout << " Finals E - Ekin /
3351 (*ktiter)->GetDefinition()-
3352 (*ktiter)->Get4Momentum().e
3353 (*ktiter)->Get4Momentum().e
3354 (*ktiter)->Get4Momentum().v
3355 pfins += (*ktiter)->Get4Momentum(
3356 }
3357
3358 G4cout << " Secondaries " << psecs << "
3359 <<" Captured " << pcpts << ", Fi
3360 <<" Sum " << psecs + ptgts + pcpts
3361 <<" Sum+PTransfer " << psecs + ptgt
3362 << G4endl<< G4endl;
3363
3364
3365 return true;
3366
3367 }
3368
3369 //-------------------------------------------
3370 G4ReactionProductVector * G4BinaryCascade::Pr
3371 //-------------------------------------------
3372 {
3373 // else
3374 // {
3375 // G4ReactionProduct * aNew=0;
3376 // // return nucleus e and p
3377 // if (fragment != 0 ) {
3378 // aNew = new G4ReactionProduct(G4
3379 // aNew->SetMomentum(fragment->Get
3380 // aNew->SetTotalEnergy(fragment->
3381 // delete fragment;
3382 // fragment=0;
3383 // } else if (products->size() == 0) {
3384 // // FixMe GF: for testing withou
3385 //#include "G4Gamma.hh"
3386 // aNew = new G4ReactionProduct(G4
3387 // aNew->SetMomentum(G4ThreeVector
3388 // aNew->SetTotalEnergy(0.01*MeV);
3389 // }
3390 // if ( aNew != 0 ) products->push_bac
3391 // }
3392 return products;
3393 }
3394
3395 //-------------------------------------------
3396