Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc (Version 11.3.0) and /processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc (Version 11.1.1)


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