source: trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc@ 1036

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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