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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 25.4 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
40G4BinaryLightIonReaction::G4BinaryLightIonReaction()
41 : G4HadronicInteraction("Binary Cascade"), theModel() ,
42 theHandler(0) , theProjectileFragmentation(0)
43{
44 theHandler= new G4ExcitationHandler;
45 SetPrecompound(new G4PreCompoundModel(theHandler));
46}
47
48G4HadFinalState *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//****************************************************************************
571G4bool 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 }
Note: See TracBrowser for help on using the repository browser.