| Geant4 Cross Reference |
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 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (May 1998)
28 //
29 //
30 // Modified:
31 // 30 June 1998 by V. Lara:
32 // -Modified the Transform method for use G4ParticleTable and
33 // therefore G4IonTable. It makes possible to convert all kind
34 // of fragments (G4Fragment) produced in deexcitation to
35 // G4DynamicParticle
36 // -It uses default algorithms for:
37 // Evaporation: G4Evaporation
38 // MultiFragmentation: G4StatMF
39 // Fermi Breakup model: G4FermiBreakUp
40 // 24 Jul 2008 by M. A. Cortes Giraldo:
41 // -Max Z,A for Fermi Break-Up turns to 9,17 by default
42 // -BreakItUp() reorganised and bug in Evaporation loop fixed
43 // -Transform() optimised
44 // (September 2008) by J. M. Quesada. External choices have been added for :
45 // -inverse cross section option (default OPTxs=3)
46 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
47 // September 2009 by J. M. Quesada:
48 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
49 // 27 Nov 2009 by V.Ivanchenko:
50 // -cleanup the logic, reduce number internal vectors, fixed memory leak.
51 // 11 May 2010 by V.Ivanchenko:
52 // -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
53 // final photon deexcitation; used check on adundance of a fragment, decay
54 // unstable fragments with A <5
55 // 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition:
56 // products of Fermi Break Up cannot be further deexcited by this model
57 // 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods
58 // to the source
59 // 23 January 2012 by V.Ivanchenko general cleanup including destruction of
60 // objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and
61 // not delete it here
62
63 #include "G4ExcitationHandler.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4LorentzVector.hh"
66 #include "G4ThreeVector.hh"
67 #include "G4ParticleTable.hh"
68 #include "G4ParticleTypes.hh"
69 #include "G4Ions.hh"
70 #include "G4Electron.hh"
71 #include "G4Lambda.hh"
72
73 #include "G4VMultiFragmentation.hh"
74 #include "G4VFermiBreakUp.hh"
75 #include "G4Element.hh"
76 #include "G4ElementTable.hh"
77
78 #include "G4VEvaporation.hh"
79 #include "G4VEvaporationChannel.hh"
80 #include "G4Evaporation.hh"
81 #include "G4PhotonEvaporation.hh"
82 #include "G4StatMF.hh"
83 #include "G4FermiBreakUpVI.hh"
84 #include "G4NuclearLevelData.hh"
85 #include "G4PhysicsModelCatalog.hh"
86
87 G4ExcitationHandler::G4ExcitationHandler()
88 : minEForMultiFrag(1.*CLHEP::TeV), minExcitation(1.*CLHEP::eV),
89 maxExcitation(100.*CLHEP::MeV)
90 {
91 thePartTable = G4ParticleTable::GetParticleTable();
92 theTableOfIons = thePartTable->GetIonTable();
93 nist = G4NistManager::Instance();
94
95 theMultiFragmentation = new G4StatMF();
96 theFermiModel = new G4FermiBreakUpVI();
97 thePhotonEvaporation = new G4PhotonEvaporation();
98 SetEvaporation(new G4Evaporation(thePhotonEvaporation), true);
99 theResults.reserve(60);
100 results.reserve(30);
101 theEvapList.reserve(30);
102
103 theElectron = G4Electron::Electron();
104 theNeutron = G4Neutron::NeutronDefinition();
105 theProton = G4Proton::ProtonDefinition();
106 theDeuteron = G4Deuteron::DeuteronDefinition();
107 theTriton = G4Triton::TritonDefinition();
108 theHe3 = G4He3::He3Definition();
109 theAlpha = G4Alpha::AlphaDefinition();
110 theLambda = G4Lambda::Lambda();
111
112 fLambdaMass = theLambda->GetPDGMass();
113
114 if(fVerbose > 1) { G4cout << "### New handler " << this << G4endl; }
115 }
116
117 G4ExcitationHandler::~G4ExcitationHandler()
118 {
119 delete theMultiFragmentation;
120 delete theFermiModel;
121 if(isEvapLocal) { delete theEvaporation; }
122 }
123
124 void G4ExcitationHandler::SetParameters()
125 {
126 G4NuclearLevelData* ndata = G4NuclearLevelData::GetInstance();
127 auto param = ndata->GetParameters();
128 isActive = true;
129 // check if de-excitation is needed
130 if (fDummy == param->GetDeexChannelsType()) {
131 isActive = false;
132 } else {
133 // upload data for elements used in geometry
134 G4int Zmax = 20;
135 const G4ElementTable* table = G4Element::GetElementTable();
136 for (auto const & elm : *table) { Zmax = std::max(Zmax, elm->GetZasInt()); }
137 ndata->UploadNuclearLevelData(Zmax+1);
138 }
139 minEForMultiFrag = param->GetMinExPerNucleounForMF();
140 minExcitation = param->GetMinExcitation();
141 maxExcitation = param->GetPrecoHighEnergy();
142
143 // allowing local debug printout
144 fVerbose = std::max(fVerbose, param->GetVerbose());
145 if (isActive) {
146 if (nullptr == thePhotonEvaporation) {
147 SetPhotonEvaporation(new G4PhotonEvaporation());
148 }
149 if (nullptr == theFermiModel) {
150 SetFermiModel(new G4FermiBreakUpVI());
151 }
152 if (nullptr == theMultiFragmentation) {
153 SetMultiFragmentation(new G4StatMF());
154 }
155 if (nullptr == theEvaporation) {
156 SetEvaporation(new G4Evaporation(thePhotonEvaporation), true);
157 }
158 }
159 theFermiModel->SetVerbose(fVerbose);
160 if(fVerbose > 1) {
161 G4cout << "G4ExcitationHandler::SetParameters() done " << this << G4endl;
162 }
163 }
164
165 void G4ExcitationHandler::Initialise()
166 {
167 if(isInitialised) { return; }
168 if(fVerbose > 1) {
169 G4cout << "G4ExcitationHandler::Initialise() started " << this << G4endl;
170 }
171 G4DeexPrecoParameters* param =
172 G4NuclearLevelData::GetInstance()->GetParameters();
173 isInitialised = true;
174 SetParameters();
175 if(isActive) {
176 theFermiModel->Initialise();
177 theEvaporation->InitialiseChannels();
178 }
179 // dump level is controlled by parameter class
180 param->Dump();
181 }
182
183 void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr, G4bool flag)
184 {
185 if(nullptr != ptr && ptr != theEvaporation) {
186 theEvaporation = ptr;
187 theEvaporation->SetPhotonEvaporation(thePhotonEvaporation);
188 theEvaporation->SetFermiBreakUp(theFermiModel);
189 isEvapLocal = flag;
190 if(fVerbose > 1) {
191 G4cout << "G4ExcitationHandler::SetEvaporation() " << ptr << " done for " << this << G4endl;
192 }
193 }
194 }
195
196 void
197 G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr)
198 {
199 if(nullptr != ptr && ptr != theMultiFragmentation) {
200 delete theMultiFragmentation;
201 theMultiFragmentation = ptr;
202 }
203 }
204
205 void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr)
206 {
207 if(nullptr != ptr && ptr != theFermiModel) {
208 delete theFermiModel;
209 theFermiModel = ptr;
210 if(nullptr != theEvaporation) {
211 theEvaporation->SetFermiBreakUp(theFermiModel);
212 }
213 }
214 }
215
216 void
217 G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
218 {
219 if(nullptr != ptr && ptr != thePhotonEvaporation) {
220 delete thePhotonEvaporation;
221 thePhotonEvaporation = ptr;
222 if(nullptr != theEvaporation) {
223 theEvaporation->SetPhotonEvaporation(ptr);
224 }
225 if(fVerbose > 1) {
226 G4cout << "G4ExcitationHandler::SetPhotonEvaporation() " << ptr
227 << " for handler " << this << G4endl;
228 }
229 }
230 }
231
232 void G4ExcitationHandler::SetDeexChannelsType(G4DeexChannelType val)
233 {
234 G4Evaporation* evap = static_cast<G4Evaporation*>(theEvaporation);
235 if(fVerbose > 1) {
236 G4cout << "G4ExcitationHandler::SetDeexChannelsType " << val
237 << " for " << this << G4endl;
238 }
239 if(val == fDummy) {
240 isActive = false;
241 return;
242 }
243 if (nullptr == evap) { return; }
244 if (val == fEvaporation) {
245 evap->SetDefaultChannel();
246 } else if (val == fCombined) {
247 evap->SetCombinedChannel();
248 } else if (val == fGEM) {
249 evap->SetGEMChannel();
250 } else if (val == fGEMVI) {
251 evap->SetGEMVIChannel();
252 }
253 evap->InitialiseChannels();
254 if (fVerbose > 1) {
255 if (G4Threading::IsMasterThread()) {
256 G4cout << "Number of de-excitation channels is changed to: "
257 << theEvaporation->GetNumberOfChannels();
258 G4cout << " " << this;
259 }
260 G4cout << G4endl;
261 }
262 }
263
264 G4VEvaporation* G4ExcitationHandler::GetEvaporation()
265 {
266 if (nullptr != theEvaporation) { SetParameters(); }
267 return theEvaporation;
268 }
269
270 G4VMultiFragmentation* G4ExcitationHandler::GetMultiFragmentation()
271 {
272 if (nullptr != theMultiFragmentation) { SetParameters(); }
273 return theMultiFragmentation;
274 }
275
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermiModel()
277 {
278 if (nullptr != theFermiModel) { SetParameters(); }
279 return theFermiModel;
280 }
281
282 G4VEvaporationChannel* G4ExcitationHandler::GetPhotonEvaporation()
283 {
284 if(nullptr != thePhotonEvaporation) { SetParameters(); }
285 return thePhotonEvaporation;
286 }
287
288 G4ReactionProductVector *
289 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState)
290 {
291 // Variables existing until end of method
292 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
293 if (fVerbose > 1) {
294 G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " << G4endl;
295 G4cout << theInitialState << G4endl;
296 }
297 if (!isInitialised) { Initialise(); }
298
299 // pointer to fragment vector which receives temporal results
300 G4FragmentVector * theTempResult = nullptr;
301
302 theResults.clear();
303 theEvapList.clear();
304
305 // Variables to describe the excited configuration
306 G4double exEnergy = theInitialState.GetExcitationEnergy();
307 G4int A = theInitialState.GetA_asInt();
308 G4int Z = theInitialState.GetZ_asInt();
309 G4int nL = theInitialState.GetNumberOfLambdas();
310
311 // too much excitation
312 if (exEnergy > A*maxExcitation && A > 0) {
313 ++fWarnings;
314 if(fWarnings < 0) {
315 G4ExceptionDescription ed;
316 ed << "High excitation Fragment Z= " << Z << " A= " << A
317 << " Eex/A(MeV)= " << exEnergy/A;
318 G4Exception("G4ExcitationHandler::BreakItUp()","had0034",JustWarning,ed,"");
319 }
320 }
321
322 // for hyper-nuclei subtract lambdas from the projectile fragment
323 G4double lambdaF = 0.0;
324 G4LorentzVector lambdaLV = theInitialStatePtr->GetMomentum();
325 if (0 < nL) {
326
327 // is it a stable hyper-nuclei?
328 if(A >= 3 && A <= 5 && nL <= 2) {
329 G4int pdg = 0;
330 if(3 == A && 1 == nL) {
331 pdg = 1010010030;
332 } else if(5 == A && 2 == Z && 1 == nL) {
333 pdg = 1010020050;
334 } else if(4 == A) {
335 if(1 == Z && 1 == nL) {
336 pdg = 1010010040;
337 } else if(2 == Z && 1 == nL) {
338 pdg = 1010020040;
339 } else if(0 == Z && 2 == nL) {
340 pdg = 1020000040;
341 } else if(1 == Z && 2 == nL) {
342 pdg = 1020010040;
343 }
344 }
345 // initial state is one of hyper-nuclei
346 if (0 < pdg) {
347 const G4ParticleDefinition* part = thePartTable->FindParticle(pdg);
348 if(nullptr != part) {
349 G4ReactionProduct* theNew = new G4ReactionProduct(part);
350 G4ThreeVector dir = G4ThreeVector( 0.0, 0.0, 0.0 );
351 if ( lambdaLV.vect().mag() > CLHEP::eV ) {
352 dir = lambdaLV.vect().unit();
353 }
354 G4double mass = part->GetPDGMass();
355 G4double etot = std::max(lambdaLV.e(), mass);
356 dir *= std::sqrt((etot - mass)*(etot + mass));
357 theNew->SetMomentum(dir);
358 theNew->SetTotalEnergy(etot);
359 theNew->SetFormationTime(theInitialState.GetCreationTime());
360 theNew->SetCreatorModelID(theInitialState.GetCreatorModelID());
361 G4ReactionProductVector* v = new G4ReactionProductVector();
362 v->push_back(theNew);
363 return v;
364 }
365 }
366 }
367 G4double mass = theInitialStatePtr->GetGroundStateMass();
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron_mass_c2)/mass;
369
370 // de-excitation with neutrons instead of lambda inside the fragment
371 theInitialStatePtr->SetZAandMomentum(lambdaLV*(1. - lambdaF), Z, A, 0);
372
373 // 4-momentum not used in de-excitation
374 lambdaLV *= lambdaF;
375 } else if (0 > nL) {
376 ++fWarnings;
377 if(fWarnings < 0) {
378 G4ExceptionDescription ed;
379 ed << "Fragment with negative L: Z=" << Z << " A=" << A << " L=" << nL
380 << " Eex/A(MeV)= " << exEnergy/A;
381 G4Exception("G4ExcitationHandler::BreakItUp()","had0034",JustWarning,ed,"");
382 }
383 }
384
385 // In case A <= 1 the fragment will not perform any nucleon emission
386 if (A <= 1 || !isActive) {
387 theResults.push_back( theInitialStatePtr );
388
389 // check if a fragment is stable
390 } else if (exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) {
391 theResults.push_back( theInitialStatePtr );
392
393 // JMQ 150909: first step in de-excitation is treated separately
394 // Fragments after the first step are stored in theEvapList
395 } else {
396 if ((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
397 || exEnergy <= minEForMultiFrag*A) {
398 theEvapList.push_back(theInitialStatePtr);
399
400 // Statistical Multifragmentation will take place only once
401 } else {
402 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
403 if (nullptr == theTempResult) {
404 theEvapList.push_back(theInitialStatePtr);
405 } else {
406 std::size_t nsec = theTempResult->size();
407
408 // no fragmentation
409 if (0 == nsec) {
410 theEvapList.push_back(theInitialStatePtr);
411
412 // secondary are produced - sort out secondary fragments
413 } else {
414 G4bool deletePrimary = true;
415 for (auto const & ptr : *theTempResult) {
416 if (ptr == theInitialStatePtr) { deletePrimary = false; }
417 SortSecondaryFragment(ptr);
418 }
419 if (deletePrimary) { delete theInitialStatePtr; }
420 }
421 delete theTempResult; // end multifragmentation
422 }
423 }
424 }
425 if (fVerbose > 2) {
426 G4cout << "## After first step of handler " << theEvapList.size()
427 << " for evap; "
428 << theResults.size() << " results. " << G4endl;
429 }
430 // -----------------------------------
431 // FermiBreakUp and De-excitation loop
432 // -----------------------------------
433
434 static const G4int countmax = 1000;
435 std::size_t kk;
436 for (kk=0; kk<theEvapList.size(); ++kk) {
437 G4Fragment* frag = theEvapList[kk];
438 if (fVerbose > 3) {
439 G4cout << "Next evaporate: " << G4endl;
440 G4cout << *frag << G4endl;
441 }
442 if (kk >= countmax) {
443 G4ExceptionDescription ed;
444 ed << "Infinite loop in the de-excitation module: " << kk
445 << " iterations \n"
446 << " Initial fragment: \n" << theInitialState
447 << "\n Current fragment: \n" << *frag;
448 G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException,
449 ed,"Stop execution");
450
451 }
452 A = frag->GetA_asInt();
453 Z = frag->GetZ_asInt();
454 results.clear();
455 if (fVerbose > 2) {
456 G4cout << "G4ExcitationHandler# " << kk << " Z= " << Z << " A= " << A
457 << " Eex(MeV)= " << frag->GetExcitationEnergy() << G4endl;
458 }
459 // Fermi Break-Up
460 if (theFermiModel->IsApplicable(Z, A, frag->GetExcitationEnergy())) {
461 theFermiModel->BreakFragment(&results, frag);
462 std::size_t nsec = results.size();
463 if (fVerbose > 2) { G4cout << "FermiBreakUp Nsec= " << nsec << G4endl; }
464
465 // FBU takes care to delete input fragment or add it to the results
466 // The secondary may be excited - photo-evaporation should be applied
467 if (1 < nsec) {
468 for (auto const & res : results) {
469 SortSecondaryFragment(res);
470 }
471 continue;
472 }
473 // evaporation will be applied
474 }
475 // apply Evaporation, residual nucleus is always added to the results
476 // photon evaporation is possible
477 theEvaporation->BreakFragment(&results, frag);
478 if (fVerbose > 3) {
479 G4cout << "Evaporation Nsec= " << results.size() << G4endl;
480 }
481 if (0 == results.size()) {
482 theResults.push_back(frag);
483 } else {
484 SortSecondaryFragment(frag);
485 }
486
487 // Sort out secondary fragments
488 for (auto const & res : results) {
489 if(fVerbose > 4) {
490 G4cout << "Evaporated product #" << *res << G4endl;
491 }
492 SortSecondaryFragment(res);
493 } // end of loop on secondary
494 } // end of the loop over theEvapList
495 if (fVerbose > 2) {
496 G4cout << "## After 2nd step of handler " << theEvapList.size()
497 << " was evap; "
498 << theResults.size() << " results. " << G4endl;
499 }
500 G4ReactionProductVector * theReactionProductVector =
501 new G4ReactionProductVector();
502
503 // MAC (24/07/08)
504 // To optimise the storing speed, we reserve space
505 // in memory for the vector
506 theReactionProductVector->reserve( theResults.size() );
507
508 if (fVerbose > 1) {
509 G4cout << "### ExcitationHandler provides " << theResults.size()
510 << " evaporated products:" << G4endl;
511 }
512 G4LorentzVector partOfLambdaLV;
513 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4double)nL;
514 for (auto const & frag : theResults) {
515 G4LorentzVector lv0 = frag->GetMomentum();
516 G4double etot = lv0.e();
517
518 // in the case of dummy de-excitation, excitation energy is transfered
519 // into kinetic energy of output ion
520 if (!isActive) {
521 G4double mass = frag->GetGroundStateMass();
522 G4double ptot = lv0.vect().mag();
523 G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0
524 : std::sqrt((etot - mass)*(etot + mass))/ptot;
525 G4LorentzVector lv((frag->GetMomentum()).px()*fac,
526 (frag->GetMomentum()).py()*fac,
527 (frag->GetMomentum()).pz()*fac, etot);
528 frag->SetMomentum(lv);
529 }
530 if (fVerbose > 3) {
531 G4cout << *frag;
532 if (frag->NuclearPolarization()) {
533 G4cout << " " << frag->NuclearPolarization();
534 }
535 G4cout << G4endl;
536 }
537
538 G4int fragmentA = frag->GetA_asInt();
539 G4int fragmentZ = frag->GetZ_asInt();
540 G4double eexc = 0.0;
541 const G4ParticleDefinition* theKindOfFragment = nullptr;
542 G4bool isHyperN = false;
543 if (fragmentA == 0) { // photon or e-
544 theKindOfFragment = frag->GetParticleDefinition();
545 } else if (fragmentA == 1 && fragmentZ == 0) { // neutron
546 theKindOfFragment = theNeutron;
547 } else if (fragmentA == 1 && fragmentZ == 1) { // proton
548 theKindOfFragment = theProton;
549 } else if (fragmentA == 2 && fragmentZ == 1) { // deuteron
550 theKindOfFragment = theDeuteron;
551 } else if (fragmentA == 3 && fragmentZ == 1) { // triton
552 theKindOfFragment = theTriton;
553 if(0 < nL) {
554 const G4ParticleDefinition* p = thePartTable->FindParticle(1010010030);
555 if(nullptr != p) {
556 theKindOfFragment = p;
557 isHyperN = true;
558 --nL;
559 }
560 }
561 } else if (fragmentA == 3 && fragmentZ == 2) { // helium3
562 theKindOfFragment = theHe3;
563 } else if (fragmentA == 4 && fragmentZ == 2) { // alpha
564 theKindOfFragment = theAlpha;
565 if (0 < nL) {
566 const G4ParticleDefinition* p = thePartTable->FindParticle(1010020040);
567 if(nullptr != p) {
568 theKindOfFragment = p;
569 isHyperN = true;
570 --nL;
571 }
572 }
573 } else {
574
575 // fragment
576 eexc = frag->GetExcitationEnergy();
577 G4int idxf = frag->GetFloatingLevelNumber();
578 if (eexc < minExcitation) {
579 eexc = 0.0;
580 idxf = 0;
581 }
582
583 theKindOfFragment = theTableOfIons->GetIon(fragmentZ, fragmentA, eexc,
584 G4Ions::FloatLevelBase(idxf));
585 if (fVerbose > 3) {
586 G4cout << "### EXCH: Find ion Z= " << fragmentZ
587 << " A= " << fragmentA
588 << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf
589 << " " << theKindOfFragment->GetParticleName()
590 << G4endl;
591 }
592 }
593 // fragment identified
594 if (nullptr != theKindOfFragment) {
595 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
596 if (isHyperN) {
597 G4LorentzVector lv = lv0 + partOfLambdaLV;
598 G4ThreeVector dir = lv.vect().unit();
599 G4double mass = theKindOfFragment->GetPDGMass();
600 etot = std::max(lv.e(), mass);
601 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
602 dir *= ptot;
603 theNew->SetMomentum(dir);
604 // remaining not compensated 4-momentum
605 lambdaLV += (lv0 - G4LorentzVector(dir, etot));
606 } else {
607 theNew->SetMomentum(lv0.vect());
608 }
609 theNew->SetTotalEnergy(etot);
610 theNew->SetFormationTime(frag->GetCreationTime());
611 theNew->SetCreatorModelID(frag->GetCreatorModelID());
612 theReactionProductVector->push_back(theNew);
613
614 // fragment not found out ground state is created
615 } else {
616 theKindOfFragment =
617 theTableOfIons->GetIon(fragmentZ,fragmentA,0.0,noFloat,0);
618 if (theKindOfFragment) {
619 G4ThreeVector mom(0.0,0.0,0.0);
620 G4double ionmass = theKindOfFragment->GetPDGMass();
621 if (etot <= ionmass) {
622 etot = ionmass;
623 } else {
624 G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
625 mom = (frag->GetMomentum().vect().unit())*ptot;
626 }
627 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
628 theNew->SetMomentum(mom);
629 theNew->SetTotalEnergy(etot);
630 theNew->SetFormationTime(frag->GetCreationTime());
631 theNew->SetCreatorModelID(frag->GetCreatorModelID());
632 theReactionProductVector->push_back(theNew);
633 if (fVerbose > 3) {
634 G4cout << " ground state, energy corrected E(MeV)= "
635 << etot << G4endl;
636 }
637 }
638 }
639 delete frag;
640 }
641 // remaining lambdas are free; conserve quantum numbers but
642 // not 4-momentum
643 if (0 < nL) {
644 G4ThreeVector dir = G4ThreeVector(0.0, 0.0, 0.0);
645 if (lambdaLV.vect().mag() > CLHEP::eV) {
646 dir = lambdaLV.vect().unit();
647 }
648 G4double etot = std::max(lambdaLV.e()/(G4double)nL, fLambdaMass);
649 dir *= std::sqrt((etot - fLambdaMass)*(etot + fLambdaMass));
650 for (G4int i=0; i<nL; ++i) {
651 G4ReactionProduct* theNew = new G4ReactionProduct(theLambda);
652 theNew->SetMomentum(dir);
653 theNew->SetTotalEnergy(etot);
654 theNew->SetFormationTime(theInitialState.GetCreationTime());
655 theNew->SetCreatorModelID(theInitialState.GetCreatorModelID());
656 theReactionProductVector->push_back(theNew);
657 }
658 }
659 if (fVerbose > 3) {
660 G4cout << "@@@@@@@@@@ End G4Excitation Handler "<< G4endl;
661 }
662 return theReactionProductVector;
663 }
664
665 void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
666 {
667 outFile << "G4ExcitationHandler description\n"
668 << "This class samples de-excitation of excited nucleus using\n"
669 << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
670 << "evaporation, fission, and photo-evaporation models. Evaporated\n"
671 << "particle may be proton, neutron, and other light fragment \n"
672 << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
673 << "or electrons due to internal conversion \n";
674 }
675
676
677
678