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