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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

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