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