| 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 #include <algorithm> 26 #include <algorithm>
27 #include <vector> 27 #include <vector>
28 #include <cmath> 28 #include <cmath>
29 #include <numeric> 29 #include <numeric>
30 30
31 #include "G4BinaryLightIonReaction.hh" 31 #include "G4BinaryLightIonReaction.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4LorentzVector.hh" 34 #include "G4LorentzVector.hh"
35 #include "G4LorentzRotation.hh" 35 #include "G4LorentzRotation.hh"
36 #include "G4ReactionProductVector.hh" 36 #include "G4ReactionProductVector.hh"
37 #include "G4ping.hh" 37 #include "G4ping.hh"
38 #include "G4Delete.hh" 38 #include "G4Delete.hh"
39 #include "G4Neutron.hh" 39 #include "G4Neutron.hh"
40 #include "G4VNuclearDensity.hh" 40 #include "G4VNuclearDensity.hh"
41 #include "G4FermiMomentum.hh" 41 #include "G4FermiMomentum.hh"
42 #include "G4HadTmpUtil.hh" 42 #include "G4HadTmpUtil.hh"
43 #include "G4PreCompoundModel.hh" 43 #include "G4PreCompoundModel.hh"
44 #include "G4HadronicInteractionRegistry.hh" 44 #include "G4HadronicInteractionRegistry.hh"
45 #include "G4Log.hh" 45 #include "G4Log.hh"
46 #include "G4PhysicsModelCatalog.hh" 46 #include "G4PhysicsModelCatalog.hh"
47 #include "G4HadronicParameters.hh" 47 #include "G4HadronicParameters.hh"
48 48
49 G4int G4BinaryLightIonReaction::theBLIR_ID = - 49 G4int G4BinaryLightIonReaction::theBLIR_ID = -1;
50 50
51 //#define debug_G4BinaryLightIonReaction 51 //#define debug_G4BinaryLightIonReaction
52 //#define debug_BLIR_finalstate 52 //#define debug_BLIR_finalstate
53 //#define debug_BLIR_result 53 //#define debug_BLIR_result
54 54
55 G4BinaryLightIonReaction::G4BinaryLightIonReac 55 G4BinaryLightIonReaction::G4BinaryLightIonReaction(G4VPreCompoundModel* ptr)
56 : G4HadronicInteraction("Binary Light Ion Casc 56 : G4HadronicInteraction("Binary Light Ion Cascade"),
57 theProjectileFragmentation(ptr), 57 theProjectileFragmentation(ptr),
58 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spect 58 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
59 projectile3dNucleus(0),target3dNucleus(0) 59 projectile3dNucleus(0),target3dNucleus(0)
60 { 60 {
61 if(!ptr) { 61 if(!ptr) {
62 G4HadronicInteraction* p = 62 G4HadronicInteraction* p =
63 G4HadronicInteractionRegistry::Instance( 63 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
64 G4VPreCompoundModel* pre = static_cast<G4V 64 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
65 if(!pre) { pre = new G4PreCompoundModel(); 65 if(!pre) { pre = new G4PreCompoundModel(); }
66 theProjectileFragmentation = pre; 66 theProjectileFragmentation = pre;
67 } 67 }
68 theModel = new G4BinaryCascade(theProjectile 68 theModel = new G4BinaryCascade(theProjectileFragmentation);
69 theHandler = theProjectileFragmentation->Get 69 theHandler = theProjectileFragmentation->GetExcitationHandler();
70 theBLIR_ID = G4PhysicsModelCatalog::GetM 70 theBLIR_ID = G4PhysicsModelCatalog::GetModelID("model_G4BinaryLightIonReaction");
71 debug_G4BinaryLightIonReactionResults = G4Ha 71 debug_G4BinaryLightIonReactionResults = G4HadronicParameters::Instance()->GetBinaryDebug();
72 } 72 }
73 73
74 G4BinaryLightIonReaction::~G4BinaryLightIonRea 74 G4BinaryLightIonReaction::~G4BinaryLightIonReaction()
75 {} 75 {}
76 76
77 void G4BinaryLightIonReaction::ModelDescriptio 77 void G4BinaryLightIonReaction::ModelDescription(std::ostream& outFile) const
78 { 78 {
79 outFile << "G4Binary Light Ion Cascade is an 79 outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
80 << "using G4BinaryCasacde to model the i 80 << "using G4BinaryCasacde to model the interaction of a light\n"
81 << "nucleus with a nucleus.\n" 81 << "nucleus with a nucleus.\n"
82 << "The lighter of the two nuclei is tre 82 << "The lighter of the two nuclei is treated like a set of projectiles\n"
83 << "which are transported simultaneously 83 << "which are transported simultaneously through the heavier nucleus.\n";
84 } 84 }
85 85
86 //-------------------------------------------- 86 //--------------------------------------------------------------------------------
87 struct ReactionProduct4Mom 87 struct ReactionProduct4Mom
88 { 88 {
89 G4LorentzVector operator()(G4LorentzVector 89 G4LorentzVector operator()(G4LorentzVector a,G4ReactionProduct* b) {return a + G4LorentzVector(b->GetMomentum(), b->GetTotalEnergy() );}
90 }; 90 };
91 91
92 G4HadFinalState *G4BinaryLightIonReaction:: 92 G4HadFinalState *G4BinaryLightIonReaction::
93 ApplyYourself(const G4HadProjectile &aTrack, G 93 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
94 { 94 {
95 if(debug_G4BinaryLightIonReactionResults) G4 95 if(debug_G4BinaryLightIonReactionResults) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl;
96 G4ping debug("debug_G4BinaryLightIonReaction 96 G4ping debug("debug_G4BinaryLightIonReaction");
97 pA=aTrack.GetDefinition()->GetBaryonNumber() 97 pA=aTrack.GetDefinition()->GetBaryonNumber();
98 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCha 98 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus);
99 tA=targetNucleus.GetA_asInt(); 99 tA=targetNucleus.GetA_asInt();
100 tZ=targetNucleus.GetZ_asInt(); 100 tZ=targetNucleus.GetZ_asInt();
101 G4double timePrimary = aTrack.GetGlobalTime( 101 G4double timePrimary = aTrack.GetGlobalTime();
102 G4LorentzVector mom(aTrack.Get4Momentum()); 102 G4LorentzVector mom(aTrack.Get4Momentum());
103 //G4cout << "proj mom : " << mom << G4endl; 103 //G4cout << "proj mom : " << mom << G4endl;
104 G4LorentzRotation toBreit(mom.boostVector()) 104 G4LorentzRotation toBreit(mom.boostVector());
105 105
106 G4bool swapped=SetLighterAsProjectile(mom, t 106 G4bool swapped=SetLighterAsProjectile(mom, toBreit);
107 //G4cout << "after swap, swapped? / mom " < 107 //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
108 G4ReactionProductVector * result = 0; 108 G4ReactionProductVector * result = 0;
109 G4ReactionProductVector * cascaders=0; //new 109 G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
110 // G4double m_nucl(0); // to check energ 110 // G4double m_nucl(0); // to check energy balance
111 111
112 112
113 // G4double m1=G4ParticleTable::GetPartic 113 // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
114 // G4cout << "Entering the decision point 114 // G4cout << "Entering the decision point "
115 // << (mom.t()-mom.mag())/pA << " 115 // << (mom.t()-mom.mag())/pA << " "
116 // << pA<<" "<< pZ<<" " 116 // << pA<<" "<< pZ<<" "
117 // << tA<<" "<< tZ<<G4endl 117 // << tA<<" "<< tZ<<G4endl
118 // << " "<<mom.t()-mom.mag()<<" " 118 // << " "<<mom.t()-mom.mag()<<" "
119 // << mom.t()- m1<<G4endl; 119 // << mom.t()- m1<<G4endl;
120 if( (mom.t()-mom.mag())/pA < 50*MeV ) 120 if( (mom.t()-mom.mag())/pA < 50*MeV )
121 { 121 {
122 // G4cout << "Using pre-compound only 122 // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
123 // m_nucl = mom.mag(); 123 // m_nucl = mom.mag();
124 cascaders=FuseNucleiAndPrompound(mom); 124 cascaders=FuseNucleiAndPrompound(mom);
125 if( !cascaders ) 125 if( !cascaders )
126 { 126 {
127 127
128 // abort!! happens for too low e 128 // abort!! happens for too low energy for nuclei to fuse
129 129
130 theResult.Clear(); 130 theResult.Clear();
131 theResult.SetStatusChange(isAliv 131 theResult.SetStatusChange(isAlive);
132 theResult.SetEnergyChange(aTrack 132 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
133 theResult.SetMomentumChange(aTra 133 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
134 return &theResult; 134 return &theResult;
135 } 135 }
136 } 136 }
137 else 137 else
138 { 138 {
139 result=Interact(mom,toBreit); 139 result=Interact(mom,toBreit);
140 140
141 if(! result ) 141 if(! result )
142 { 142 {
143 // abort!! 143 // abort!!
144 144
145 G4cerr << "G4BinaryLightIonReaction 145 G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
146 G4cerr << " Primary " << aTrack.Get 146 G4cerr << " Primary " << aTrack.GetDefinition()
147 << ", (A,Z)=(" << aTrack.GetDefi 147 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
148 << "," << aTrack.GetDefinition() 148 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
149 << ", kinetic energy " << aTrack 149 << ", kinetic energy " << aTrack.GetKineticEnergy()
150 << G4endl; 150 << G4endl;
151 G4cerr << " Target nucleus (A,Z)=(" 151 G4cerr << " Target nucleus (A,Z)=("
152 << (swapped?pA:tA) << "," 152 << (swapped?pA:tA) << ","
153 << (swapped?pZ:tZ) << ")" << 153 << (swapped?pZ:tZ) << ")" << G4endl;
154 G4cerr << " if frequent, please sub 154 G4cerr << " if frequent, please submit above information as bug report"
155 << G4endl << G4endl; 155 << G4endl << G4endl;
156 156
157 theResult.Clear(); 157 theResult.Clear();
158 theResult.SetStatusChange(isAlive); 158 theResult.SetStatusChange(isAlive);
159 theResult.SetEnergyChange(aTrack.Ge 159 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
160 theResult.SetMomentumChange(aTrack. 160 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
161 return &theResult; 161 return &theResult;
162 } 162 }
163 163
164 // Calculate excitation energy, 164 // Calculate excitation energy,
165 G4double theStatisticalExEnergy = GetProj 165 G4double theStatisticalExEnergy = GetProjectileExcitation();
166 166
167 167
168 pInitialState = mom; 168 pInitialState = mom;
169 //G4cout << "BLIC: pInitialState from 169 //G4cout << "BLIC: pInitialState from aTrack : " << pInitialState;
170 pInitialState.setT(pInitialState.getT() + 170 pInitialState.setT(pInitialState.getT() +
171 G4ParticleTable::GetParticleTable()->GetIo 171 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA));
172 //G4cout << "BLIC: target nucleus added 172 //G4cout << "BLIC: target nucleus added : " << pInitialState << G4endl;
173 173
174 delete target3dNucleus;target3dNucleus=0; 174 delete target3dNucleus;target3dNucleus=0;
175 delete projectile3dNucleus;projectile3dNu 175 delete projectile3dNucleus;projectile3dNucleus=0;
176 176
177 G4ReactionProductVector * spectators= new 177 G4ReactionProductVector * spectators= new G4ReactionProductVector;
178 178
179 cascaders = new G4ReactionProductVector; 179 cascaders = new G4ReactionProductVector;
180 180
181 G4LorentzVector pspectators=SortResult(re 181 G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
182 // this also sets spectatorA and 182 // this also sets spectatorA and spectatorZ
183 183
184 // pFinalState=std::accumulate(casca 184 // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
185 185
186 std::vector<G4ReactionProduct *>::iterato 186 std::vector<G4ReactionProduct *>::iterator iter;
187 187
188 // G4cout << "pInitialState, pFin 188 // G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
189 // if ( spectA-spectatorA !=0 || spec 189 // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
190 // { 190 // {
191 // G4cout << "spect Nucl != spect 191 // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
192 // spectatorA <<" "<< spectatorZ << 192 // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
193 // } 193 // }
194 delete result; 194 delete result;
195 result=0; 195 result=0;
196 G4LorentzVector momentum(pInitialState-pF 196 G4LorentzVector momentum(pInitialState-pFinalState);
197 G4int loopcount(0); 197 G4int loopcount(0);
198 //G4cout << "BLIC: momentum, pspectato 198 //G4cout << "BLIC: momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
199 while (std::abs(momentum.e()-pspectators. 199 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) /* Loop checking, 31.08.2015, G.Folger */
200 200 // see if on loopcount
201 { 201 {
202 G4LorentzVector pCorrect(pInitialState- 202 G4LorentzVector pCorrect(pInitialState-pspectators);
203 //G4cout << "BLIC:: BIC nonconservatio 203 //G4cout << "BLIC:: BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
204 // Correct outgoing casacde particles.. 204 // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
205 G4bool EnergyIsCorrect=EnergyAndMomentu 205 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
206 if ( ! EnergyIsCorrect && debug_G4Binar 206 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
207 { 207 {
208 G4cout << "Warning - G4BinaryLightIon 208 G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
209 } 209 }
210 pFinalState=G4LorentzVector(0,0,0,0); 210 pFinalState=G4LorentzVector(0,0,0,0);
211 for(iter=cascaders->begin(); iter!=casc 211 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
212 { 212 {
213 pFinalState += G4LorentzVector( (*ite 213 pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() );
214 } 214 }
215 momentum=pInitialState-pFinalState; 215 momentum=pInitialState-pFinalState;
216 if (++loopcount > 10 ) 216 if (++loopcount > 10 )
217 { 217 {
218 break; 218 break;
219 } 219 }
220 } 220 }
221 221
222 // Check if Energy/Momemtum is now ok, if 222 // Check if Energy/Momemtum is now ok, if not return initial state
223 if ( std::abs(momentum.e()-pspectators.e( 223 if ( std::abs(momentum.e()-pspectators.e()) > 10*MeV )
224 { 224 {
225 for (iter=spectators->begin();iter!=specta 225 for (iter=spectators->begin();iter!=spectators->end();iter++)
226 { 226 {
227 delete *iter; 227 delete *iter;
228 } 228 }
229 delete spectators; 229 delete spectators;
230 for(iter=cascaders->begin(); iter!=cascad 230 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
231 { 231 {
232 delete *iter; 232 delete *iter;
233 } 233 }
234 delete cascaders; 234 delete cascaders;
235 235
236 G4cout << "G4BinaryLightIonReaction.cc: m 236 G4cout << "G4BinaryLightIonReaction.cc: mom check: " << G4endl
237 << " initial - final " << momentum 237 << " initial - final " << momentum << " 3.mag "<< momentum.vect().mag() << G4endl
238 << " .. pInitialState/pFinalState/s 238 << " .. pInitialState/pFinalState/spectators " << G4endl
239 << pInitialState << G4endl 239 << pInitialState << G4endl
240 << pFinalState << G4endl 240 << pFinalState << G4endl
241 << pspectators << G4endl 241 << pspectators << G4endl
242 << " .. A,Z " << spectatorA <<" "<< 242 << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
243 G4cout << "G4BinaryLightIonReaction inval 243 G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
244 G4cout << " Primary " << aTrack.GetDefini 244 G4cout << " Primary " << aTrack.GetDefinition()
245 << ", (A,Z)=(" << aTrack.GetDefin 245 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
246 << "," << aTrack.GetDefinition()- 246 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
247 << ", kinetic energy " << aTrack. 247 << ", kinetic energy " << aTrack.GetKineticEnergy()
248 << G4endl; 248 << G4endl;
249 G4cout << " Target nucleus (A,Z)=(" << t 249 G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt()
250 << "," << targetNucleus.GetZ_ 250 << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
251 G4cout << " if frequent, please submit ab 251 G4cout << " if frequent, please submit above information as bug report"
252 << G4endl << G4endl; 252 << G4endl << G4endl;
253 #ifdef debug_G4BinaryLightIonReaction 253 #ifdef debug_G4BinaryLightIonReaction
254 G4ExceptionDescription ed; 254 G4ExceptionDescription ed;
255 ed << "G4BinaryLightIonreaction: Ter 255 ed << "G4BinaryLightIonreaction: Terminate for above error" << G4endl;
256 G4Exception("G4BinaryLightIonreactio 256 G4Exception("G4BinaryLightIonreaction::ApplyYourSelf()", "BLIC001", FatalException,
257 ed); 257 ed);
258 258
259 #endif 259 #endif
260 theResult.Clear(); 260 theResult.Clear();
261 theResult.SetStatusChange(isAlive); 261 theResult.SetStatusChange(isAlive);
262 theResult.SetEnergyChange(aTrack.GetKinet 262 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
263 theResult.SetMomentumChange(aTrack.Get4Mo 263 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
264 return &theResult; 264 return &theResult;
265 265
266 } 266 }
267 if (spectatorA > 0 ) 267 if (spectatorA > 0 )
268 { 268 {
269 // DeExciteSpectatorNucleus() also 269 // DeExciteSpectatorNucleus() also handles also case of A=1, Z=0,1
270 DeExciteSpectatorNucleus(specta 270 DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
271 } else { // no spectators 271 } else { // no spectators
272 delete spectators; 272 delete spectators;
273 } 273 }
274 } 274 }
275 // Rotate to lab 275 // Rotate to lab
276 G4LorentzRotation toZ; 276 G4LorentzRotation toZ;
277 toZ.rotateZ(-1*mom.phi()); 277 toZ.rotateZ(-1*mom.phi());
278 toZ.rotateY(-1*mom.theta()); 278 toZ.rotateY(-1*mom.theta());
279 G4LorentzRotation toLab(toZ.inverse()); 279 G4LorentzRotation toLab(toZ.inverse());
280 280
281 // Fill the particle change, while rotating. 281 // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
282 // theResult.Clear(); 282 // theResult.Clear();
283 theResult.Clear(); 283 theResult.Clear();
284 theResult.SetStatusChange(stopAndKill); 284 theResult.SetStatusChange(stopAndKill);
285 G4LorentzVector ptot(0); 285 G4LorentzVector ptot(0);
286 #ifdef debug_BLIR_result 286 #ifdef debug_BLIR_result
287 G4LorentzVector p_raw; 287 G4LorentzVector p_raw;
288 #endif 288 #endif
289 //G4int i=0; 289 //G4int i=0;
290 290
291 G4ReactionProductVector::iterator iter 291 G4ReactionProductVector::iterator iter;
292 for(iter=cascaders->begin(); iter!=cascaders 292 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
293 { 293 {
294 if((*iter)->GetNewlyAdded()) 294 if((*iter)->GetNewlyAdded())
295 { 295 {
296 G4DynamicParticle * aNewDP = 296 G4DynamicParticle * aNewDP =
297 new G4DynamicParticle((*iter)->GetDe 297 new G4DynamicParticle((*iter)->GetDefinition(),
298 (*iter)->GetTotalEnergy(), 298 (*iter)->GetTotalEnergy(),
299 (*iter)->GetMomentum() ); 299 (*iter)->GetMomentum() );
300 G4LorentzVector tmp = aNewDP->Get4Moment 300 G4LorentzVector tmp = aNewDP->Get4Momentum();
301 #ifdef debug_BLIR_result 301 #ifdef debug_BLIR_result
302 p_raw+= tmp; 302 p_raw+= tmp;
303 #endif 303 #endif
304 if(swapped) 304 if(swapped)
305 { 305 {
306 tmp*=toBreit.inverse(); 306 tmp*=toBreit.inverse();
307 tmp.setVect(-tmp.vect()); 307 tmp.setVect(-tmp.vect());
308 } 308 }
309 tmp *= toLab; 309 tmp *= toLab;
310 aNewDP->Set4Momentum(tmp); 310 aNewDP->Set4Momentum(tmp);
311 G4HadSecondary aNew = G4HadSecondary(aNe 311 G4HadSecondary aNew = G4HadSecondary(aNewDP);
312 G4double time = 0; 312 G4double time = 0; //(*iter)->GetCreationTime();
313 //if(time < 0.0) { time = 0.0; } 313 //if(time < 0.0) { time = 0.0; }
314 aNew.SetTime(timePrimary + time); 314 aNew.SetTime(timePrimary + time);
315 //aNew.SetCreatorModelID((*iter)-> 315 //aNew.SetCreatorModelID((*iter)->GetCreatorModelID()); //AR-02Aug2021 : For some reasons, it does NOT work!
316 aNew.SetCreatorModelID(theBLIR_ID) 316 aNew.SetCreatorModelID(theBLIR_ID);
317 317
318 theResult.AddSecondary(aNew); 318 theResult.AddSecondary(aNew);
319 ptot += tmp; 319 ptot += tmp;
320 //G4cout << "BLIC: Secondary " < 320 //G4cout << "BLIC: Secondary " << aNew->GetDefinition()->GetParticleName()
321 // <<" "<< aNew->GetMomen 321 // <<" "<< aNew->GetMomentum()<<" "<< aNew->GetTotalEnergy() << G4endl;
322 } 322 }
323 delete *iter; 323 delete *iter;
324 } 324 }
325 delete cascaders; 325 delete cascaders;
326 326
327 #ifdef debug_BLIR_result 327 #ifdef debug_BLIR_result
328 //G4cout << "Result analysis, secondaries " 328 //G4cout << "Result analysis, secondaries " << theResult.GetNumberOfSecondaries() << G4endl;
329 //G4cout << "p_tot_raw " << p_raw << " sum p 329 //G4cout << "p_tot_raw " << p_raw << " sum p final " << ptot << G4endl;
330 G4double m_nucl= G4ParticleTable::GetPartic 330 G4double m_nucl= G4ParticleTable::GetParticleTable()->GetIonTable()->
331 GetIonMass(targetNucleus.GetZ_asInt( 331 GetIonMass(targetNucleus.GetZ_asInt(),targetNucleus.GetA_asInt());
332 // delete? tZ=targetNucleus.GetZ_asInt(); 332 // delete? tZ=targetNucleus.GetZ_asInt();
333 333
334 //G4cout << "BLIC Energy conservation initia 334 //G4cout << "BLIC Energy conservation initial/primary/nucleus/final/delta(init-final) "
335 // << aTrack.GetTotalEnergy() + m_nuc 335 // << aTrack.GetTotalEnergy() + m_nucl <<" "<< aTrack.GetTotalEnergy() <<" "<< m_nucl <<" "<<ptot.e()
336 // <<" "<< aTrack.GetTotalEnergy() + m_ 336 // <<" "<< aTrack.GetTotalEnergy() + m_nucl - ptot.e() << G4endl;
337 G4cout << "BLIC momentum conservation " << a 337 G4cout << "BLIC momentum conservation " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl)
338 << " ptot " << ptot << " delta " << aTra 338 << " ptot " << ptot << " delta " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot
339 << " 3mom.mag() " << (aTrack.Get4 339 << " 3mom.mag() " << (aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot).vect().mag() << G4endl;
340 #endif 340 #endif
341 341
342 if(debug_G4BinaryLightIonReactionResults) G4 342 if(debug_G4BinaryLightIonReactionResults) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl;
343 343
344 return &theResult; 344 return &theResult;
345 } 345 }
346 346
347 //-------------------------------------------- 347 //--------------------------------------------------------------------------------
348 348
349 //******************************************** 349 //****************************************************************************
350 G4bool G4BinaryLightIonReaction::EnergyAndMome 350 G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector(
351 G4ReactionProductVector* Output, G4Lorentz 351 G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom)
352 //******************************************** 352 //****************************************************************************
353 { 353 {
354 const int nAttemptScale = 2500; 354 const int nAttemptScale = 2500;
355 const double ErrLimit = 1.E-6; 355 const double ErrLimit = 1.E-6;
356 if (Output->empty()) 356 if (Output->empty())
357 return TRUE; 357 return TRUE;
358 G4LorentzVector SumMom(0,0,0,0); 358 G4LorentzVector SumMom(0,0,0,0);
359 G4double SumMass = 0; 359 G4double SumMass = 0;
360 G4double TotalCollisionMass = TotalCo 360 G4double TotalCollisionMass = TotalCollisionMom.m();
361 size_t i = 0; 361 size_t i = 0;
362 // Calculate sum hadron 4-momenta and summin 362 // Calculate sum hadron 4-momenta and summing hadron mass
363 for(i = 0; i < Output->size(); i++) 363 for(i = 0; i < Output->size(); i++)
364 { 364 {
365 SumMom += G4LorentzVector((*Output)[i]->G 365 SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
366 SumMass += (*Output)[i]->GetDefinition()-> 366 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
367 } 367 }
368 // G4cout << " E/P corrector, SumMass, Sum 368 // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass "
369 // << SumMass <<" "<< SumMom.m2() << 369 // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl;
370 if (SumMass > TotalCollisionMass) return FAL 370 if (SumMass > TotalCollisionMass) return FALSE;
371 SumMass = SumMom.m2(); 371 SumMass = SumMom.m2();
372 if (SumMass < 0) return FALSE; 372 if (SumMass < 0) return FALSE;
373 SumMass = std::sqrt(SumMass); 373 SumMass = std::sqrt(SumMass);
374 374
375 // Compute c.m.s. hadron velocity and boost 375 // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
376 G4ThreeVector Beta = -SumMom.boostVector(); 376 G4ThreeVector Beta = -SumMom.boostVector();
377 //G4cout << " == pre boost 2 "<< SumMo 377 //G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl;
378 //--old Output->Boost(Beta); 378 //--old Output->Boost(Beta);
379 for(i = 0; i < Output->size(); i++) 379 for(i = 0; i < Output->size(); i++)
380 { 380 {
381 G4LorentzVector mom = G4LorentzVector((*Ou 381 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
382 mom *= Beta; 382 mom *= Beta;
383 (*Output)[i]->SetMomentum(mom.vect()); 383 (*Output)[i]->SetMomentum(mom.vect());
384 (*Output)[i]->SetTotalEnergy(mom.e()); 384 (*Output)[i]->SetTotalEnergy(mom.e());
385 } 385 }
386 386
387 // Scale total c.m.s. hadron energy (hadron 387 // Scale total c.m.s. hadron energy (hadron system mass).
388 // It should be equal interaction mass 388 // It should be equal interaction mass
389 G4double Scale = 0,OldScale=0; 389 G4double Scale = 0,OldScale=0;
390 G4double factor = 1.; 390 G4double factor = 1.;
391 G4int cAttempt = 0; 391 G4int cAttempt = 0;
392 G4double Sum = 0; 392 G4double Sum = 0;
393 G4bool success = false; 393 G4bool success = false;
394 for(cAttempt = 0; cAttempt < nAttemptScale; 394 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
395 { 395 {
396 Sum = 0; 396 Sum = 0;
397 for(i = 0; i < Output->size(); i++) 397 for(i = 0; i < Output->size(); i++)
398 { 398 {
399 G4LorentzVector HadronMom = G4LorentzVec 399 G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
400 HadronMom.setVect(HadronMom.vect()+ fact 400 HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
401 G4double E = std::sqrt(HadronMom.vect(). 401 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
402 HadronMom.setE(E); 402 HadronMom.setE(E);
403 (*Output)[i]->SetMomentum(HadronMom.vect 403 (*Output)[i]->SetMomentum(HadronMom.vect());
404 (*Output)[i]->SetTotalEnergy(HadronMom.e 404 (*Output)[i]->SetTotalEnergy(HadronMom.e());
405 Sum += E; 405 Sum += E;
406 } 406 }
407 OldScale=Scale; 407 OldScale=Scale;
408 Scale = TotalCollisionMass/Sum - 1; 408 Scale = TotalCollisionMass/Sum - 1;
409 // G4cout << "E/P corr - " << cAttempt << 409 // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl;
410 if (std::abs(Scale) <= ErrLimit 410 if (std::abs(Scale) <= ErrLimit
411 || OldScale == Scale) // protect ' 411 || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below
412 { 412 {
413 if (debug_G4BinaryLightIonReactionResult 413 if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl;
414 success = true; 414 success = true;
415 break; 415 break;
416 } 416 }
417 if ( cAttempt > 10 ) 417 if ( cAttempt > 10 )
418 { 418 {
419 // G4cout << " speed it up? " << 419 // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl;
420 factor=std::max(1.,G4Log(std::abs(OldSca 420 factor=std::max(1.,G4Log(std::abs(OldScale/(OldScale-Scale))));
421 // G4cout << " ? factor ? " << factor 421 // G4cout << " ? factor ? " << factor << G4endl;
422 } 422 }
423 } 423 }
424 424
425 if( (!success) && debug_G4BinaryLightIonRea 425 if( (!success) && debug_G4BinaryLightIonReactionResults)
426 { 426 {
427 G4cout << "G4G4BinaryLightIonReaction::Ene 427 G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl;
428 G4cout << " Scale not unity at end of it 428 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
429 G4cout << " Increase number of attempts 429 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
430 } 430 }
431 431
432 // Compute c.m.s. interaction velocity and K 432 // Compute c.m.s. interaction velocity and KTV back boost
433 Beta = TotalCollisionMom.boostVector(); 433 Beta = TotalCollisionMom.boostVector();
434 //--old Output->Boost(Beta); 434 //--old Output->Boost(Beta);
435 for(i = 0; i < Output->size(); i++) 435 for(i = 0; i < Output->size(); i++)
436 { 436 {
437 G4LorentzVector mom = G4LorentzVector((*Ou 437 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
438 mom *= Beta; 438 mom *= Beta;
439 (*Output)[i]->SetMomentum(mom.vect()); 439 (*Output)[i]->SetMomentum(mom.vect());
440 (*Output)[i]->SetTotalEnergy(mom.e()); 440 (*Output)[i]->SetTotalEnergy(mom.e());
441 } 441 }
442 return TRUE; 442 return TRUE;
443 } 443 }
444 G4bool G4BinaryLightIonReaction::SetLighterAsP 444 G4bool G4BinaryLightIonReaction::SetLighterAsProjectile(G4LorentzVector & mom,const G4LorentzRotation & toBreit)
445 { 445 {
446 G4bool swapped = false; 446 G4bool swapped = false;
447 if(tA<pA) 447 if(tA<pA)
448 { 448 {
449 swapped = true; 449 swapped = true;
450 G4int tmp(0); 450 G4int tmp(0);
451 tmp = tA; tA=pA; pA=tmp; 451 tmp = tA; tA=pA; pA=tmp;
452 tmp = tZ; tZ=pZ; pZ=tmp; 452 tmp = tZ; tZ=pZ; pZ=tmp;
453 G4double m1=G4ParticleTable::GetParticle 453 G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
454 G4LorentzVector it(m1, G4ThreeVector(0,0 454 G4LorentzVector it(m1, G4ThreeVector(0,0,0));
455 mom = toBreit*it; 455 mom = toBreit*it;
456 } 456 }
457 return swapped; 457 return swapped;
458 } 458 }
459 G4ReactionProductVector * G4BinaryLightIonReac 459 G4ReactionProductVector * G4BinaryLightIonReaction::FuseNucleiAndPrompound(const G4LorentzVector & mom)
460 { 460 {
461 // Check if kinematically nuclei can fuse. 461 // Check if kinematically nuclei can fuse.
462 G4double mFused=G4ParticleTable::GetParticl 462 G4double mFused=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ+tZ,pA+tA);
463 G4double mTarget=G4ParticleTable::GetPartic 463 G4double mTarget=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA);
464 G4LorentzVector pCompound(mom.e()+mTarget,m 464 G4LorentzVector pCompound(mom.e()+mTarget,mom.vect());
465 G4double m2Compound=pCompound.m2(); 465 G4double m2Compound=pCompound.m2();
466 if (m2Compound < sqr(mFused) ) { 466 if (m2Compound < sqr(mFused) ) {
467 //G4cout << "G4BLIC: projectile p, mTarge 467 //G4cout << "G4BLIC: projectile p, mTarget, mFused, mCompound, delta: " <<mom << " " << mTarget << " " << mFused
468 // << " " << sqrt(m2Compound)<< " " 468 // << " " << sqrt(m2Compound)<< " " << sqrt(m2Compound) - mFused << G4endl;
469 return 0; 469 return 0;
470 } 470 }
471 471
472 G4Fragment aPreFrag; 472 G4Fragment aPreFrag;
473 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA); 473 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA);
474 aPreFrag.SetNumberOfParticles(pA); 474 aPreFrag.SetNumberOfParticles(pA);
475 aPreFrag.SetNumberOfCharged(pZ); 475 aPreFrag.SetNumberOfCharged(pZ);
476 aPreFrag.SetNumberOfHoles(0); 476 aPreFrag.SetNumberOfHoles(0);
477 //GF FIXME: whyusing plop in z direction? t 477 //GF FIXME: whyusing plop in z direction? this will not conserve momentum?
478 //G4ThreeVector plop(0.,0., mom.vect().mag( 478 //G4ThreeVector plop(0.,0., mom.vect().mag());
479 //G4LorentzVector aL(mom.t()+mTarget, plop) 479 //G4LorentzVector aL(mom.t()+mTarget, plop);
480 G4LorentzVector aL(mom.t()+mTarget,mom.vect 480 G4LorentzVector aL(mom.t()+mTarget,mom.vect());
481 aPreFrag.SetMomentum(aL); 481 aPreFrag.SetMomentum(aL);
482 482
483 483
484 //G4cout << "Fragment INFO "<< pA+tA 484 //G4cout << "Fragment INFO "<< pA+tA <<" "<<pZ+tZ<<" "
485 // << aL <<" "<<G4endl << aPreF 485 // << aL <<" "<<G4endl << aPreFrag << G4endl;
486 G4ReactionProductVector * cascaders = thePr 486 G4ReactionProductVector * cascaders = theProjectileFragmentation->DeExcite(aPreFrag);
487 //G4double tSum = 0; 487 //G4double tSum = 0;
488 for(size_t count = 0; count<cascaders->size 488 for(size_t count = 0; count<cascaders->size(); count++)
489 { 489 {
490 cascaders->operator[](count)->SetNewlyAd 490 cascaders->operator[](count)->SetNewlyAdded(true);
491 //tSum += cascaders->operator[](count)-> 491 //tSum += cascaders->operator[](count)->GetKineticEnergy();
492 } 492 }
493 // G4cout << "Exiting pre-compound on 493 // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl;
494 return cascaders; 494 return cascaders;
495 } 495 }
496 G4ReactionProductVector * G4BinaryLightIonReac 496 G4ReactionProductVector * G4BinaryLightIonReaction::Interact(G4LorentzVector & mom, const G4LorentzRotation & toBreit)
497 { 497 {
498 G4ReactionProductVector * result = 0; 498 G4ReactionProductVector * result = 0;
499 G4double projectileMass(0); 499 G4double projectileMass(0);
500 G4LorentzVector it; 500 G4LorentzVector it;
501 501
502 G4int tryCount(0); 502 G4int tryCount(0);
503 do 503 do
504 { 504 {
505 ++tryCount; 505 ++tryCount;
506 projectile3dNucleus = new G4Fancy3DNu 506 projectile3dNucleus = new G4Fancy3DNucleus;
507 projectile3dNucleus->Init(pA, pZ); 507 projectile3dNucleus->Init(pA, pZ);
508 projectile3dNucleus->CenterNucleons() 508 projectile3dNucleus->CenterNucleons();
509 projectileMass=G4ParticleTable::GetPa 509 projectileMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
510 projectile3dNucleus->GetCharge( 510 projectile3dNucleus->GetCharge(),projectile3dNucleus->GetMassNumber());
511 it=toBreit * G4LorentzVector(projecti 511 it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0));
512 512
513 target3dNucleus = new G4Fancy3DNucleu 513 target3dNucleus = new G4Fancy3DNucleus;
514 target3dNucleus->Init(tA, tZ); 514 target3dNucleus->Init(tA, tZ);
515 G4double impactMax = target3dNucleus- 515 G4double impactMax = target3dNucleus->GetOuterRadius()+projectile3dNucleus->GetOuterRadius();
516 // G4cout << "out radius - nuc 516 // G4cout << "out radius - nucleus - projectile " << target3dNucleus->GetOuterRadius()/fermi << " - " << projectile3dNucleus->GetOuterRadius()/fermi << G4endl;
517 G4double aX=(2.*G4UniformRand()-1.)*i 517 G4double aX=(2.*G4UniformRand()-1.)*impactMax;
518 G4double aY=(2.*G4UniformRand()-1.)*i 518 G4double aY=(2.*G4UniformRand()-1.)*impactMax;
519 G4ThreeVector pos(aX, aY, -2.*impactM 519 G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi);
520 520
521 G4KineticTrackVector * initalState = 521 G4KineticTrackVector * initalState = new G4KineticTrackVector;
522 projectile3dNucleus->StartLoop(); 522 projectile3dNucleus->StartLoop();
523 G4Nucleon * aNuc; 523 G4Nucleon * aNuc;
524 G4LorentzVector tmpV(0,0,0,0); 524 G4LorentzVector tmpV(0,0,0,0);
525 #ifdef debug_BLIR_finalstate 525 #ifdef debug_BLIR_finalstate
526 G4LorentzVector pinitial; 526 G4LorentzVector pinitial;
527 #endif 527 #endif
528 G4LorentzVector nucleonMom(1./pA*mom) 528 G4LorentzVector nucleonMom(1./pA*mom);
529 nucleonMom.setZ(nucleonMom.vect().mag 529 nucleonMom.setZ(nucleonMom.vect().mag());
530 nucleonMom.setX(0); 530 nucleonMom.setX(0);
531 nucleonMom.setY(0); 531 nucleonMom.setY(0);
532 theFermi.Init(pA,pZ); 532 theFermi.Init(pA,pZ);
533 while( (aNuc=projectile3dNucleus->Get 533 while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */
534 { 534 {
535 G4LorentzVector p4 = aNuc->GetMome 535 G4LorentzVector p4 = aNuc->GetMomentum();
536 tmpV+=p4; 536 tmpV+=p4;
537 G4ThreeVector nucleonPosition(aNuc 537 G4ThreeVector nucleonPosition(aNuc->GetPosition());
538 G4double density=(projectile3dNucl 538 G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition);
539 nucleonPosition += pos; 539 nucleonPosition += pos;
540 G4KineticTrack * it1 = new G4Kinet 540 G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom );
541 it1->SetState(G4KineticTrack::outs 541 it1->SetState(G4KineticTrack::outside);
542 G4double pfermi= theFermi.GetFermi 542 G4double pfermi= theFermi.GetFermiMomentum(density);
543 G4double mass = aNuc->GetDefinitio 543 G4double mass = aNuc->GetDefinition()->GetPDGMass();
544 G4double Efermi= std::sqrt( sqr(ma 544 G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass;
545 it1->SetProjectilePotential(-Eferm 545 it1->SetProjectilePotential(-Efermi);
546 initalState->push_back(it1); 546 initalState->push_back(it1);
547 #ifdef debug_BLIR_finalstate 547 #ifdef debug_BLIR_finalstate
548 pinitial += it1->Get4Momentum() 548 pinitial += it1->Get4Momentum();
549 #endif 549 #endif
550 } 550 }
551 551
552 result=theModel->Propagate(initalStat 552 result=theModel->Propagate(initalState, target3dNucleus);
553 #ifdef debug_BLIR_finalstate 553 #ifdef debug_BLIR_finalstate
554 if( result && result->size()>0) 554 if( result && result->size()>0)
555 { 555 {
556 G4cout << " Cascade result " << G4endl; 556 G4cout << " Cascade result " << G4endl;
557 G4LorentzVector presult; 557 G4LorentzVector presult;
558 G4ReactionProductVector::iterator 558 G4ReactionProductVector::iterator iter;
559 G4ReactionProduct xp; 559 G4ReactionProduct xp;
560 for (iter=result->begin(); iter ! 560 for (iter=result->begin(); iter !=result->end(); ++iter)
561 { 561 {
562 presult += G4LorentzVector((*ite 562 presult += G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy());
563 G4cout << (*iter)->GetDefinition()->Ge 563 G4cout << (*iter)->GetDefinition()->GetParticleName() << " : "
564 << "("<< (*iter)->GetMomentum().x()<<" 564 << "("<< (*iter)->GetMomentum().x()<<","
565 << (*iter)->GetMomentum().y()<<"," 565 << (*iter)->GetMomentum().y()<<","
566 << (*iter)->GetMomentum().z()<<";" 566 << (*iter)->GetMomentum().z()<<";"
567 << (*iter)->GetTotalEnergy() <<")"< 567 << (*iter)->GetTotalEnergy() <<")"<< G4endl;
568 } 568 }
569 569
570 G4cout << "BLIC check result : in 570 G4cout << "BLIC check result : initial " << pinitial << " mass tgt " << target3dNucleus->GetMass()
571 << " final " << presult 571 << " final " << presult
572 << " IF - FF " << pinitial +G 572 << " IF - FF " << pinitial +G4LorentzVector(target3dNucleus->GetMass()) - presult << G4endl;
573 573
574 } 574 }
575 #endif 575 #endif
576 if( result && result->size()==0) 576 if( result && result->size()==0)
577 { 577 {
578 delete result; 578 delete result;
579 result=0; 579 result=0;
580 } 580 }
581 if ( ! result ) 581 if ( ! result )
582 { 582 {
583 delete target3dNucleus; 583 delete target3dNucleus;
584 delete projectile3dNucleus; 584 delete projectile3dNucleus;
585 } 585 }
586 586
587 // std::for_each(initalState->begin() 587 // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>());
588 // delete initalState; 588 // delete initalState;
589 589
590 } while (! result && tryCount< 150); / 590 } while (! result && tryCount< 150); /* Loop checking, 31.08.2015, G.Folger */
591 return result; 591 return result;
592 } 592 }
593 G4double G4BinaryLightIonReaction::GetProjecti 593 G4double G4BinaryLightIonReaction::GetProjectileExcitation()
594 { 594 {
595 595
596 G4Nucleon * aNuc; 596 G4Nucleon * aNuc;
597 // the projectileNucleus excitation ener 597 // the projectileNucleus excitation energy estimate...
598 G4double theStatisticalExEnergy = 0; 598 G4double theStatisticalExEnergy = 0;
599 projectile3dNucleus->StartLoop(); 599 projectile3dNucleus->StartLoop();
600 while( (aNuc=projectile3dNucleus->GetNex 600 while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */
601 { 601 {
602 //G4cout << " Nucleon : " << a 602 //G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl;
603 if(aNuc->AreYouHit()) { 603 if(aNuc->AreYouHit()) {
604 G4ThreeVector aPosition(aNuc->GetP 604 G4ThreeVector aPosition(aNuc->GetPosition());
605 G4double localDensity = projectile 605 G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition);
606 G4double localPfermi = theFermi.Ge 606 G4double localPfermi = theFermi.GetFermiMomentum(localDensity);
607 G4double nucMass = aNuc->GetDefini 607 G4double nucMass = aNuc->GetDefinition()->GetPDGMass();
608 G4double localFermiEnergy = std::s 608 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
609 G4double deltaE = localFermiEnergy 609 G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag());
610 theStatisticalExEnergy += deltaE; 610 theStatisticalExEnergy += deltaE;
611 } 611 }
612 } 612 }
613 return theStatisticalExEnergy; 613 return theStatisticalExEnergy;
614 } 614 }
615 615
616 G4LorentzVector G4BinaryLightIonReaction::Sort 616 G4LorentzVector G4BinaryLightIonReaction::SortResult(G4ReactionProductVector * result, G4ReactionProductVector * spectators,G4ReactionProductVector * cascaders)
617 { 617 {
618 unsigned int i(0); 618 unsigned int i(0);
619 spectatorA=spectatorZ=0; 619 spectatorA=spectatorZ=0;
620 G4LorentzVector pspectators(0,0,0,0); 620 G4LorentzVector pspectators(0,0,0,0);
621 pFinalState=G4LorentzVector(0,0,0,0); 621 pFinalState=G4LorentzVector(0,0,0,0);
622 for(i=0; i<result->size(); i++) 622 for(i=0; i<result->size(); i++)
623 { 623 {
624 if( (*result)[i]->GetNewlyAdded() ) 624 if( (*result)[i]->GetNewlyAdded() )
625 { 625 {
626 pFinalState += G4LorentzVector( (*res 626 pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
627 cascaders->push_back((*result)[i]); 627 cascaders->push_back((*result)[i]);
628 } 628 }
629 else { 629 else {
630 // G4cout <<" spectator ... 630 // G4cout <<" spectator ... ";
631 pspectators += G4LorentzVector( (*res 631 pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
632 spectators->push_back((*result)[i]); 632 spectators->push_back((*result)[i]);
633 spectatorA++; 633 spectatorA++;
634 spectatorZ+= G4lrint((*result)[i]->Ge 634 spectatorZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()/eplus);
635 } 635 }
636 636
637 // G4cout << (*result)[i]<< " " 637 // G4cout << (*result)[i]<< " "
638 // << (*result)[i]->GetDefinition 638 // << (*result)[i]->GetDefinition()->GetParticleName() << " "
639 // << (*result)[i]->GetMomentum() 639 // << (*result)[i]->GetMomentum()<< " "
640 // << (*result)[i]->GetTotalEnerg 640 // << (*result)[i]->GetTotalEnergy() << G4endl;
641 } 641 }
642 //G4cout << "pFinalState / pspectators, 642 //G4cout << "pFinalState / pspectators, (A,Z), p " << pFinalState << " / " << spectators->size()
643 // << " (" << spectatorA << ", "<< sp 643 // << " (" << spectatorA << ", "<< spectatorZ << "), 4-mom: " << pspectators << G4endl;
644 644
645 return pspectators; 645 return pspectators;
646 } 646 }
647 647
648 void G4BinaryLightIonReaction::DeExciteSpectat 648 void G4BinaryLightIonReaction::DeExciteSpectatorNucleus(G4ReactionProductVector * spectators, G4ReactionProductVector * cascaders,
649 649 G4double theStatisticalExEnergy, G4LorentzVector & pSpectators)
650 { 650 {
651 // call precompound model 651 // call precompound model
652 G4ReactionProductVector * proFrag = 0; 652 G4ReactionProductVector * proFrag = 0;
653 G4LorentzVector pFragment(0.,0.,0.,0.); 653 G4LorentzVector pFragment(0.,0.,0.,0.);
654 // G4cout << " == pre boost 1 "<< mome 654 // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
655 G4LorentzRotation boost_fragments; 655 G4LorentzRotation boost_fragments;
656 // G4cout << " == post boost 1 "<< mom 656 // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
657 // G4LorentzRotation boost_spectator_mom 657 // G4LorentzRotation boost_spectator_mom(-momentum.boostVector());
658 // G4cout << "- momentum " << boost_spe 658 // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl;
659 G4LorentzVector pFragments(0,0,0,0); 659 G4LorentzVector pFragments(0,0,0,0);
660 660
661 if(spectatorZ>0 && spectatorA>1) 661 if(spectatorZ>0 && spectatorA>1)
662 { 662 {
663 // Make the fragment 663 // Make the fragment
664 G4Fragment aProRes; 664 G4Fragment aProRes;
665 aProRes.SetZandA_asInt(spectatorZ, spect 665 aProRes.SetZandA_asInt(spectatorZ, spectatorA);
666 aProRes.SetNumberOfParticles(0); 666 aProRes.SetNumberOfParticles(0);
667 aProRes.SetNumberOfCharged(0); 667 aProRes.SetNumberOfCharged(0);
668 aProRes.SetNumberOfHoles(pA-spectatorA); 668 aProRes.SetNumberOfHoles(pA-spectatorA);
669 G4double mFragment=G4ParticleTable::GetP 669 G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(spectatorZ,spectatorA);
670 pFragment=G4LorentzVector(0,0,0,mFragmen 670 pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) );
671 aProRes.SetMomentum(pFragment); 671 aProRes.SetMomentum(pFragment);
672 672
673 proFrag = theHandler->BreakItUp(aProRes) 673 proFrag = theHandler->BreakItUp(aProRes);
674 674
675 boost_fragments = G4LorentzRotation(pSpe 675 boost_fragments = G4LorentzRotation(pSpectators.boostVector());
676 676
677 // G4cout << " Fragment a,z, Mass Fr 677 // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE "
678 // << spectatorA <<" "<< spectator 678 // << spectatorA <<" "<< spectatorZ <<" "<< mFragment <<" "
679 // << momentum.mag() <<" "<< momen 679 // << momentum.mag() <<" "<< momentum.mag() - mFragment
680 // << " "<<theStatisticalExEnergy 680 // << " "<<theStatisticalExEnergy
681 // << " "<< boost_fragments*pFragm 681 // << " "<< boost_fragments*pFragment<< G4endl;
682 G4ReactionProductVector::iterator ispect 682 G4ReactionProductVector::iterator ispectator;
683 for (ispectator=spectators->begin();ispe 683 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
684 { 684 {
685 delete *ispectator; 685 delete *ispectator;
686 } 686 }
687 } 687 }
688 else if(spectatorA!=0) 688 else if(spectatorA!=0)
689 { 689 {
690 G4ReactionProductVector::iterator ispecta 690 G4ReactionProductVector::iterator ispectator;
691 for (ispectator=spectators->begin();ispec 691 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
692 { 692 {
693 (*ispectator)->SetNewlyAdded(true); 693 (*ispectator)->SetNewlyAdded(true);
694 cascaders->push_back(*ispectator); 694 cascaders->push_back(*ispectator);
695 pFinalState+=G4LorentzVector((*ispect 695 pFinalState+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
696 //G4cout << "BLIC: spectator 696 //G4cout << "BLIC: spectatorA>0, Z=0 from spectator "
697 // << (*ispectator)->GetDefi 697 // << (*ispectator)->GetDefinition()->GetParticleName() << " "
698 // << (*ispectator)->GetMome 698 // << (*ispectator)->GetMomentum()<< " "
699 // << (*ispectator)->GetTota 699 // << (*ispectator)->GetTotalEnergy() << G4endl;
700 } 700 }
701 701
702 } 702 }
703 // / if (spectators) 703 // / if (spectators)
704 delete spectators; 704 delete spectators;
705 705
706 // collect the evaporation part and boost t 706 // collect the evaporation part and boost to spectator frame
707 G4ReactionProductVector::iterator ii; 707 G4ReactionProductVector::iterator ii;
708 if(proFrag) 708 if(proFrag)
709 { 709 {
710 for(ii=proFrag->begin(); ii!=proFrag->en 710 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
711 { 711 {
712 (*ii)->SetNewlyAdded(true); 712 (*ii)->SetNewlyAdded(true);
713 G4LorentzVector tmp((*ii)->GetMomentu 713 G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy());
714 tmp *= boost_fragments; 714 tmp *= boost_fragments;
715 (*ii)->SetMomentum(tmp.vect()); 715 (*ii)->SetMomentum(tmp.vect());
716 (*ii)->SetTotalEnergy(tmp.e()); 716 (*ii)->SetTotalEnergy(tmp.e());
717 // result->push_back(*ii); 717 // result->push_back(*ii);
718 pFragments += tmp; 718 pFragments += tmp;
719 } 719 }
720 } 720 }
721 721
722 // G4cout << "Fragmented p, momentum, de 722 // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum
723 // <<" "<< pFragments-momentum < 723 // <<" "<< pFragments-momentum << G4endl;
724 724
725 // correct p/E of Cascade secondaries 725 // correct p/E of Cascade secondaries
726 G4LorentzVector pCas=pInitialState - pFragm 726 G4LorentzVector pCas=pInitialState - pFragments;
727 727
728 //G4cout <<"BLIC: Going to correct from 728 //G4cout <<"BLIC: Going to correct from " << pFinalState << " to " << pCas << G4endl;
729 // the creation of excited fragment did vi 729 // the creation of excited fragment did violate E/p, so correct cascaders to get overall conservation.
730 G4bool EnergyIsCorrect=EnergyAndMomentumCor 730 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas);
731 if ( ! EnergyIsCorrect && debug_G4BinaryLig 731 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
732 { 732 {
733 G4cout << "G4BinaryLightIonReaction E/P 733 G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl;
734 } 734 }
735 735
736 // Add deexcitation secondaries 736 // Add deexcitation secondaries
737 if(proFrag) 737 if(proFrag)
738 { 738 {
739 for(ii=proFrag->begin(); ii!=proFrag->en 739 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
740 { 740 {
741 cascaders->push_back(*ii); 741 cascaders->push_back(*ii);
742 } 742 }
743 delete proFrag; 743 delete proFrag;
744 } 744 }
745 //G4cout << "EnergyIsCorrect? " << Energ 745 //G4cout << "EnergyIsCorrect? " << EnergyIsCorrect << G4endl;
746 if ( ! EnergyIsCorrect ) 746 if ( ! EnergyIsCorrect )
747 { 747 {
748 // G4cout <<" ! EnergyIsCorrect " << 748 // G4cout <<" ! EnergyIsCorrect " << pFinalState << " to " << pInitialState << G4endl;
749 if (! EnergyAndMomentumCorrector(cascade 749 if (! EnergyAndMomentumCorrector(cascaders,pInitialState))
750 { 750 {
751 if(debug_G4BinaryLightIonReactionResu 751 if(debug_G4BinaryLightIonReactionResults)
752 G4cout << "G4BinaryLightIonReactio 752 G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl;
753 } 753 }
754 } 754 }
755 755
756 } 756 }
757 757
758 758