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

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

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