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 10.3.p3)


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