| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include <algorithm>
27 #include <vector>
28 #include <cmath>
29 #include <numeric>
30
31 #include "G4BinaryLightIonReaction.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LorentzVector.hh"
35 #include "G4LorentzRotation.hh"
36 #include "G4ReactionProductVector.hh"
37 #include "G4ping.hh"
38 #include "G4Delete.hh"
39 #include "G4Neutron.hh"
40 #include "G4VNuclearDensity.hh"
41 #include "G4FermiMomentum.hh"
42 #include "G4HadTmpUtil.hh"
43 #include "G4PreCompoundModel.hh"
44 #include "G4HadronicInteractionRegistry.hh"
45 #include "G4Log.hh"
46 #include "G4PhysicsModelCatalog.hh"
47 #include "G4HadronicParameters.hh"
48
49 G4int G4BinaryLightIonReaction::theBLIR_ID = -
50
51 //#define debug_G4BinaryLightIonReaction
52 //#define debug_BLIR_finalstate
53 //#define debug_BLIR_result
54
55 G4BinaryLightIonReaction::G4BinaryLightIonReac
56 : G4HadronicInteraction("Binary Light Ion Casc
57 theProjectileFragmentation(ptr),
58 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spect
59 projectile3dNucleus(0),target3dNucleus(0)
60 {
61 if(!ptr) {
62 G4HadronicInteraction* p =
63 G4HadronicInteractionRegistry::Instance(
64 G4VPreCompoundModel* pre = static_cast<G4V
65 if(!pre) { pre = new G4PreCompoundModel();
66 theProjectileFragmentation = pre;
67 }
68 theModel = new G4BinaryCascade(theProjectile
69 theHandler = theProjectileFragmentation->Get
70 theBLIR_ID = G4PhysicsModelCatalog::GetM
71 debug_G4BinaryLightIonReactionResults = G4Ha
72 }
73
74 G4BinaryLightIonReaction::~G4BinaryLightIonRea
75 {}
76
77 void G4BinaryLightIonReaction::ModelDescriptio
78 {
79 outFile << "G4Binary Light Ion Cascade is an
80 << "using G4BinaryCasacde to model the i
81 << "nucleus with a nucleus.\n"
82 << "The lighter of the two nuclei is tre
83 << "which are transported simultaneously
84 }
85
86 //--------------------------------------------
87 struct ReactionProduct4Mom
88 {
89 G4LorentzVector operator()(G4LorentzVector
90 };
91
92 G4HadFinalState *G4BinaryLightIonReaction::
93 ApplyYourself(const G4HadProjectile &aTrack, G
94 {
95 if(debug_G4BinaryLightIonReactionResults) G4
96 G4ping debug("debug_G4BinaryLightIonReaction
97 pA=aTrack.GetDefinition()->GetBaryonNumber()
98 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCha
99 tA=targetNucleus.GetA_asInt();
100 tZ=targetNucleus.GetZ_asInt();
101 G4double timePrimary = aTrack.GetGlobalTime(
102 G4LorentzVector mom(aTrack.Get4Momentum());
103 //G4cout << "proj mom : " << mom << G4endl;
104 G4LorentzRotation toBreit(mom.boostVector())
105
106 G4bool swapped=SetLighterAsProjectile(mom, t
107 //G4cout << "after swap, swapped? / mom " <
108 G4ReactionProductVector * result = 0;
109 G4ReactionProductVector * cascaders=0; //new
110 // G4double m_nucl(0); // to check energ
111
112
113 // G4double m1=G4ParticleTable::GetPartic
114 // G4cout << "Entering the decision point
115 // << (mom.t()-mom.mag())/pA << "
116 // << pA<<" "<< pZ<<" "
117 // << tA<<" "<< tZ<<G4endl
118 // << " "<<mom.t()-mom.mag()<<" "
119 // << mom.t()- m1<<G4endl;
120 if( (mom.t()-mom.mag())/pA < 50*MeV )
121 {
122 // G4cout << "Using pre-compound only
123 // m_nucl = mom.mag();
124 cascaders=FuseNucleiAndPrompound(mom);
125 if( !cascaders )
126 {
127
128 // abort!! happens for too low e
129
130 theResult.Clear();
131 theResult.SetStatusChange(isAliv
132 theResult.SetEnergyChange(aTrack
133 theResult.SetMomentumChange(aTra
134 return &theResult;
135 }
136 }
137 else
138 {
139 result=Interact(mom,toBreit);
140
141 if(! result )
142 {
143 // abort!!
144
145 G4cerr << "G4BinaryLightIonReaction
146 G4cerr << " Primary " << aTrack.Get
147 << ", (A,Z)=(" << aTrack.GetDefi
148 << "," << aTrack.GetDefinition()
149 << ", kinetic energy " << aTrack
150 << G4endl;
151 G4cerr << " Target nucleus (A,Z)=("
152 << (swapped?pA:tA) << ","
153 << (swapped?pZ:tZ) << ")" <<
154 G4cerr << " if frequent, please sub
155 << G4endl << G4endl;
156
157 theResult.Clear();
158 theResult.SetStatusChange(isAlive);
159 theResult.SetEnergyChange(aTrack.Ge
160 theResult.SetMomentumChange(aTrack.
161 return &theResult;
162 }
163
164 // Calculate excitation energy,
165 G4double theStatisticalExEnergy = GetProj
166
167
168 pInitialState = mom;
169 //G4cout << "BLIC: pInitialState from
170 pInitialState.setT(pInitialState.getT() +
171 G4ParticleTable::GetParticleTable()->GetIo
172 //G4cout << "BLIC: target nucleus added
173
174 delete target3dNucleus;target3dNucleus=0;
175 delete projectile3dNucleus;projectile3dNu
176
177 G4ReactionProductVector * spectators= new
178
179 cascaders = new G4ReactionProductVector;
180
181 G4LorentzVector pspectators=SortResult(re
182 // this also sets spectatorA and
183
184 // pFinalState=std::accumulate(casca
185
186 std::vector<G4ReactionProduct *>::iterato
187
188 // G4cout << "pInitialState, pFin
189 // if ( spectA-spectatorA !=0 || spec
190 // {
191 // G4cout << "spect Nucl != spect
192 // spectatorA <<" "<< spectatorZ <<
193 // }
194 delete result;
195 result=0;
196 G4LorentzVector momentum(pInitialState-pF
197 G4int loopcount(0);
198 //G4cout << "BLIC: momentum, pspectato
199 while (std::abs(momentum.e()-pspectators.
200
201 {
202 G4LorentzVector pCorrect(pInitialState-
203 //G4cout << "BLIC:: BIC nonconservatio
204 // Correct outgoing casacde particles..
205 G4bool EnergyIsCorrect=EnergyAndMomentu
206 if ( ! EnergyIsCorrect && debug_G4Binar
207 {
208 G4cout << "Warning - G4BinaryLightIon
209 }
210 pFinalState=G4LorentzVector(0,0,0,0);
211 for(iter=cascaders->begin(); iter!=casc
212 {
213 pFinalState += G4LorentzVector( (*ite
214 }
215 momentum=pInitialState-pFinalState;
216 if (++loopcount > 10 )
217 {
218 break;
219 }
220 }
221
222 // Check if Energy/Momemtum is now ok, if
223 if ( std::abs(momentum.e()-pspectators.e(
224 {
225 for (iter=spectators->begin();iter!=specta
226 {
227 delete *iter;
228 }
229 delete spectators;
230 for(iter=cascaders->begin(); iter!=cascad
231 {
232 delete *iter;
233 }
234 delete cascaders;
235
236 G4cout << "G4BinaryLightIonReaction.cc: m
237 << " initial - final " << momentum
238 << " .. pInitialState/pFinalState/s
239 << pInitialState << G4endl
240 << pFinalState << G4endl
241 << pspectators << G4endl
242 << " .. A,Z " << spectatorA <<" "<<
243 G4cout << "G4BinaryLightIonReaction inval
244 G4cout << " Primary " << aTrack.GetDefini
245 << ", (A,Z)=(" << aTrack.GetDefin
246 << "," << aTrack.GetDefinition()-
247 << ", kinetic energy " << aTrack.
248 << G4endl;
249 G4cout << " Target nucleus (A,Z)=(" << t
250 << "," << targetNucleus.GetZ_
251 G4cout << " if frequent, please submit ab
252 << G4endl << G4endl;
253 #ifdef debug_G4BinaryLightIonReaction
254 G4ExceptionDescription ed;
255 ed << "G4BinaryLightIonreaction: Ter
256 G4Exception("G4BinaryLightIonreactio
257 ed);
258
259 #endif
260 theResult.Clear();
261 theResult.SetStatusChange(isAlive);
262 theResult.SetEnergyChange(aTrack.GetKinet
263 theResult.SetMomentumChange(aTrack.Get4Mo
264 return &theResult;
265
266 }
267 if (spectatorA > 0 )
268 {
269 // DeExciteSpectatorNucleus() also
270 DeExciteSpectatorNucleus(specta
271 } else { // no spectators
272 delete spectators;
273 }
274 }
275 // Rotate to lab
276 G4LorentzRotation toZ;
277 toZ.rotateZ(-1*mom.phi());
278 toZ.rotateY(-1*mom.theta());
279 G4LorentzRotation toLab(toZ.inverse());
280
281 // Fill the particle change, while rotating.
282 // theResult.Clear();
283 theResult.Clear();
284 theResult.SetStatusChange(stopAndKill);
285 G4LorentzVector ptot(0);
286 #ifdef debug_BLIR_result
287 G4LorentzVector p_raw;
288 #endif
289 //G4int i=0;
290
291 G4ReactionProductVector::iterator iter
292 for(iter=cascaders->begin(); iter!=cascaders
293 {
294 if((*iter)->GetNewlyAdded())
295 {
296 G4DynamicParticle * aNewDP =
297 new G4DynamicParticle((*iter)->GetDe
298 (*iter)->GetTotalEnergy(),
299 (*iter)->GetMomentum() );
300 G4LorentzVector tmp = aNewDP->Get4Moment
301 #ifdef debug_BLIR_result
302 p_raw+= tmp;
303 #endif
304 if(swapped)
305 {
306 tmp*=toBreit.inverse();
307 tmp.setVect(-tmp.vect());
308 }
309 tmp *= toLab;
310 aNewDP->Set4Momentum(tmp);
311 G4HadSecondary aNew = G4HadSecondary(aNe
312 G4double time = 0;
313 //if(time < 0.0) { time = 0.0; }
314 aNew.SetTime(timePrimary + time);
315 //aNew.SetCreatorModelID((*iter)->
316 aNew.SetCreatorModelID(theBLIR_ID)
317
318 theResult.AddSecondary(aNew);
319 ptot += tmp;
320 //G4cout << "BLIC: Secondary " <
321 // <<" "<< aNew->GetMomen
322 }
323 delete *iter;
324 }
325 delete cascaders;
326
327 #ifdef debug_BLIR_result
328 //G4cout << "Result analysis, secondaries "
329 //G4cout << "p_tot_raw " << p_raw << " sum p
330 G4double m_nucl= G4ParticleTable::GetPartic
331 GetIonMass(targetNucleus.GetZ_asInt(
332 // delete? tZ=targetNucleus.GetZ_asInt();
333
334 //G4cout << "BLIC Energy conservation initia
335 // << aTrack.GetTotalEnergy() + m_nuc
336 // <<" "<< aTrack.GetTotalEnergy() + m_
337 G4cout << "BLIC momentum conservation " << a
338 << " ptot " << ptot << " delta " << aTra
339 << " 3mom.mag() " << (aTrack.Get4
340 #endif
341
342 if(debug_G4BinaryLightIonReactionResults) G4
343
344 return &theResult;
345 }
346
347 //--------------------------------------------
348
349 //********************************************
350 G4bool G4BinaryLightIonReaction::EnergyAndMome
351 G4ReactionProductVector* Output, G4Lorentz
352 //********************************************
353 {
354 const int nAttemptScale = 2500;
355 const double ErrLimit = 1.E-6;
356 if (Output->empty())
357 return TRUE;
358 G4LorentzVector SumMom(0,0,0,0);
359 G4double SumMass = 0;
360 G4double TotalCollisionMass = TotalCo
361 size_t i = 0;
362 // Calculate sum hadron 4-momenta and summin
363 for(i = 0; i < Output->size(); i++)
364 {
365 SumMom += G4LorentzVector((*Output)[i]->G
366 SumMass += (*Output)[i]->GetDefinition()->
367 }
368 // G4cout << " E/P corrector, SumMass, Sum
369 // << SumMass <<" "<< SumMom.m2() <<
370 if (SumMass > TotalCollisionMass) return FAL
371 SumMass = SumMom.m2();
372 if (SumMass < 0) return FALSE;
373 SumMass = std::sqrt(SumMass);
374
375 // Compute c.m.s. hadron velocity and boost
376 G4ThreeVector Beta = -SumMom.boostVector();
377 //G4cout << " == pre boost 2 "<< SumMo
378 //--old Output->Boost(Beta);
379 for(i = 0; i < Output->size(); i++)
380 {
381 G4LorentzVector mom = G4LorentzVector((*Ou
382 mom *= Beta;
383 (*Output)[i]->SetMomentum(mom.vect());
384 (*Output)[i]->SetTotalEnergy(mom.e());
385 }
386
387 // Scale total c.m.s. hadron energy (hadron
388 // It should be equal interaction mass
389 G4double Scale = 0,OldScale=0;
390 G4double factor = 1.;
391 G4int cAttempt = 0;
392 G4double Sum = 0;
393 G4bool success = false;
394 for(cAttempt = 0; cAttempt < nAttemptScale;
395 {
396 Sum = 0;
397 for(i = 0; i < Output->size(); i++)
398 {
399 G4LorentzVector HadronMom = G4LorentzVec
400 HadronMom.setVect(HadronMom.vect()+ fact
401 G4double E = std::sqrt(HadronMom.vect().
402 HadronMom.setE(E);
403 (*Output)[i]->SetMomentum(HadronMom.vect
404 (*Output)[i]->SetTotalEnergy(HadronMom.e
405 Sum += E;
406 }
407 OldScale=Scale;
408 Scale = TotalCollisionMass/Sum - 1;
409 // G4cout << "E/P corr - " << cAttempt <<
410 if (std::abs(Scale) <= ErrLimit
411 || OldScale == Scale) // protect '
412 {
413 if (debug_G4BinaryLightIonReactionResult
414 success = true;
415 break;
416 }
417 if ( cAttempt > 10 )
418 {
419 // G4cout << " speed it up? " <<
420 factor=std::max(1.,G4Log(std::abs(OldSca
421 // G4cout << " ? factor ? " << factor
422 }
423 }
424
425 if( (!success) && debug_G4BinaryLightIonRea
426 {
427 G4cout << "G4G4BinaryLightIonReaction::Ene
428 G4cout << " Scale not unity at end of it
429 G4cout << " Increase number of attempts
430 }
431
432 // Compute c.m.s. interaction velocity and K
433 Beta = TotalCollisionMom.boostVector();
434 //--old Output->Boost(Beta);
435 for(i = 0; i < Output->size(); i++)
436 {
437 G4LorentzVector mom = G4LorentzVector((*Ou
438 mom *= Beta;
439 (*Output)[i]->SetMomentum(mom.vect());
440 (*Output)[i]->SetTotalEnergy(mom.e());
441 }
442 return TRUE;
443 }
444 G4bool G4BinaryLightIonReaction::SetLighterAsP
445 {
446 G4bool swapped = false;
447 if(tA<pA)
448 {
449 swapped = true;
450 G4int tmp(0);
451 tmp = tA; tA=pA; pA=tmp;
452 tmp = tZ; tZ=pZ; pZ=tmp;
453 G4double m1=G4ParticleTable::GetParticle
454 G4LorentzVector it(m1, G4ThreeVector(0,0
455 mom = toBreit*it;
456 }
457 return swapped;
458 }
459 G4ReactionProductVector * G4BinaryLightIonReac
460 {
461 // Check if kinematically nuclei can fuse.
462 G4double mFused=G4ParticleTable::GetParticl
463 G4double mTarget=G4ParticleTable::GetPartic
464 G4LorentzVector pCompound(mom.e()+mTarget,m
465 G4double m2Compound=pCompound.m2();
466 if (m2Compound < sqr(mFused) ) {
467 //G4cout << "G4BLIC: projectile p, mTarge
468 // << " " << sqrt(m2Compound)<< " "
469 return 0;
470 }
471
472 G4Fragment aPreFrag;
473 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA);
474 aPreFrag.SetNumberOfParticles(pA);
475 aPreFrag.SetNumberOfCharged(pZ);
476 aPreFrag.SetNumberOfHoles(0);
477 //GF FIXME: whyusing plop in z direction? t
478 //G4ThreeVector plop(0.,0., mom.vect().mag(
479 //G4LorentzVector aL(mom.t()+mTarget, plop)
480 G4LorentzVector aL(mom.t()+mTarget,mom.vect
481 aPreFrag.SetMomentum(aL);
482
483
484 //G4cout << "Fragment INFO "<< pA+tA
485 // << aL <<" "<<G4endl << aPreF
486 G4ReactionProductVector * cascaders = thePr
487 //G4double tSum = 0;
488 for(size_t count = 0; count<cascaders->size
489 {
490 cascaders->operator[](count)->SetNewlyAd
491 //tSum += cascaders->operator[](count)->
492 }
493 // G4cout << "Exiting pre-compound on
494 return cascaders;
495 }
496 G4ReactionProductVector * G4BinaryLightIonReac
497 {
498 G4ReactionProductVector * result = 0;
499 G4double projectileMass(0);
500 G4LorentzVector it;
501
502 G4int tryCount(0);
503 do
504 {
505 ++tryCount;
506 projectile3dNucleus = new G4Fancy3DNu
507 projectile3dNucleus->Init(pA, pZ);
508 projectile3dNucleus->CenterNucleons()
509 projectileMass=G4ParticleTable::GetPa
510 projectile3dNucleus->GetCharge(
511 it=toBreit * G4LorentzVector(projecti
512
513 target3dNucleus = new G4Fancy3DNucleu
514 target3dNucleus->Init(tA, tZ);
515 G4double impactMax = target3dNucleus-
516 // G4cout << "out radius - nuc
517 G4double aX=(2.*G4UniformRand()-1.)*i
518 G4double aY=(2.*G4UniformRand()-1.)*i
519 G4ThreeVector pos(aX, aY, -2.*impactM
520
521 G4KineticTrackVector * initalState =
522 projectile3dNucleus->StartLoop();
523 G4Nucleon * aNuc;
524 G4LorentzVector tmpV(0,0,0,0);
525 #ifdef debug_BLIR_finalstate
526 G4LorentzVector pinitial;
527 #endif
528 G4LorentzVector nucleonMom(1./pA*mom)
529 nucleonMom.setZ(nucleonMom.vect().mag
530 nucleonMom.setX(0);
531 nucleonMom.setY(0);
532 theFermi.Init(pA,pZ);
533 while( (aNuc=projectile3dNucleus->Get
534 {
535 G4LorentzVector p4 = aNuc->GetMome
536 tmpV+=p4;
537 G4ThreeVector nucleonPosition(aNuc
538 G4double density=(projectile3dNucl
539 nucleonPosition += pos;
540 G4KineticTrack * it1 = new G4Kinet
541 it1->SetState(G4KineticTrack::outs
542 G4double pfermi= theFermi.GetFermi
543 G4double mass = aNuc->GetDefinitio
544 G4double Efermi= std::sqrt( sqr(ma
545 it1->SetProjectilePotential(-Eferm
546 initalState->push_back(it1);
547 #ifdef debug_BLIR_finalstate
548 pinitial += it1->Get4Momentum()
549 #endif
550 }
551
552 result=theModel->Propagate(initalStat
553 #ifdef debug_BLIR_finalstate
554 if( result && result->size()>0)
555 {
556 G4cout << " Cascade result " << G4endl;
557 G4LorentzVector presult;
558 G4ReactionProductVector::iterator
559 G4ReactionProduct xp;
560 for (iter=result->begin(); iter !
561 {
562 presult += G4LorentzVector((*ite
563 G4cout << (*iter)->GetDefinition()->Ge
564 << "("<< (*iter)->GetMomentum().x()<<"
565 << (*iter)->GetMomentum().y()<<","
566 << (*iter)->GetMomentum().z()<<";"
567 << (*iter)->GetTotalEnergy() <<")"<
568 }
569
570 G4cout << "BLIC check result : in
571 << " final " << presult
572 << " IF - FF " << pinitial +G
573
574 }
575 #endif
576 if( result && result->size()==0)
577 {
578 delete result;
579 result=0;
580 }
581 if ( ! result )
582 {
583 delete target3dNucleus;
584 delete projectile3dNucleus;
585 }
586
587 // std::for_each(initalState->begin()
588 // delete initalState;
589
590 } while (! result && tryCount< 150); /
591 return result;
592 }
593 G4double G4BinaryLightIonReaction::GetProjecti
594 {
595
596 G4Nucleon * aNuc;
597 // the projectileNucleus excitation ener
598 G4double theStatisticalExEnergy = 0;
599 projectile3dNucleus->StartLoop();
600 while( (aNuc=projectile3dNucleus->GetNex
601 {
602 //G4cout << " Nucleon : " << a
603 if(aNuc->AreYouHit()) {
604 G4ThreeVector aPosition(aNuc->GetP
605 G4double localDensity = projectile
606 G4double localPfermi = theFermi.Ge
607 G4double nucMass = aNuc->GetDefini
608 G4double localFermiEnergy = std::s
609 G4double deltaE = localFermiEnergy
610 theStatisticalExEnergy += deltaE;
611 }
612 }
613 return theStatisticalExEnergy;
614 }
615
616 G4LorentzVector G4BinaryLightIonReaction::Sort
617 {
618 unsigned int i(0);
619 spectatorA=spectatorZ=0;
620 G4LorentzVector pspectators(0,0,0,0);
621 pFinalState=G4LorentzVector(0,0,0,0);
622 for(i=0; i<result->size(); i++)
623 {
624 if( (*result)[i]->GetNewlyAdded() )
625 {
626 pFinalState += G4LorentzVector( (*res
627 cascaders->push_back((*result)[i]);
628 }
629 else {
630 // G4cout <<" spectator ...
631 pspectators += G4LorentzVector( (*res
632 spectators->push_back((*result)[i]);
633 spectatorA++;
634 spectatorZ+= G4lrint((*result)[i]->Ge
635 }
636
637 // G4cout << (*result)[i]<< " "
638 // << (*result)[i]->GetDefinition
639 // << (*result)[i]->GetMomentum()
640 // << (*result)[i]->GetTotalEnerg
641 }
642 //G4cout << "pFinalState / pspectators,
643 // << " (" << spectatorA << ", "<< sp
644
645 return pspectators;
646 }
647
648 void G4BinaryLightIonReaction::DeExciteSpectat
649
650 {
651 // call precompound model
652 G4ReactionProductVector * proFrag = 0;
653 G4LorentzVector pFragment(0.,0.,0.,0.);
654 // G4cout << " == pre boost 1 "<< mome
655 G4LorentzRotation boost_fragments;
656 // G4cout << " == post boost 1 "<< mom
657 // G4LorentzRotation boost_spectator_mom
658 // G4cout << "- momentum " << boost_spe
659 G4LorentzVector pFragments(0,0,0,0);
660
661 if(spectatorZ>0 && spectatorA>1)
662 {
663 // Make the fragment
664 G4Fragment aProRes;
665 aProRes.SetZandA_asInt(spectatorZ, spect
666 aProRes.SetNumberOfParticles(0);
667 aProRes.SetNumberOfCharged(0);
668 aProRes.SetNumberOfHoles(pA-spectatorA);
669 G4double mFragment=G4ParticleTable::GetP
670 pFragment=G4LorentzVector(0,0,0,mFragmen
671 aProRes.SetMomentum(pFragment);
672
673 proFrag = theHandler->BreakItUp(aProRes)
674
675 boost_fragments = G4LorentzRotation(pSpe
676
677 // G4cout << " Fragment a,z, Mass Fr
678 // << spectatorA <<" "<< spectator
679 // << momentum.mag() <<" "<< momen
680 // << " "<<theStatisticalExEnergy
681 // << " "<< boost_fragments*pFragm
682 G4ReactionProductVector::iterator ispect
683 for (ispectator=spectators->begin();ispe
684 {
685 delete *ispectator;
686 }
687 }
688 else if(spectatorA!=0)
689 {
690 G4ReactionProductVector::iterator ispecta
691 for (ispectator=spectators->begin();ispec
692 {
693 (*ispectator)->SetNewlyAdded(true);
694 cascaders->push_back(*ispectator);
695 pFinalState+=G4LorentzVector((*ispect
696 //G4cout << "BLIC: spectator
697 // << (*ispectator)->GetDefi
698 // << (*ispectator)->GetMome
699 // << (*ispectator)->GetTota
700 }
701
702 }
703 // / if (spectators)
704 delete spectators;
705
706 // collect the evaporation part and boost t
707 G4ReactionProductVector::iterator ii;
708 if(proFrag)
709 {
710 for(ii=proFrag->begin(); ii!=proFrag->en
711 {
712 (*ii)->SetNewlyAdded(true);
713 G4LorentzVector tmp((*ii)->GetMomentu
714 tmp *= boost_fragments;
715 (*ii)->SetMomentum(tmp.vect());
716 (*ii)->SetTotalEnergy(tmp.e());
717 // result->push_back(*ii);
718 pFragments += tmp;
719 }
720 }
721
722 // G4cout << "Fragmented p, momentum, de
723 // <<" "<< pFragments-momentum <
724
725 // correct p/E of Cascade secondaries
726 G4LorentzVector pCas=pInitialState - pFragm
727
728 //G4cout <<"BLIC: Going to correct from
729 // the creation of excited fragment did vi
730 G4bool EnergyIsCorrect=EnergyAndMomentumCor
731 if ( ! EnergyIsCorrect && debug_G4BinaryLig
732 {
733 G4cout << "G4BinaryLightIonReaction E/P
734 }
735
736 // Add deexcitation secondaries
737 if(proFrag)
738 {
739 for(ii=proFrag->begin(); ii!=proFrag->en
740 {
741 cascaders->push_back(*ii);
742 }
743 delete proFrag;
744 }
745 //G4cout << "EnergyIsCorrect? " << Energ
746 if ( ! EnergyIsCorrect )
747 {
748 // G4cout <<" ! EnergyIsCorrect " <<
749 if (! EnergyAndMomentumCorrector(cascade
750 {
751 if(debug_G4BinaryLightIonReactionResu
752 G4cout << "G4BinaryLightIonReactio
753 }
754 }
755
756 }
757
758