1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | |
39 | |
40 | |
41 | |
42 | |
43 | |
44 | |
45 | |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | |
52 | |
53 | |
54 | |
55 | |
56 | |
57 | |
58 | |
59 | |
60 | |
61 | |
62 | |
63 | |
64 | |
65 | #include "G4ExcitationHandler.hh" |
66 | #include "G4SystemOfUnits.hh" |
67 | #include "G4LorentzVector.hh" |
68 | #include "G4NistManager.hh" |
69 | #include "G4ParticleTable.hh" |
70 | #include "G4ParticleTypes.hh" |
71 | #include "G4Ions.hh" |
72 | #include "G4Electron.hh" |
73 | |
74 | #include "G4VMultiFragmentation.hh" |
75 | #include "G4VFermiBreakUp.hh" |
76 | |
77 | #include "G4VEvaporation.hh" |
78 | #include "G4VEvaporationChannel.hh" |
79 | #include "G4Evaporation.hh" |
80 | #include "G4PhotonEvaporation.hh" |
81 | #include "G4StatMF.hh" |
82 | #include "G4FermiBreakUpVI.hh" |
83 | #include "G4NuclearLevelData.hh" |
84 | #include "G4Pow.hh" |
85 | |
86 | G4ExcitationHandler::G4ExcitationHandler() |
87 | : maxZForFermiBreakUp(9),maxAForFermiBreakUp(17), |
88 | fVerbose(0),isInitialised(false),isEvapLocal(true) |
89 | { |
90 | theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable(); |
91 | nist = G4NistManager::Instance(); |
92 | |
93 | theMultiFragmentation = nullptr; |
94 | theFermiModel = nullptr; |
95 | G4Pow::GetInstance(); |
96 | theEvaporation = new G4Evaporation(); |
97 | thePhotonEvaporation = theEvaporation->GetPhotonEvaporation(); |
98 | theResults.reserve(60); |
99 | results.reserve(30); |
100 | theEvapList.reserve(30); |
101 | thePhotoEvapList.reserve(10); |
102 | |
103 | SetParameters(); |
104 | electron = G4Electron::Electron(); |
105 | |
106 | if(fVerbose > 0) { G4cout(*G4cout_p) << "### New handler " << this << G4endlstd::endl; } |
107 | } |
108 | |
109 | G4ExcitationHandler::~G4ExcitationHandler() |
110 | { |
111 | |
112 | delete theMultiFragmentation; |
113 | delete theFermiModel; |
114 | if(isEvapLocal) { delete theEvaporation; } |
115 | } |
116 | |
117 | void G4ExcitationHandler::SetParameters() |
118 | { |
119 | G4DeexPrecoParameters* param = |
120 | G4NuclearLevelData::GetInstance()->GetParameters(); |
121 | isActive = true; |
122 | if(fDummy == param->GetDeexChannelsType()) { isActive = false; } |
123 | minEForMultiFrag = param->GetMinExPerNucleounForMF(); |
124 | minExcitation = param->GetMinExcitation(); |
125 | icID = param->GetInternalConversionID(); |
126 | if(isActive) { |
127 | if(!thePhotonEvaporation) { |
128 | SetPhotonEvaporation(new G4PhotonEvaporation()); |
129 | } |
130 | if(!theFermiModel) { SetFermiModel(new G4FermiBreakUpVI()); } |
131 | if(!theMultiFragmentation) { SetMultiFragmentation(new G4StatMF()); } |
132 | } |
133 | } |
134 | |
135 | void G4ExcitationHandler::Initialise() |
136 | { |
137 | if(isInitialised) { return; } |
138 | if(fVerbose > 0) { |
139 | G4cout(*G4cout_p) << "G4ExcitationHandler::Initialise() started " << this << G4endlstd::endl; |
140 | } |
141 | G4DeexPrecoParameters* param = |
| Value stored to 'param' during its initialization is never read |
142 | G4NuclearLevelData::GetInstance()->GetParameters(); |
143 | isInitialised = true; |
144 | SetParameters(); |
145 | if(isActive) { |
146 | theFermiModel->Initialise(); |
147 | theEvaporation->InitialiseChannels(); |
148 | } |
149 | if (1==0) |
150 | param->Dump(); |
151 | } |
152 | |
153 | void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr, G4bool flag) |
154 | { |
155 | if(ptr && ptr != theEvaporation) { |
156 | delete theEvaporation; |
157 | theEvaporation = ptr; |
158 | thePhotonEvaporation = ptr->GetPhotonEvaporation(); |
159 | theEvaporation->SetFermiBreakUp(theFermiModel); |
160 | isEvapLocal = flag; |
161 | } |
162 | } |
163 | |
164 | void |
165 | G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr) |
166 | { |
167 | if(ptr && ptr != theMultiFragmentation) { |
168 | delete theMultiFragmentation; |
169 | theMultiFragmentation = ptr; |
170 | } |
171 | } |
172 | |
173 | void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr) |
174 | { |
175 | if(ptr && ptr != theFermiModel) { |
176 | delete theFermiModel; |
177 | theFermiModel = ptr; |
178 | theEvaporation->SetFermiBreakUp(theFermiModel); |
179 | } |
180 | } |
181 | |
182 | void |
183 | G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr) |
184 | { |
185 | if(ptr && ptr != thePhotonEvaporation) { |
186 | thePhotonEvaporation = ptr; |
187 | theEvaporation->SetPhotonEvaporation(ptr); |
188 | } |
189 | } |
190 | |
191 | void G4ExcitationHandler::SetDeexChannelsType(G4DeexChannelType val) |
192 | { |
193 | G4Evaporation* evap = static_cast<G4Evaporation*>(theEvaporation); |
194 | if(val == fDummy) { |
195 | isActive = false; |
196 | return; |
197 | } |
198 | if(!evap) { return; } |
199 | if(val == fEvaporation) { |
200 | evap->SetDefaultChannel(); |
201 | } else if(val == fCombined) { |
202 | evap->SetCombinedChannel(); |
203 | } else if(val == fGEM) { |
204 | evap->SetGEMChannel(); |
205 | } |
206 | evap->InitialiseChannels(); |
207 | if(G4Threading::IsMasterThread()) { |
208 | G4cout(*G4cout_p) << "Number of de-excitation channels is changed to " |
209 | << theEvaporation->GetNumberOfChannels(); |
210 | if(fVerbose > 0) { G4cout(*G4cout_p) << " " << this; } |
211 | G4cout(*G4cout_p) << G4endlstd::endl; |
212 | } |
213 | } |
214 | |
215 | G4ReactionProductVector * |
216 | G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) |
217 | { |
218 | |
219 | G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState); |
220 | if(fVerbose > 1) { |
221 | G4cout(*G4cout_p) << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " << G4endlstd::endl; |
222 | G4cout(*G4cout_p) << theInitialState << G4endlstd::endl; |
223 | } |
224 | if(!isInitialised) { Initialise(); } |
225 | |
226 | |
227 | G4FragmentVector * theTempResult = nullptr; |
228 | |
229 | theResults.clear(); |
230 | thePhotoEvapList.clear(); |
231 | theEvapList.clear(); |
232 | |
233 | |
234 | G4double exEnergy = theInitialState.GetExcitationEnergy(); |
235 | G4int A = theInitialState.GetA_asInt(); |
236 | G4int Z = theInitialState.GetZ_asInt(); |
237 | |
238 | |
239 | if (A <= 1 || !isActive) { |
240 | theResults.push_back( theInitialStatePtr ); |
241 | |
242 | |
243 | } else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) { |
244 | theResults.push_back( theInitialStatePtr ); |
245 | |
246 | |
247 | |
248 | } else { |
249 | if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp) |
250 | || exEnergy <= minEForMultiFrag*A) { |
251 | theEvapList.push_back(theInitialStatePtr); |
252 | |
253 | |
254 | } else { |
255 | theTempResult = theMultiFragmentation->BreakItUp(theInitialState); |
256 | if(!theTempResult) { |
257 | theEvapList.push_back(theInitialStatePtr); |
258 | } else { |
259 | size_t nsec = theTempResult->size(); |
260 | |
261 | |
262 | if(0 == nsec) { |
263 | theEvapList.push_back(theInitialStatePtr); |
264 | |
265 | |
266 | } else { |
267 | G4bool deletePrimary = true; |
268 | G4FragmentVector::iterator j; |
269 | for (j = theTempResult->begin(); j != theTempResult->end(); ++j) { |
270 | |
271 | if((*j) == theInitialStatePtr) { deletePrimary = false; } |
272 | A = (*j)->GetA_asInt(); |
273 | |
274 | |
275 | if(A <= 1) { |
276 | theResults.push_back(*j); |
277 | |
278 | |
279 | } else { |
280 | G4double exEnergy1 = (*j)->GetExcitationEnergy(); |
281 | |
282 | |
283 | if(exEnergy1 < minExcitation) { |
284 | Z = (*j)->GetZ_asInt(); |
285 | if(nist->GetIsotopeAbundance(Z, A) > 0.0) { |
286 | theResults.push_back(*j); |
287 | } else { |
288 | theEvapList.push_back(*j); |
289 | } |
290 | |
291 | } else { |
292 | theEvapList.push_back(*j); |
293 | } |
294 | } |
295 | } |
296 | if( deletePrimary ) { delete theInitialStatePtr; } |
297 | } |
298 | delete theTempResult; |
299 | } |
300 | } |
301 | } |
302 | if(fVerbose > 2) { |
303 | G4cout(*G4cout_p) << "## After first step " << theEvapList.size() << " for evap; " |
304 | << thePhotoEvapList.size() << " for photo-evap; " |
305 | << theResults.size() << " results. " << G4endlstd::endl; |
306 | } |
307 | |
308 | |
309 | |
310 | |
311 | static const G4int countmax = 1000; |
312 | G4Fragment* frag; |
313 | size_t kk; |
314 | for (kk=0; kk<theEvapList.size(); ++kk) { |
315 | frag = theEvapList[kk]; |
316 | if(fVerbose > 2) { |
317 | G4cout(*G4cout_p) << "Next evaporate: " << G4endlstd::endl; |
318 | G4cout(*G4cout_p) << *frag << G4endlstd::endl; |
319 | } |
320 | if(kk >= countmax) { |
321 | G4ExceptionDescription ed; |
322 | ed << "Infinite loop in the de-excitation module: " << kk |
323 | << " iterations \n" |
324 | << " Initial fragment: \n" << theInitialState |
325 | << "\n Current fragment: \n" << *frag; |
326 | G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException, |
327 | ed,"Stop execution"); |
328 | |
329 | } |
330 | A = frag->GetA_asInt(); |
331 | Z = frag->GetZ_asInt(); |
332 | results.clear(); |
333 | |
334 | |
335 | if(theFermiModel->IsApplicable(Z, A, frag->GetExcitationEnergy())) { |
336 | theFermiModel->BreakFragment(&results, frag); |
337 | size_t nsec = results.size(); |
338 | if(fVerbose > 2) { G4cout(*G4cout_p) << "FermiBreakUp Nsec= " << nsec << G4endlstd::endl; } |
339 | |
340 | |
341 | |
342 | for(size_t j=0; j<nsec; ++j) { |
343 | exEnergy = results[j]->GetExcitationEnergy(); |
344 | if(exEnergy < minExcitation) { theResults.push_back(results[j]); } |
345 | else { thePhotoEvapList.push_back(results[j]); } |
346 | } |
347 | continue; |
348 | } |
349 | |
350 | theEvaporation->BreakFragment(&results, frag); |
351 | size_t nsec = results.size(); |
352 | if(fVerbose > 2) { G4cout(*G4cout_p) << "Evaporation Nsec= " << nsec << G4endlstd::endl; } |
353 | |
354 | |
355 | if(1 >= nsec) { |
356 | theResults.push_back(frag); |
357 | continue; |
358 | } |
359 | |
360 | |
361 | for (size_t j = 0; j<nsec; ++j) { |
362 | if(fVerbose > 3) { |
363 | G4cout(*G4cout_p) << "Evaporated product #" << j << G4endlstd::endl; |
364 | G4cout(*G4cout_p) << results[j] << G4endlstd::endl; |
365 | } |
366 | A = results[j]->GetA_asInt(); |
367 | if(A <= 1) { |
368 | theResults.push_back(results[j]); |
369 | continue; |
370 | } |
371 | exEnergy = results[j]->GetExcitationEnergy(); |
372 | |
373 | |
374 | if(exEnergy >= minExcitation) { |
375 | theEvapList.push_back(results[j]); |
376 | |
377 | |
378 | } else { |
379 | Z = results[j]->GetZ_asInt(); |
380 | |
381 | |
382 | if(nist->GetIsotopeAbundance(Z, A) > 0.0) { |
383 | theResults.push_back(results[j]); |
384 | |
385 | } else { |
386 | theEvapList.push_back(results[j]); |
387 | } |
388 | } |
389 | } |
390 | } |
391 | if(fVerbose > 2) { |
392 | G4cout(*G4cout_p) << "## After 2nd step " << theEvapList.size() << " was evap; " |
393 | << thePhotoEvapList.size() << " for photo-evap; " |
394 | << theResults.size() << " results. " << G4endlstd::endl; |
395 | } |
396 | |
397 | |
398 | |
399 | |
400 | |
401 | size_t kkmax = thePhotoEvapList.size(); |
402 | for (kk=0; kk<kkmax; ++kk) { |
403 | frag = thePhotoEvapList[kk]; |
404 | if(fVerbose > 2) { |
405 | G4cout(*G4cout_p) << "Next photon evaporate: " << thePhotonEvaporation << G4endlstd::endl; |
406 | G4cout(*G4cout_p) << *frag << G4endlstd::endl; |
407 | } |
408 | exEnergy = frag->GetExcitationEnergy(); |
409 | |
410 | |
411 | if(exEnergy > minExcitation) { |
412 | thePhotonEvaporation->BreakUpChain(&theResults, frag); |
413 | } |
414 | |
415 | |
416 | theResults.push_back(frag); |
417 | |
418 | } |
419 | if(fVerbose > 2) { |
420 | G4cout(*G4cout_p) << "## After 3d step " << theEvapList.size() << " was evap; " |
421 | << thePhotoEvapList.size() << " was photo-evap; " |
422 | << theResults.size() << " results. " << G4endlstd::endl; |
423 | } |
424 | G4ReactionProductVector * theReactionProductVector = |
425 | new G4ReactionProductVector(); |
426 | |
427 | |
428 | |
429 | theReactionProductVector->reserve( theResults.size() ); |
430 | |
431 | G4int theFragmentA, theFragmentZ; |
432 | |
433 | if(fVerbose > 1) { |
434 | G4cout(*G4cout_p) << "### ExcitationHandler provides " << theResults.size() |
435 | << " evaporated products:" << G4endlstd::endl; |
436 | } |
437 | kkmax = theResults.size(); |
438 | for (kk=0; kk<kkmax; ++kk) { |
439 | frag = theResults[kk]; |
440 | |
441 | |
442 | |
443 | if(!isActive && 0 == kk) { |
444 | G4double mass = frag->GetGroundStateMass(); |
445 | G4double ptot = (frag->GetMomentum()).vect().mag(); |
446 | G4double etot = (frag->GetMomentum()).e(); |
447 | G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0 |
448 | : std::sqrt((etot - mass)*(etot + mass))/ptot; |
449 | G4LorentzVector lv((frag->GetMomentum()).px()*fac, |
450 | (frag->GetMomentum()).py()*fac, |
451 | (frag->GetMomentum()).pz()*fac, etot); |
452 | frag->SetMomentum(lv); |
453 | } |
454 | if(fVerbose > 1) { |
455 | G4cout(*G4cout_p) << kk << "-th fragment " << frag; |
456 | if(frag->NuclearPolarization()) { |
457 | G4cout(*G4cout_p) << " " << frag->NuclearPolarization(); |
458 | } |
459 | G4cout(*G4cout_p) << G4endlstd::endl; |
460 | G4cout(*G4cout_p) << *frag << G4endlstd::endl; |
461 | } |
462 | |
463 | theFragmentA = frag->GetA_asInt(); |
464 | theFragmentZ = frag->GetZ_asInt(); |
465 | G4double etot= frag->GetMomentum().e(); |
466 | G4double eexc = 0.0; |
467 | const G4ParticleDefinition* theKindOfFragment = nullptr; |
468 | if (theFragmentA == 0) { |
469 | theKindOfFragment = frag->GetParticleDefinition(); |
470 | } else if (theFragmentA == 1 && theFragmentZ == 0) { |
471 | theKindOfFragment = G4Neutron::NeutronDefinition(); |
472 | } else if (theFragmentA == 1 && theFragmentZ == 1) { |
473 | theKindOfFragment = G4Proton::ProtonDefinition(); |
474 | } else if (theFragmentA == 2 && theFragmentZ == 1) { |
475 | theKindOfFragment = G4Deuteron::DeuteronDefinition(); |
476 | } else if (theFragmentA == 3 && theFragmentZ == 1) { |
477 | theKindOfFragment = G4Triton::TritonDefinition(); |
478 | } else if (theFragmentA == 3 && theFragmentZ == 2) { |
479 | theKindOfFragment = G4He3::He3Definition(); |
480 | } else if (theFragmentA == 4 && theFragmentZ == 2) { |
481 | theKindOfFragment = G4Alpha::AlphaDefinition();; |
482 | } else { |
483 | |
484 | |
485 | eexc = frag->GetExcitationEnergy(); |
486 | G4int idxf = frag->GetFloatingLevelNumber(); |
487 | if(eexc < minExcitation) { |
488 | eexc = 0.0; |
489 | idxf = 0; |
490 | } |
491 | |
492 | theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc, |
493 | G4Ions::FloatLevelBase(idxf)); |
494 | if(fVerbose > 1) { |
495 | G4cout(*G4cout_p) << "### EXCH: Find ion Z= " << theFragmentZ << " A= " << theFragmentA |
496 | << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf |
497 | << " " << theKindOfFragment << G4endlstd::endl; |
498 | } |
499 | } |
500 | |
501 | if(theKindOfFragment) { |
502 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
503 | theNew->SetMomentum(frag->GetMomentum().vect()); |
504 | theNew->SetTotalEnergy(etot); |
505 | theNew->SetFormationTime(frag->GetCreationTime()); |
506 | if(theKindOfFragment == electron) { theNew->SetCreatorModel(icID); } |
507 | theReactionProductVector->push_back(theNew); |
508 | |
509 | |
510 | } else { |
511 | theKindOfFragment = |
512 | theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0,noFloatG4Ions::G4FloatLevelBase::no_Float,0); |
513 | if(theKindOfFragment) { |
514 | G4ThreeVector mom(0.0,0.0,0.0); |
515 | G4double ionmass = theKindOfFragment->GetPDGMass(); |
516 | if(etot <= ionmass) { |
517 | etot = ionmass; |
518 | } else { |
519 | G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass)); |
520 | mom = (frag->GetMomentum().vect().unit())*ptot; |
521 | } |
522 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
523 | theNew->SetMomentum(mom); |
524 | theNew->SetTotalEnergy(etot); |
525 | theNew->SetFormationTime(frag->GetCreationTime()); |
526 | theReactionProductVector->push_back(theNew); |
527 | if(fVerbose > 2) { |
528 | G4cout(*G4cout_p) << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA |
529 | << " ground state, energy corrected E(MeV)= " << etot << G4endlstd::endl; |
530 | } |
531 | } |
532 | } |
533 | delete frag; |
534 | if(fVerbose > 1) { G4cout(*G4cout_p) << "G4Fragment #" << kk << " is deleted" << G4endlstd::endl; } |
535 | } |
536 | if(fVerbose > 2) { |
537 | G4cout(*G4cout_p) << "@@@@@@@@@@ End G4Excitation Handler "<< G4endlstd::endl; |
538 | } |
539 | return theReactionProductVector; |
540 | } |
541 | |
542 | void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const |
543 | { |
544 | outFile << "G4ExcitationHandler description\n" |
545 | << "This class samples de-excitation of excited nucleus using\n" |
546 | << "Fermi Break-up model for light fragments (Z < 9, A < 17), " |
547 | << "evaporation, fission, and photo-evaporation models. Evaporated\n" |
548 | << "particle may be proton, neutron, and other light fragment \n" |
549 | << "(Z < 13, A < 29). During photon evaporation produced gamma \n" |
550 | << "or electrons due to internal conversion \n"; |
551 | } |
552 | |
553 | |
554 | |
555 | |
556 | |
557 | |