1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | #include "G4BinaryLightIonReaction.hh" |
---|
27 | #include "G4LorentzVector.hh" |
---|
28 | #include "G4LorentzRotation.hh" |
---|
29 | #include <algorithm> |
---|
30 | #include "G4ReactionProductVector.hh" |
---|
31 | #include <vector> |
---|
32 | #include "G4ping.hh" |
---|
33 | #include "G4Delete.hh" |
---|
34 | #include "G4Neutron.hh" |
---|
35 | #include "G4VNuclearDensity.hh" |
---|
36 | #include "G4FermiMomentum.hh" |
---|
37 | #include "G4HadTmpUtil.hh" |
---|
38 | #include <cmath> |
---|
39 | |
---|
40 | G4BinaryLightIonReaction::G4BinaryLightIonReaction() |
---|
41 | : G4HadronicInteraction("Binary Cascade"), theModel() , |
---|
42 | theHandler(0) , theProjectileFragmentation(0) |
---|
43 | { |
---|
44 | theHandler= new G4ExcitationHandler; |
---|
45 | SetPrecompound(new G4PreCompoundModel(theHandler)); |
---|
46 | } |
---|
47 | |
---|
48 | G4HadFinalState *G4BinaryLightIonReaction:: |
---|
49 | ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus ) |
---|
50 | { |
---|
51 | static G4int eventcounter=0; |
---|
52 | eventcounter++; |
---|
53 | if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number starts ######### "<<eventcounter<<G4endl; |
---|
54 | G4ping debug("debug_G4BinaryLightIonReaction"); |
---|
55 | G4double a1=aTrack.GetDefinition()->GetBaryonNumber(); |
---|
56 | G4double z1=aTrack.GetDefinition()->GetPDGCharge(); |
---|
57 | G4double a2=targetNucleus.GetN(); |
---|
58 | G4double z2=targetNucleus.GetZ(); |
---|
59 | debug.push_back(a1); |
---|
60 | debug.push_back(z1); |
---|
61 | debug.push_back(a2); |
---|
62 | debug.push_back(z2); |
---|
63 | // debug.push_back(m2); |
---|
64 | G4LorentzVector mom(aTrack.Get4Momentum()); |
---|
65 | debug.push_back(mom); |
---|
66 | debug.dump(); |
---|
67 | G4LorentzRotation toBreit(mom.boostVector()); |
---|
68 | |
---|
69 | G4bool swapped = false; |
---|
70 | if(a2<a1) |
---|
71 | { |
---|
72 | debug.push_back("swapping...."); |
---|
73 | swapped = true; |
---|
74 | G4double tmp(0); |
---|
75 | tmp = a2; a2=a1; a1=tmp; |
---|
76 | tmp = z2; z2=z1; z1=tmp; |
---|
77 | G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(G4lrint(z1),G4lrint(a1)); |
---|
78 | G4LorentzVector it(m1, G4ThreeVector(0,0,0)); |
---|
79 | mom = toBreit*it; |
---|
80 | } |
---|
81 | debug.push_back("After swap"); |
---|
82 | debug.push_back(a1); |
---|
83 | debug.push_back(z1); |
---|
84 | debug.push_back(a2); |
---|
85 | debug.push_back(z2); |
---|
86 | debug.push_back(mom); |
---|
87 | debug.dump(); |
---|
88 | |
---|
89 | G4ReactionProductVector * result = NULL; |
---|
90 | G4ReactionProductVector * cascaders= new G4ReactionProductVector; |
---|
91 | G4double m_nucl(0); // to check energy balance |
---|
92 | |
---|
93 | |
---|
94 | // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(G4lrint(z1),G4lrint(a1)); |
---|
95 | // G4cout << "Entering the decision point " |
---|
96 | // << (mom.t()-mom.mag())/a1 << " " |
---|
97 | // << a1<<" "<< z1<<" " |
---|
98 | // << a2<<" "<< z2<<G4endl |
---|
99 | // << " "<<mom.t()-mom.mag()<<" " |
---|
100 | // << mom.t()- m1<<G4endl; |
---|
101 | if( (mom.t()-mom.mag())/a1 < 50*MeV ) |
---|
102 | { |
---|
103 | // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl; |
---|
104 | // m_nucl = mom.mag(); |
---|
105 | delete cascaders; |
---|
106 | G4Fragment aPreFrag; |
---|
107 | aPreFrag.SetA(a1+a2); |
---|
108 | aPreFrag.SetZ(z1+z2); |
---|
109 | aPreFrag.SetNumberOfParticles(G4lrint(a1)); |
---|
110 | aPreFrag.SetNumberOfCharged(G4lrint(z1)); |
---|
111 | aPreFrag.SetNumberOfHoles(0); |
---|
112 | G4ThreeVector plop(0.,0., mom.vect().mag()); |
---|
113 | G4double m2=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(G4lrint(z2),G4lrint(a2)); |
---|
114 | m_nucl=m2; |
---|
115 | G4LorentzVector aL(mom.t()+m2, plop); |
---|
116 | aPreFrag.SetMomentum(aL); |
---|
117 | G4ParticleDefinition * preFragDef; |
---|
118 | preFragDef = G4ParticleTable::GetParticleTable() |
---|
119 | ->FindIon(G4lrint(z1+z2),G4lrint(a1+a2),0,G4lrint(z1+z2)); |
---|
120 | aPreFrag.SetParticleDefinition(preFragDef); |
---|
121 | |
---|
122 | // G4cout << "Fragment INFO "<< a1+a2 <<" "<<z1+z2<<" " |
---|
123 | // << aL <<" "<<preFragDef->GetParticleName()<<G4endl; |
---|
124 | cascaders = theProjectileFragmentation->DeExcite(aPreFrag); |
---|
125 | G4double tSum = 0; |
---|
126 | for(size_t count = 0; count<cascaders->size(); count++) |
---|
127 | { |
---|
128 | cascaders->operator[](count)->SetNewlyAdded(true); |
---|
129 | tSum += cascaders->operator[](count)->GetKineticEnergy(); |
---|
130 | } |
---|
131 | // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl; |
---|
132 | } |
---|
133 | else |
---|
134 | { |
---|
135 | |
---|
136 | |
---|
137 | G4V3DNucleus * fancyNucleus = NULL; |
---|
138 | G4Fancy3DNucleus * projectile = NULL; |
---|
139 | G4double m1(0) ,m2(0); |
---|
140 | G4LorentzVector it; |
---|
141 | |
---|
142 | G4FermiMomentum theFermi; |
---|
143 | G4int tryCount(0); |
---|
144 | while(!result) |
---|
145 | { |
---|
146 | projectile = new G4Fancy3DNucleus; |
---|
147 | projectile->Init(a1, z1); |
---|
148 | projectile->CenterNucleons(); |
---|
149 | m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( |
---|
150 | projectile->GetCharge(),projectile->GetMassNumber()); |
---|
151 | it=toBreit * G4LorentzVector(m1,G4ThreeVector(0,0,0)); |
---|
152 | fancyNucleus = new G4Fancy3DNucleus; |
---|
153 | fancyNucleus->Init(a2, z2); |
---|
154 | m2=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( |
---|
155 | fancyNucleus->GetCharge(),fancyNucleus->GetMassNumber()); |
---|
156 | m_nucl = ( swapped ) ? m1 : m2; |
---|
157 | // G4cout << " mass table, nucleus, delta : " << m2 <<" "<< fancyNucleus->GetMass() |
---|
158 | // <<" "<<m2-fancyNucleus->GetMass() << G4endl; |
---|
159 | G4double impactMax = fancyNucleus->GetOuterRadius()+projectile->GetOuterRadius(); |
---|
160 | // G4cout << "out radius - nucleus - projectile " << fancyNucleus->GetOuterRadius()/fermi << " - " << projectile->GetOuterRadius()/fermi << G4endl; |
---|
161 | G4double aX=(2.*G4UniformRand()-1.)*impactMax; |
---|
162 | G4double aY=(2.*G4UniformRand()-1.)*impactMax; |
---|
163 | G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi); |
---|
164 | debug.push_back("Impact parameter"); |
---|
165 | debug.push_back(aX); |
---|
166 | debug.push_back(aY); |
---|
167 | debug.push_back(-2.*impactMax); |
---|
168 | debug.dump(); |
---|
169 | |
---|
170 | G4KineticTrackVector * initalState = new G4KineticTrackVector; |
---|
171 | projectile->StartLoop(); |
---|
172 | G4Nucleon * aNuc; |
---|
173 | G4LorentzVector tmpV(0,0,0,0); |
---|
174 | G4LorentzVector nucleonMom(1./a1*mom); |
---|
175 | nucleonMom.setZ(nucleonMom.vect().mag()); |
---|
176 | nucleonMom.setX(0); |
---|
177 | nucleonMom.setY(0); |
---|
178 | debug.push_back(" projectile nucleon momentum"); |
---|
179 | debug.push_back(nucleonMom); |
---|
180 | debug.dump(); |
---|
181 | theFermi.Init(a1,z1); |
---|
182 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
183 | { |
---|
184 | G4LorentzVector p4 = aNuc->GetMomentum(); |
---|
185 | tmpV+=p4; |
---|
186 | G4ThreeVector nucleonPosition(aNuc->GetPosition()); |
---|
187 | G4double density=(projectile->GetNuclearDensity())->GetDensity(nucleonPosition); |
---|
188 | nucleonPosition += pos; |
---|
189 | G4KineticTrack * it = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom ); |
---|
190 | it->SetState(G4KineticTrack::outside); |
---|
191 | G4double pfermi= theFermi.GetFermiMomentum(density); |
---|
192 | G4double mass = aNuc->GetDefinition()->GetPDGMass(); |
---|
193 | G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass; |
---|
194 | it->SetProjectilePotential(-Efermi); |
---|
195 | initalState->push_back(it); |
---|
196 | } |
---|
197 | debug.push_back(" Sum of proj. nucleon momentum"); |
---|
198 | debug.push_back(tmpV); |
---|
199 | debug.dump(); |
---|
200 | |
---|
201 | result=theModel.Propagate(initalState, fancyNucleus); |
---|
202 | debug.push_back("################# Result size"); |
---|
203 | if (result) { |
---|
204 | debug.push_back(result->size()); |
---|
205 | } else { |
---|
206 | debug.push_back(" -none-"); |
---|
207 | } |
---|
208 | debug.dump(); |
---|
209 | // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>()); |
---|
210 | // delete initalState; |
---|
211 | |
---|
212 | if(! result || result->size()==0) |
---|
213 | { |
---|
214 | if (result) {delete result; result=0;} |
---|
215 | delete fancyNucleus; |
---|
216 | delete projectile; |
---|
217 | if (++tryCount > 200) |
---|
218 | { |
---|
219 | // abort!! |
---|
220 | |
---|
221 | G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; |
---|
222 | G4cerr << " Primary " << aTrack.GetDefinition() |
---|
223 | << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() |
---|
224 | << "," << aTrack.GetDefinition()->GetPDGCharge() << ") " |
---|
225 | << ", kinetic energy " << aTrack.GetKineticEnergy() |
---|
226 | << G4endl; |
---|
227 | G4cerr << " Target nucleus (A,Z)=(" << targetNucleus.GetN() |
---|
228 | << "," << targetNucleus.GetZ() << G4endl; |
---|
229 | G4cerr << " if frequent, please submit above information as bug report" |
---|
230 | << G4endl << G4endl; |
---|
231 | |
---|
232 | theResult.Clear(); |
---|
233 | theResult.SetStatusChange(isAlive); |
---|
234 | theResult.SetEnergyChange(aTrack.GetKineticEnergy()); |
---|
235 | theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); |
---|
236 | return &theResult; |
---|
237 | |
---|
238 | } |
---|
239 | } |
---|
240 | else |
---|
241 | { |
---|
242 | break; |
---|
243 | } |
---|
244 | } |
---|
245 | debug.push_back(" Attempts to create final state"); |
---|
246 | debug.push_back(tryCount); |
---|
247 | debug.dump(); |
---|
248 | debug.push_back("################# Through the loop ? "); debug.dump(); |
---|
249 | //inverse transformation in case we swapped. |
---|
250 | G4int resA(0), resZ(0); |
---|
251 | G4Nucleon * aNuc; |
---|
252 | // fancyNucleus->StartLoop(); |
---|
253 | // while( (aNuc=fancyNucleus->GetNextNucleon()) ) |
---|
254 | // { |
---|
255 | // G4cout << " tgt Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl; |
---|
256 | // } |
---|
257 | G4ReactionProductVector * spectators= new G4ReactionProductVector; |
---|
258 | debug.push_back("getting at the hits"); debug.dump(); |
---|
259 | // the projectile excitation energy estimate... |
---|
260 | G4double theStatisticalExEnergy = 0; |
---|
261 | projectile->StartLoop(); |
---|
262 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
263 | { |
---|
264 | // G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl; |
---|
265 | debug.push_back("getting the hits"); debug.dump(); |
---|
266 | if(!aNuc->AreYouHit()) |
---|
267 | { |
---|
268 | resA++; |
---|
269 | resZ+=G4lrint(aNuc->GetDefinition()->GetPDGCharge()); |
---|
270 | } |
---|
271 | else |
---|
272 | { |
---|
273 | debug.push_back(" ##### a hit ##### "); debug.dump(); |
---|
274 | G4ThreeVector aPosition(aNuc->GetPosition()); |
---|
275 | G4double localDensity = projectile->GetNuclearDensity()->GetDensity(aPosition); |
---|
276 | G4double localPfermi = theFermi.GetFermiMomentum(localDensity); |
---|
277 | G4double nucMass = aNuc->GetDefinition()->GetPDGMass(); |
---|
278 | G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass; |
---|
279 | G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag()); |
---|
280 | theStatisticalExEnergy += deltaE; |
---|
281 | } |
---|
282 | debug.push_back("collected a hit"); |
---|
283 | debug.push_back(aNuc->GetMomentum()); |
---|
284 | debug.dump(); |
---|
285 | } |
---|
286 | delete fancyNucleus; |
---|
287 | delete projectile; |
---|
288 | G4ping debug("debug_G4BinaryLightIonReaction_1"); |
---|
289 | debug.push_back("have the hits. A,Z, excitE"); |
---|
290 | debug.push_back(resA); |
---|
291 | debug.push_back(resZ); |
---|
292 | debug.push_back(theStatisticalExEnergy); |
---|
293 | debug.dump(); |
---|
294 | // Calculate excitation energy |
---|
295 | G4LorentzVector iState = mom; |
---|
296 | iState.setT(iState.getT()+m2); |
---|
297 | |
---|
298 | G4LorentzVector fState(0,0,0,0); |
---|
299 | G4LorentzVector pspectators(0,0,0,0); |
---|
300 | unsigned int i(0); |
---|
301 | // G4int spectA(0),spectZ(0); |
---|
302 | for(i=0; i<result->size(); i++) |
---|
303 | { |
---|
304 | if( (*result)[i]->GetNewlyAdded() ) |
---|
305 | { |
---|
306 | fState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); |
---|
307 | cascaders->push_back((*result)[i]); |
---|
308 | // G4cout <<" secondary ... "; |
---|
309 | debug.push_back("secondary "); |
---|
310 | debug.push_back((*result)[i]->GetDefinition()->GetParticleName()); |
---|
311 | debug.push_back(G4LorentzVector((*result)[i]->GetMomentum(),(*result)[i]->GetTotalEnergy())); |
---|
312 | debug.dump(); |
---|
313 | } |
---|
314 | else { |
---|
315 | // G4cout <<" spectator ... "; |
---|
316 | pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); |
---|
317 | spectators->push_back((*result)[i]); |
---|
318 | debug.push_back("spectator "); |
---|
319 | debug.push_back((*result)[i]->GetDefinition()->GetParticleName()); |
---|
320 | debug.push_back(G4LorentzVector((*result)[i]->GetMomentum(),(*result)[i]->GetTotalEnergy())); |
---|
321 | debug.dump(); |
---|
322 | // spectA++; |
---|
323 | // spectZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()); |
---|
324 | } |
---|
325 | |
---|
326 | // G4cout << (*result)[i]<< " " |
---|
327 | // << (*result)[i]->GetDefinition()->GetParticleName() << " " |
---|
328 | // << (*result)[i]->GetMomentum()<< " " |
---|
329 | // << (*result)[i]->GetTotalEnergy() << G4endl; |
---|
330 | } |
---|
331 | // if ( spectA-resA !=0 || spectZ-resZ !=0) |
---|
332 | // { |
---|
333 | // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" << |
---|
334 | // resA <<" "<< resZ <<" ; " << spectA <<" "<< spectZ << G4endl; |
---|
335 | // } |
---|
336 | delete result; |
---|
337 | debug.push_back(" iState - (fState+pspectators) "); |
---|
338 | debug.push_back(iState-fState-pspectators); |
---|
339 | debug.dump(); |
---|
340 | G4LorentzVector momentum(iState-fState); |
---|
341 | G4int loopcount(0); |
---|
342 | while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) |
---|
343 | { |
---|
344 | debug.push_back("the momentum balance"); |
---|
345 | debug.push_back(iState); |
---|
346 | debug.push_back(fState); |
---|
347 | debug.push_back(momentum-pspectators); |
---|
348 | debug.push_back(momentum); |
---|
349 | debug.dump(); |
---|
350 | G4LorentzVector pCorrect(iState-pspectators); |
---|
351 | G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect); |
---|
352 | if ( ! EnergyIsCorrect && getenv("debug_G4BinaryLightIonReactionResults")) |
---|
353 | { |
---|
354 | G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl; |
---|
355 | } |
---|
356 | fState=G4LorentzVector(); |
---|
357 | for(i=0; i<cascaders->size(); i++) |
---|
358 | { |
---|
359 | fState += G4LorentzVector( (*cascaders)[i]->GetMomentum(), (*cascaders)[i]->GetTotalEnergy() ); |
---|
360 | } |
---|
361 | momentum=iState-fState; |
---|
362 | debug.push_back("the momentum balance after correction"); |
---|
363 | debug.push_back(iState); |
---|
364 | debug.push_back(fState); |
---|
365 | debug.push_back(momentum-pspectators); |
---|
366 | debug.push_back(momentum); |
---|
367 | debug.dump(); |
---|
368 | if (++loopcount > 10 ) |
---|
369 | { |
---|
370 | if ( momentum.vect().mag() > momentum.e() ) |
---|
371 | { |
---|
372 | G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl; |
---|
373 | throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()"); |
---|
374 | } else { |
---|
375 | break; |
---|
376 | } |
---|
377 | |
---|
378 | } |
---|
379 | } |
---|
380 | |
---|
381 | |
---|
382 | |
---|
383 | |
---|
384 | // call precompound model |
---|
385 | G4ReactionProductVector * proFrag = NULL; |
---|
386 | G4LorentzVector pFragment; |
---|
387 | // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; |
---|
388 | G4LorentzRotation boost_fragments; |
---|
389 | // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; |
---|
390 | // G4LorentzRotation boost_spectator_mom(-momentum.boostVector()); |
---|
391 | // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl; |
---|
392 | G4LorentzVector pFragments(0); |
---|
393 | if(resZ>0 && resA>1) |
---|
394 | { |
---|
395 | // Make the fragment |
---|
396 | G4Fragment aProRes; |
---|
397 | aProRes.SetA(resA); |
---|
398 | aProRes.SetZ(resZ); |
---|
399 | aProRes.SetNumberOfParticles(0); |
---|
400 | aProRes.SetNumberOfCharged(0); |
---|
401 | aProRes.SetNumberOfHoles(G4lrint(a1)-resA); |
---|
402 | G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(resZ,resA); |
---|
403 | G4LorentzVector pFragment(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) ); |
---|
404 | aProRes.SetMomentum(pFragment); |
---|
405 | G4ParticleDefinition * resDef; |
---|
406 | resDef = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ); |
---|
407 | aProRes.SetParticleDefinition(resDef); |
---|
408 | proFrag = theHandler->BreakItUp(aProRes); |
---|
409 | if ( momentum.vect().mag() > momentum.e() ) |
---|
410 | { |
---|
411 | G4cerr << "mom check: " << momentum |
---|
412 | << " 3.mag "<< momentum.vect().mag() << G4endl |
---|
413 | << " .. iState/fState/spectators " << iState <<" " |
---|
414 | << fState << " " << pspectators << G4endl |
---|
415 | << " .. A,Z " << resA <<" "<< resZ << G4endl; |
---|
416 | G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; |
---|
417 | G4cerr << " Primary " << aTrack.GetDefinition() |
---|
418 | << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() |
---|
419 | << "," << aTrack.GetDefinition()->GetPDGCharge() << ") " |
---|
420 | << ", kinetic energy " << aTrack.GetKineticEnergy() |
---|
421 | << G4endl; |
---|
422 | G4cerr << " Target nucleus (A,Z)=(" << targetNucleus.GetN() |
---|
423 | << "," << targetNucleus.GetZ() << G4endl; |
---|
424 | G4cerr << " if frequent, please submit above information as bug report" |
---|
425 | << G4endl << G4endl; |
---|
426 | |
---|
427 | theResult.Clear(); |
---|
428 | theResult.SetStatusChange(isAlive); |
---|
429 | theResult.SetEnergyChange(aTrack.GetKineticEnergy()); |
---|
430 | theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); |
---|
431 | return &theResult; |
---|
432 | } |
---|
433 | |
---|
434 | G4LorentzRotation boost_fragments_here(momentum.boostVector()); |
---|
435 | boost_fragments = boost_fragments_here; |
---|
436 | |
---|
437 | // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE " |
---|
438 | // << resA <<" "<< resZ <<" "<< mFragment <<" " |
---|
439 | // << momentum.mag() <<" "<< momentum.mag() - mFragment |
---|
440 | // << " "<<theStatisticalExEnergy |
---|
441 | // << " "<< boost_fragments*pFragment<< G4endl; |
---|
442 | G4ReactionProductVector::iterator ispectator; |
---|
443 | for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) |
---|
444 | { |
---|
445 | delete *ispectator; |
---|
446 | } |
---|
447 | } |
---|
448 | else if(resA!=0) |
---|
449 | { |
---|
450 | G4ReactionProductVector::iterator ispectator; |
---|
451 | for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) |
---|
452 | { |
---|
453 | (*ispectator)->SetNewlyAdded(true); |
---|
454 | cascaders->push_back(*ispectator); |
---|
455 | pFragments+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy()); |
---|
456 | // G4cout << "from spectator " |
---|
457 | // << (*ispectator)->GetDefinition()->GetParticleName() << " " |
---|
458 | // << (*ispectator)->GetMomentum()<< " " |
---|
459 | // << (*ispectator)->GetTotalEnergy() << G4endl; |
---|
460 | } |
---|
461 | } |
---|
462 | if (spectators) delete spectators; |
---|
463 | |
---|
464 | // collect the evaporation part |
---|
465 | debug.push_back("the nucleon count balance"); |
---|
466 | debug.push_back(resA); |
---|
467 | debug.push_back(resZ); |
---|
468 | if(proFrag) debug.push_back(proFrag->size()); |
---|
469 | debug.dump(); |
---|
470 | G4ReactionProductVector::iterator ii; |
---|
471 | if(proFrag) for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) |
---|
472 | { |
---|
473 | (*ii)->SetNewlyAdded(true); |
---|
474 | G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy()); |
---|
475 | tmp *= boost_fragments; |
---|
476 | (*ii)->SetMomentum(tmp.vect()); |
---|
477 | (*ii)->SetTotalEnergy(tmp.e()); |
---|
478 | // result->push_back(*ii); |
---|
479 | pFragments += tmp; |
---|
480 | } |
---|
481 | |
---|
482 | // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum |
---|
483 | // <<" "<< pFragments-momentum << G4endl; |
---|
484 | debug.push_back("################# done with evaporation"); debug.dump(); |
---|
485 | |
---|
486 | // correct p/E of Cascade secondaries |
---|
487 | G4LorentzVector pCas=iState - pFragments; |
---|
488 | // G4cout <<" Going to correct from " << fState << " to " << pCas << G4endl; |
---|
489 | G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas); |
---|
490 | if ( ! EnergyIsCorrect ) |
---|
491 | { |
---|
492 | if(getenv("debug_G4BinaryLightIonReactionResults")) |
---|
493 | G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl; |
---|
494 | } |
---|
495 | |
---|
496 | // Add deexcitation secondaries |
---|
497 | if(proFrag) for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) |
---|
498 | { |
---|
499 | cascaders->push_back(*ii); |
---|
500 | } |
---|
501 | if (proFrag) delete proFrag; |
---|
502 | |
---|
503 | if ( ! EnergyIsCorrect ) |
---|
504 | { |
---|
505 | if (! EnergyAndMomentumCorrector(cascaders,iState)) |
---|
506 | { |
---|
507 | if(getenv("debug_G4BinaryLightIonReactionResults")) |
---|
508 | G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl; |
---|
509 | } |
---|
510 | } |
---|
511 | |
---|
512 | } |
---|
513 | |
---|
514 | // Rotate to lab |
---|
515 | G4LorentzRotation toZ; |
---|
516 | toZ.rotateZ(-1*mom.phi()); |
---|
517 | toZ.rotateY(-1*mom.theta()); |
---|
518 | G4LorentzRotation toLab(toZ.inverse()); |
---|
519 | |
---|
520 | // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped. |
---|
521 | // theResult.Clear(); |
---|
522 | theResult.Clear(); |
---|
523 | theResult.SetStatusChange(stopAndKill); |
---|
524 | G4double Etot(0); |
---|
525 | size_t i=0; |
---|
526 | for(i=0; i<cascaders->size(); i++) |
---|
527 | { |
---|
528 | if((*cascaders)[i]->GetNewlyAdded()) |
---|
529 | { |
---|
530 | G4DynamicParticle * aNew = |
---|
531 | new G4DynamicParticle((*cascaders)[i]->GetDefinition(), |
---|
532 | (*cascaders)[i]->GetTotalEnergy(), |
---|
533 | (*cascaders)[i]->GetMomentum() ); |
---|
534 | G4LorentzVector tmp = aNew->Get4Momentum(); |
---|
535 | if(swapped) |
---|
536 | { |
---|
537 | tmp*=toBreit.inverse(); |
---|
538 | tmp.setVect(-tmp.vect()); |
---|
539 | } |
---|
540 | tmp *= toLab; |
---|
541 | aNew->Set4Momentum(tmp); |
---|
542 | theResult.AddSecondary(aNew); |
---|
543 | Etot += tmp.e(); |
---|
544 | // G4cout << "LIBIC: Secondary " << aNew->GetDefinition()->GetParticleName() |
---|
545 | // <<" "<< aNew->GetMomentum() |
---|
546 | // <<" "<< aNew->GetTotalEnergy() |
---|
547 | // << G4endl; |
---|
548 | } |
---|
549 | delete (*cascaders)[i]; |
---|
550 | } |
---|
551 | if(cascaders) delete cascaders; |
---|
552 | |
---|
553 | G4ping debug1("debug_G4BinaryLightIonReactionResults"); |
---|
554 | debug1.push_back("Result analysis, secondaries"); |
---|
555 | debug1.push_back(theResult.GetNumberOfSecondaries()); |
---|
556 | debug1.dump(); |
---|
557 | debug1.push_back(" Energy conservation initial/final/delta(init-final) "); |
---|
558 | debug1.push_back(aTrack.GetTotalEnergy() + m_nucl); |
---|
559 | debug1.push_back(aTrack.GetTotalEnergy()); |
---|
560 | debug1.push_back(m_nucl); |
---|
561 | debug1.push_back(Etot); |
---|
562 | debug1.push_back(aTrack.GetTotalEnergy() + m_nucl - Etot); |
---|
563 | debug1.dump(); |
---|
564 | |
---|
565 | if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### "<<eventcounter<<G4endl; |
---|
566 | |
---|
567 | return &theResult; |
---|
568 | } |
---|
569 | |
---|
570 | //**************************************************************************** |
---|
571 | G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector( |
---|
572 | G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom) |
---|
573 | //**************************************************************************** |
---|
574 | { |
---|
575 | const int nAttemptScale = 2500; |
---|
576 | const double ErrLimit = 1.E-6; |
---|
577 | if (Output->empty()) |
---|
578 | return TRUE; |
---|
579 | G4LorentzVector SumMom(0); |
---|
580 | G4double SumMass = 0; |
---|
581 | G4double TotalCollisionMass = TotalCollisionMom.m(); |
---|
582 | size_t i = 0; |
---|
583 | // Calculate sum hadron 4-momenta and summing hadron mass |
---|
584 | for(i = 0; i < Output->size(); i++) |
---|
585 | { |
---|
586 | SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
587 | SumMass += (*Output)[i]->GetDefinition()->GetPDGMass(); |
---|
588 | } |
---|
589 | // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass " |
---|
590 | // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl; |
---|
591 | if (SumMass > TotalCollisionMass) return FALSE; |
---|
592 | SumMass = SumMom.m2(); |
---|
593 | if (SumMass < 0) return FALSE; |
---|
594 | SumMass = std::sqrt(SumMass); |
---|
595 | |
---|
596 | // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s. |
---|
597 | G4ThreeVector Beta = -SumMom.boostVector(); |
---|
598 | // G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl; |
---|
599 | //--old Output->Boost(Beta); |
---|
600 | for(i = 0; i < Output->size(); i++) |
---|
601 | { |
---|
602 | G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
603 | mom *= Beta; |
---|
604 | (*Output)[i]->SetMomentum(mom.vect()); |
---|
605 | (*Output)[i]->SetTotalEnergy(mom.e()); |
---|
606 | } |
---|
607 | |
---|
608 | // Scale total c.m.s. hadron energy (hadron system mass). |
---|
609 | // It should be equal interaction mass |
---|
610 | G4double Scale = 0,OldScale=0; |
---|
611 | G4double factor = 1.; |
---|
612 | G4int cAttempt = 0; |
---|
613 | G4double Sum = 0; |
---|
614 | G4bool success = false; |
---|
615 | for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++) |
---|
616 | { |
---|
617 | Sum = 0; |
---|
618 | for(i = 0; i < Output->size(); i++) |
---|
619 | { |
---|
620 | G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
621 | HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect()); |
---|
622 | G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass())); |
---|
623 | HadronMom.setE(E); |
---|
624 | (*Output)[i]->SetMomentum(HadronMom.vect()); |
---|
625 | (*Output)[i]->SetTotalEnergy(HadronMom.e()); |
---|
626 | Sum += E; |
---|
627 | } |
---|
628 | OldScale=Scale; |
---|
629 | Scale = TotalCollisionMass/Sum - 1; |
---|
630 | // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl; |
---|
631 | if (std::abs(Scale) <= ErrLimit |
---|
632 | || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below |
---|
633 | { |
---|
634 | if (getenv("debug_G4BinaryLightIonReactionResults")) G4cout << "E/p corrector: " << cAttempt << G4endl; |
---|
635 | success = true; |
---|
636 | break; |
---|
637 | } |
---|
638 | if ( cAttempt > 10 ) |
---|
639 | { |
---|
640 | // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl; |
---|
641 | factor=std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale)))); |
---|
642 | // G4cout << " ? factor ? " << factor << G4endl; |
---|
643 | } |
---|
644 | } |
---|
645 | |
---|
646 | if( (!success) && getenv("debug_G4BinaryLightIonReactionResults")) |
---|
647 | { |
---|
648 | G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl; |
---|
649 | G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl; |
---|
650 | G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl; |
---|
651 | } |
---|
652 | |
---|
653 | // Compute c.m.s. interaction velocity and KTV back boost |
---|
654 | Beta = TotalCollisionMom.boostVector(); |
---|
655 | //--old Output->Boost(Beta); |
---|
656 | for(i = 0; i < Output->size(); i++) |
---|
657 | { |
---|
658 | G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
659 | mom *= Beta; |
---|
660 | (*Output)[i]->SetMomentum(mom.vect()); |
---|
661 | (*Output)[i]->SetTotalEnergy(mom.e()); |
---|
662 | } |
---|
663 | return TRUE; |
---|
664 | } |
---|