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