| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4ExcitationHandler.cc 79200 2014-02-20 11:22:39Z gcosmo $
>> 27 //
26 // Hadronic Process: Nuclear De-excitations 28 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (May 1998) 29 // by V. Lara (May 1998)
28 // 30 //
29 // 31 //
30 // Modified: 32 // Modified:
31 // 30 June 1998 by V. Lara: 33 // 30 June 1998 by V. Lara:
32 // -Modified the Transform method for use 34 // -Modified the Transform method for use G4ParticleTable and
33 // therefore G4IonTable. It makes possib 35 // therefore G4IonTable. It makes possible to convert all kind
34 // of fragments (G4Fragment) produced in 36 // of fragments (G4Fragment) produced in deexcitation to
35 // G4DynamicParticle 37 // G4DynamicParticle
36 // -It uses default algorithms for: 38 // -It uses default algorithms for:
37 // Evaporation: G4Evaporation 39 // Evaporation: G4Evaporation
38 // MultiFragmentation: G4StatMF 40 // MultiFragmentation: G4StatMF
39 // Fermi Breakup model: G4FermiBr 41 // Fermi Breakup model: G4FermiBreakUp
40 // 24 Jul 2008 by M. A. Cortes Giraldo: 42 // 24 Jul 2008 by M. A. Cortes Giraldo:
41 // -Max Z,A for Fermi Break-Up turns to 9 43 // -Max Z,A for Fermi Break-Up turns to 9,17 by default
42 // -BreakItUp() reorganised and bug in Ev 44 // -BreakItUp() reorganised and bug in Evaporation loop fixed
43 // -Transform() optimised 45 // -Transform() optimised
44 // (September 2008) by J. M. Quesada. External 46 // (September 2008) by J. M. Quesada. External choices have been added for :
45 // -inverse cross section option (default 47 // -inverse cross section option (default OPTxs=3)
46 // -superimposed Coulomb barrier (if useS 48 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
47 // September 2009 by J. M. Quesada: 49 // September 2009 by J. M. Quesada:
48 // -according to Igor Pshenichnov, SMM wi 50 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
49 // 27 Nov 2009 by V.Ivanchenko: 51 // 27 Nov 2009 by V.Ivanchenko:
50 // -cleanup the logic, reduce number inte 52 // -cleanup the logic, reduce number internal vectors, fixed memory leak.
51 // 11 May 2010 by V.Ivanchenko: 53 // 11 May 2010 by V.Ivanchenko:
52 // -FermiBreakUp activated, used integer 54 // -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
53 // final photon deexcitation; used check 55 // final photon deexcitation; used check on adundance of a fragment, decay
54 // unstable fragments with A <5 56 // unstable fragments with A <5
55 // 22 March 2011 by V.Ivanchenko: general clea 57 // 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition:
56 // products of Fermi Break Up cannot be 58 // products of Fermi Break Up cannot be further deexcited by this model
57 // 30 March 2011 by V.Ivanchenko removed priva 59 // 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods
58 // to the source 60 // to the source
59 // 23 January 2012 by V.Ivanchenko general cle 61 // 23 January 2012 by V.Ivanchenko general cleanup including destruction of
60 // objects, propagate G4PhotonEvaporation p 62 // objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and
61 // not delete it here 63 // not delete it here
62 64
>> 65 #include <list>
>> 66
63 #include "G4ExcitationHandler.hh" 67 #include "G4ExcitationHandler.hh"
64 #include "G4SystemOfUnits.hh" 68 #include "G4SystemOfUnits.hh"
65 #include "G4LorentzVector.hh" 69 #include "G4LorentzVector.hh"
66 #include "G4ThreeVector.hh" << 70 #include "G4NistManager.hh"
67 #include "G4ParticleTable.hh" 71 #include "G4ParticleTable.hh"
68 #include "G4ParticleTypes.hh" 72 #include "G4ParticleTypes.hh"
69 #include "G4Ions.hh" 73 #include "G4Ions.hh"
70 #include "G4Electron.hh" <<
71 #include "G4Lambda.hh" <<
72 74
73 #include "G4VMultiFragmentation.hh" 75 #include "G4VMultiFragmentation.hh"
74 #include "G4VFermiBreakUp.hh" 76 #include "G4VFermiBreakUp.hh"
75 #include "G4Element.hh" << 77 #include "G4VFermiFragment.hh"
76 #include "G4ElementTable.hh" <<
77 78
78 #include "G4VEvaporation.hh" 79 #include "G4VEvaporation.hh"
79 #include "G4VEvaporationChannel.hh" 80 #include "G4VEvaporationChannel.hh"
>> 81 #include "G4VPhotonEvaporation.hh"
80 #include "G4Evaporation.hh" 82 #include "G4Evaporation.hh"
81 #include "G4PhotonEvaporation.hh" <<
82 #include "G4StatMF.hh" 83 #include "G4StatMF.hh"
83 #include "G4FermiBreakUpVI.hh" << 84 #include "G4PhotonEvaporation.hh"
84 #include "G4NuclearLevelData.hh" << 85 #include "G4FermiBreakUp.hh"
85 #include "G4PhysicsModelCatalog.hh" << 86 #include "G4FermiFragmentsPool.hh"
86 << 87 #include "G4Pow.hh"
87 G4ExcitationHandler::G4ExcitationHandler() << 88
88 : minEForMultiFrag(1.*CLHEP::TeV), minExcita << 89 G4ExcitationHandler::G4ExcitationHandler():
89 maxExcitation(100.*CLHEP::MeV) << 90 maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4*GeV),
90 { << 91 minExcitation(keV),OPTxs(3),useSICB(false),isEvapLocal(true)
91 thePartTable = G4ParticleTable::GetParticleT << 92 {
92 theTableOfIons = thePartTable->GetIonTable() << 93 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
93 nist = G4NistManager::Instance(); <<
94 94
95 theMultiFragmentation = new G4StatMF(); << 95 theMultiFragmentation = new G4StatMF;
96 theFermiModel = new G4FermiBreakUpVI(); << 96 theFermiModel = new G4FermiBreakUp;
97 thePhotonEvaporation = new G4PhotonEvaporati << 97 thePhotonEvaporation = new G4PhotonEvaporation("ExcitationHandler",fDelayedEmission);
98 SetEvaporation(new G4Evaporation(thePhotonEv << 98 theEvaporation = new G4Evaporation(thePhotonEvaporation);
99 theResults.reserve(60); << 99 thePool = G4FermiFragmentsPool::Instance();
100 results.reserve(30); << 100 SetParameters();
101 theEvapList.reserve(30); << 101 G4Pow::GetInstance();
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 handle <<
115 } 102 }
116 103
117 G4ExcitationHandler::~G4ExcitationHandler() 104 G4ExcitationHandler::~G4ExcitationHandler()
118 { 105 {
>> 106 if(isEvapLocal) { delete theEvaporation; }
119 delete theMultiFragmentation; 107 delete theMultiFragmentation;
120 delete theFermiModel; 108 delete theFermiModel;
121 if(isEvapLocal) { delete theEvaporation; } <<
122 } 109 }
123 110
124 void G4ExcitationHandler::SetParameters() 111 void G4ExcitationHandler::SetParameters()
125 { 112 {
126 G4NuclearLevelData* ndata = G4NuclearLevelDa << 113 //for inverse cross section choice
127 auto param = ndata->GetParameters(); << 114 theEvaporation->SetOPTxs(OPTxs);
128 isActive = true; << 115 //for the choice of superimposed Coulomb Barrier for inverse cross sections
129 // check if de-excitation is needed << 116 theEvaporation->UseSICB(useSICB);
130 if (fDummy == param->GetDeexChannelsType()) << 117 theEvaporation->Initialise();
131 isActive = false; << 118 }
132 } else { << 119
133 // upload data for elements used in geomet << 120 G4ReactionProductVector *
134 G4int Zmax = 20; << 121 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
135 const G4ElementTable* table = G4Element::G << 122 {
136 for (auto const & elm : *table) { Zmax = s << 123 //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl;
137 ndata->UploadNuclearLevelData(Zmax+1); << 124
138 } << 125 // Variables existing until end of method
139 minEForMultiFrag = param->GetMinExPerNucleou << 126 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
140 minExcitation = param->GetMinExcitation(); << 127
141 maxExcitation = param->GetPrecoHighEnergy(); << 128 // pointer to fragment vector which receives temporal results
142 << 129 G4FragmentVector * theTempResult = 0;
143 // allowing local debug printout << 130
144 fVerbose = std::max(fVerbose, param->GetVerb << 131 // list of fragments to apply Evaporation or Fermi Break-Up
145 if (isActive) { << 132 std::list<G4Fragment*> theEvapList;
146 if (nullptr == thePhotonEvaporation) { << 133
147 SetPhotonEvaporation(new G4PhotonEvapora << 134 // list of fragments to apply PhotonEvaporation
148 } << 135 std::list<G4Fragment*> thePhotoEvapList;
149 if (nullptr == theFermiModel) { << 136
150 SetFermiModel(new G4FermiBreakUpVI()); << 137 // list of fragments to store final result
151 } << 138 std::list<G4Fragment*> theResults;
152 if (nullptr == theMultiFragmentation) { << 139 //
153 SetMultiFragmentation(new G4StatMF()); << 140 //G4cout << theInitialState << G4endl;
>> 141
>> 142 // Variables to describe the excited configuration
>> 143 G4double exEnergy = theInitialState.GetExcitationEnergy();
>> 144 G4int A = theInitialState.GetA_asInt();
>> 145 G4int Z = theInitialState.GetZ_asInt();
>> 146
>> 147 G4NistManager* nist = G4NistManager::Instance();
>> 148
>> 149 // In case A <= 1 the fragment will not perform any nucleon emission
>> 150 if (A <= 1)
>> 151 {
>> 152 theResults.push_back( theInitialStatePtr );
>> 153 }
>> 154 // check if a fragment is stable
>> 155 else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
>> 156 {
>> 157 theResults.push_back( theInitialStatePtr );
>> 158 }
>> 159 else
>> 160 {
>> 161 // JMQ 150909: first step in de-excitation is treated separately
>> 162 // Fragments after the first step are stored in theEvapList
>> 163 // Statistical Multifragmentation will take place only once
>> 164 //
>> 165 // move to evaporation loop
>> 166 if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
>> 167 || exEnergy <= minEForMultiFrag*A)
>> 168 {
>> 169 theEvapList.push_back(theInitialStatePtr);
>> 170 }
>> 171 else
>> 172 {
>> 173 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
>> 174 if(!theTempResult) { theEvapList.push_back(theInitialStatePtr); }
>> 175 else {
>> 176 size_t nsec = theTempResult->size();
>> 177 if(0 == nsec) { theEvapList.push_back(theInitialStatePtr); }
>> 178 else {
>> 179 // secondary are produced
>> 180 // Sort out secondary fragments
>> 181 G4bool deletePrimary = true;
>> 182 G4FragmentVector::iterator j;
>> 183 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
>> 184 if((*j) == theInitialStatePtr) { deletePrimary = false; }
>> 185 A = (*j)->GetA_asInt();
>> 186
>> 187 // gamma, p, n
>> 188 if(A <= 1) { theResults.push_back(*j); }
>> 189
>> 190 // Analyse fragment A > 1
>> 191 else {
>> 192 G4double exEnergy1 = (*j)->GetExcitationEnergy();
>> 193
>> 194 // cold fragments
>> 195 if(exEnergy1 < minExcitation) {
>> 196 Z = (*j)->GetZ_asInt();
>> 197 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
>> 198 theResults.push_back(*j); // stable fragment
>> 199
>> 200 } else {
>> 201
>> 202 // check if the cold fragment is from FBU pool
>> 203 const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
>> 204 if(ffrag) {
>> 205 if(ffrag->IsStable()) { theResults.push_back(*j); }
>> 206 else { theEvapList.push_back(*j); }
>> 207
>> 208 // cold fragment may be unstable
>> 209 } else {
>> 210 theEvapList.push_back(*j);
>> 211 }
>> 212 }
>> 213
>> 214 // hot fragments are unstable
>> 215 } else { theEvapList.push_back(*j); }
>> 216 }
>> 217 }
>> 218 if( deletePrimary ) { delete theInitialStatePtr; }
>> 219 }
>> 220 delete theTempResult;
>> 221 }
>> 222 }
154 } 223 }
155 if (nullptr == theEvaporation) { << 224 /*
156 SetEvaporation(new G4Evaporation(thePhot << 225 G4cout << "## After first step " << theEvapList.size() << " for evap; "
>> 226 << thePhotoEvapList.size() << " for photo-evap; "
>> 227 << theResults.size() << " results. " << G4endl;
>> 228 */
>> 229 // -----------------------------------
>> 230 // FermiBreakUp and De-excitation loop
>> 231 // -----------------------------------
>> 232
>> 233 std::list<G4Fragment*>::iterator iList;
>> 234 for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
>> 235 {
>> 236 //G4cout << "Next evaporate: " << G4endl;
>> 237 //G4cout << *iList << G4endl;
>> 238 A = (*iList)->GetA_asInt();
>> 239 Z = (*iList)->GetZ_asInt();
>> 240
>> 241 // Fermi Break-Up
>> 242 G4bool wasFBU = false;
>> 243 if (A < maxAForFermiBreakUp && Z < maxZForFermiBreakUp)
>> 244 {
>> 245 theTempResult = theFermiModel->BreakItUp(*(*iList));
>> 246 wasFBU = true;
>> 247 // if initial fragment returned unchanged try to evaporate it
>> 248 if(1 == theTempResult->size()) {
>> 249 delete *(theTempResult->begin());
>> 250 delete theTempResult;
>> 251 theTempResult = theEvaporation->BreakItUp(*(*iList));
>> 252 }
>> 253 }
>> 254 else // apply Evaporation in another case
>> 255 {
>> 256 theTempResult = theEvaporation->BreakItUp(*(*iList));
>> 257 }
>> 258
>> 259 G4bool deletePrimary = true;
>> 260 size_t nsec = theTempResult->size();
>> 261 //G4cout << "Nproducts= " << nsec << G4endl;
>> 262
>> 263 // Sort out secondary fragments
>> 264 if ( nsec > 0 ) {
>> 265 G4FragmentVector::iterator j;
>> 266 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
>> 267 if((*j) == (*iList)) { deletePrimary = false; }
>> 268
>> 269 //G4cout << *j << G4endl;
>> 270 A = (*j)->GetA_asInt();
>> 271 exEnergy = (*j)->GetExcitationEnergy();
>> 272
>> 273 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n
>> 274
>> 275 // evaporation is not possible
>> 276 else if(1 == nsec) {
>> 277 if(exEnergy < minExcitation) { theResults.push_back(*j); }
>> 278 else { thePhotoEvapList.push_back(*j); }
>> 279
>> 280 } else { // Analyse fragment
>> 281
>> 282 // cold fragment
>> 283 if(exEnergy < minExcitation) {
>> 284 Z = (*j)->GetZ_asInt();
>> 285
>> 286 // natural isotope
>> 287 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
>> 288 theResults.push_back(*j); // stable fragment
>> 289
>> 290 } else {
>> 291 const G4VFermiFragment* ffrag = thePool->GetFragment(Z, A);
>> 292
>> 293 // isotope from FBU pool
>> 294 if(ffrag) {
>> 295 if(ffrag->IsStable()) { theResults.push_back(*j); }
>> 296 else { theEvapList.push_back(*j); }
>> 297
>> 298 // isotope may be unstable
>> 299 } else {
>> 300 theEvapList.push_back(*j);
>> 301 }
>> 302 }
>> 303
>> 304 // hot fragment
>> 305 } else if (wasFBU) {
>> 306 thePhotoEvapList.push_back(*j); // FBU applied only once
>> 307 } else {
>> 308 theEvapList.push_back(*j);
>> 309 }
>> 310 }
>> 311 }
>> 312 }
>> 313 if( deletePrimary ) { delete (*iList); }
>> 314 delete theTempResult;
>> 315 } // end of the loop over theEvapList
>> 316
>> 317 //G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
>> 318 // << thePhotoEvapList.size() << " for photo-evap; "
>> 319 // << theResults.size() << " results. " << G4endl;
>> 320
>> 321 // -----------------------
>> 322 // Photon-Evaporation loop
>> 323 // -----------------------
>> 324
>> 325 // at this point only photon evaporation is possible
>> 326 for(iList = thePhotoEvapList.begin(); iList != thePhotoEvapList.end(); ++iList)
>> 327 {
>> 328 //G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
>> 329 //G4cout << *iList << G4endl;
>> 330 exEnergy = (*iList)->GetExcitationEnergy();
>> 331
>> 332 // only hot fragments
>> 333 if(exEnergy > minExcitation) {
>> 334 theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);
>> 335 size_t nsec = theTempResult->size();
>> 336 //G4cout << "Nproducts= " << nsec << G4endl;
>> 337
>> 338 // if there is a gamma emission then
>> 339 if (nsec > 0)
>> 340 {
>> 341 G4FragmentVector::iterator j;
>> 342 for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
>> 343 {
>> 344 theResults.push_back(*j);
>> 345 }
>> 346 }
>> 347 delete theTempResult;
>> 348 }
>> 349
>> 350 // priamry fragment is kept
>> 351 theResults.push_back(*iList);
>> 352
>> 353 } // end of photon-evaporation loop
>> 354
>> 355 //G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
>> 356 // << thePhotoEvapList.size() << " was photo-evap; "
>> 357 // << theResults.size() << " results. " << G4endl;
>> 358
>> 359 G4ReactionProductVector * theReactionProductVector =
>> 360 new G4ReactionProductVector();
>> 361
>> 362 // MAC (24/07/08)
>> 363 // To optimise the storing speed, we reserve space in memory for the vector
>> 364 theReactionProductVector->reserve( theResults.size() );
>> 365
>> 366 G4int theFragmentA, theFragmentZ;
>> 367
>> 368 std::list<G4Fragment*>::iterator i;
>> 369 for (i = theResults.begin(); i != theResults.end(); ++i)
>> 370 {
>> 371 theFragmentA = (*i)->GetA_asInt();
>> 372 theFragmentZ = (*i)->GetZ_asInt();
>> 373 G4double etot= (*i)->GetMomentum().e();
>> 374 G4ParticleDefinition* theKindOfFragment = 0;
>> 375 if (theFragmentA == 0) { // photon or e-
>> 376 theKindOfFragment = (*i)->GetParticleDefinition();
>> 377 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
>> 378 theKindOfFragment = G4Neutron::NeutronDefinition();
>> 379 } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
>> 380 theKindOfFragment = G4Proton::ProtonDefinition();
>> 381 } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
>> 382 theKindOfFragment = G4Deuteron::DeuteronDefinition();
>> 383 } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
>> 384 theKindOfFragment = G4Triton::TritonDefinition();
>> 385 } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
>> 386 theKindOfFragment = G4He3::He3Definition();
>> 387 } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
>> 388 theKindOfFragment = G4Alpha::AlphaDefinition();;
>> 389 } else {
>> 390
>> 391 // ground state by default
>> 392 G4double eexc = (*i)->GetExcitationEnergy();
>> 393 G4double excitation = eexc;
>> 394
>> 395 G4int level = 0;
>> 396 theKindOfFragment =
>> 397 theTableOfIons->GetIon(theFragmentZ,theFragmentA,level);
>> 398 /*
>> 399 G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
>> 400 << " Eexc(MeV)= " << excitation/MeV << " "
>> 401 << theKindOfFragment << G4endl;
>> 402 */
>> 403 // production of an isomer
>> 404 if(eexc > minExcitation) {
>> 405 G4double elevel1 = 0.0;
>> 406 G4double elevel2 = 0.0;
>> 407 G4ParticleDefinition* ion = 0;
>> 408 for(level=1; level<9; ++level) {
>> 409 ion = theTableOfIons->GetIon(theFragmentZ,theFragmentA,level);
>> 410 //G4cout << level << " " << ion << G4endl;
>> 411 if(ion) {
>> 412 G4Ions* ip = dynamic_cast<G4Ions*>(ion);
>> 413 if(ip) {
>> 414 elevel2 = ip->GetExcitationEnergy();
>> 415 //G4cout<<" Level "<<level<<" E(MeV)= "<<elevel2/MeV<<G4endl;
>> 416 // close level
>> 417 if(std::fabs(eexc - elevel2) < minExcitation) {
>> 418 excitation = eexc - elevel2;
>> 419 theKindOfFragment = ion;
>> 420 break;
>> 421 // previous level was closer
>> 422 } else if(elevel2 - eexc >= eexc - elevel1) {
>> 423 excitation = eexc - elevel1;
>> 424 break;
>> 425 // will check next level and save current
>> 426 } else {
>> 427 theKindOfFragment = ion;
>> 428 excitation = eexc - elevel2;
>> 429 elevel1 = elevel2;
>> 430 }
>> 431 }
>> 432 } else {
>> 433 break;
>> 434 }
>> 435 }
>> 436 }
>> 437 // correction of total energy for ground state isotopes
>> 438 etot += excitation;
>> 439 G4double ionmass = theKindOfFragment->GetPDGMass();
>> 440 if(etot < ionmass) { etot = ionmass; }
>> 441 }
>> 442 if (theKindOfFragment != 0)
>> 443 {
>> 444 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
>> 445 theNew->SetMomentum((*i)->GetMomentum().vect());
>> 446 theNew->SetTotalEnergy(etot);
>> 447 theNew->SetFormationTime((*i)->GetCreationTime());
>> 448 theReactionProductVector->push_back(theNew);
>> 449 }
>> 450 delete (*i);
157 } 451 }
158 } <<
159 theFermiModel->SetVerbose(fVerbose); <<
160 if(fVerbose > 1) { <<
161 G4cout << "G4ExcitationHandler::SetParamet <<
162 } <<
163 } <<
164 452
165 void G4ExcitationHandler::Initialise() << 453 return theReactionProductVector;
166 { <<
167 if(isInitialised) { return; } <<
168 if(fVerbose > 1) { <<
169 G4cout << "G4ExcitationHandler::Initialise <<
170 } <<
171 G4DeexPrecoParameters* param = <<
172 G4NuclearLevelData::GetInstance()->GetPara <<
173 isInitialised = true; <<
174 SetParameters(); <<
175 if(isActive) { <<
176 theFermiModel->Initialise(); <<
177 theEvaporation->InitialiseChannels(); <<
178 } <<
179 // dump level is controlled by parameter cla <<
180 param->Dump(); <<
181 } 454 }
182 455
183 void G4ExcitationHandler::SetEvaporation(G4VEv << 456 void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr)
184 { 457 {
185 if(nullptr != ptr && ptr != theEvaporation) << 458 if(ptr && ptr != theEvaporation) {
>> 459 delete theEvaporation;
186 theEvaporation = ptr; 460 theEvaporation = ptr;
187 theEvaporation->SetPhotonEvaporation(thePh << 461 thePhotonEvaporation = ptr->GetPhotonEvaporation();
188 theEvaporation->SetFermiBreakUp(theFermiMo << 462 SetParameters();
189 isEvapLocal = flag; << 463 isEvapLocal = false;
190 if(fVerbose > 1) { <<
191 G4cout << "G4ExcitationHandler::SetEvapo <<
192 } <<
193 } 464 }
194 } 465 }
195 466
196 void 467 void
197 G4ExcitationHandler::SetMultiFragmentation(G4V 468 G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr)
198 { 469 {
199 if(nullptr != ptr && ptr != theMultiFragment << 470 if(ptr && ptr != theMultiFragmentation) {
200 delete theMultiFragmentation; 471 delete theMultiFragmentation;
201 theMultiFragmentation = ptr; 472 theMultiFragmentation = ptr;
202 } 473 }
203 } 474 }
204 475
205 void G4ExcitationHandler::SetFermiModel(G4VFer 476 void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr)
206 { 477 {
207 if(nullptr != ptr && ptr != theFermiModel) { << 478 if(ptr && ptr != theFermiModel) {
208 delete theFermiModel; 479 delete theFermiModel;
209 theFermiModel = ptr; 480 theFermiModel = ptr;
210 if(nullptr != theEvaporation) { <<
211 theEvaporation->SetFermiBreakUp(theFermi <<
212 } <<
213 } 481 }
214 } 482 }
215 483
216 void 484 void
217 G4ExcitationHandler::SetPhotonEvaporation(G4VE 485 G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
218 { 486 {
219 if(nullptr != ptr && ptr != thePhotonEvapora << 487 if(ptr && ptr != thePhotonEvaporation) {
220 delete thePhotonEvaporation; <<
221 thePhotonEvaporation = ptr; 488 thePhotonEvaporation = ptr;
222 if(nullptr != theEvaporation) { << 489 theEvaporation->SetPhotonEvaporation(ptr);
223 theEvaporation->SetPhotonEvaporation(ptr <<
224 } <<
225 if(fVerbose > 1) { <<
226 G4cout << "G4ExcitationHandler::SetPhoto <<
227 << " for handler " << this << G4e <<
228 } <<
229 } 490 }
230 } 491 }
231 492
232 void G4ExcitationHandler::SetDeexChannelsType( << 493 void G4ExcitationHandler::SetMaxZForFermiBreakUp(G4int aZ)
233 { 494 {
234 G4Evaporation* evap = static_cast<G4Evaporat << 495 maxZForFermiBreakUp = aZ;
235 if(fVerbose > 1) { <<
236 G4cout << "G4ExcitationHandler::SetDeexCha <<
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 chann <<
257 << theEvaporation->GetNumberOfChannels( <<
258 G4cout << " " << this; <<
259 } <<
260 G4cout << G4endl; <<
261 } <<
262 } <<
263 <<
264 G4VEvaporation* G4ExcitationHandler::GetEvapor <<
265 { <<
266 if (nullptr != theEvaporation) { SetParamete <<
267 return theEvaporation; <<
268 } 496 }
269 497
270 G4VMultiFragmentation* G4ExcitationHandler::Ge << 498 void G4ExcitationHandler::SetMaxAForFermiBreakUp(G4int anA)
271 { 499 {
272 if (nullptr != theMultiFragmentation) { SetP << 500 maxAForFermiBreakUp = anA;
273 return theMultiFragmentation; <<
274 } 501 }
275 502
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermi << 503 void G4ExcitationHandler::SetMaxAandZForFermiBreakUp(G4int anA, G4int aZ)
277 { 504 {
278 if (nullptr != theFermiModel) { SetParameter << 505 SetMaxAForFermiBreakUp(anA);
279 return theFermiModel; << 506 SetMaxZForFermiBreakUp(aZ);
280 } 507 }
281 508
282 G4VEvaporationChannel* G4ExcitationHandler::Ge << 509 void G4ExcitationHandler::SetMinEForMultiFrag(G4double anE)
283 { 510 {
284 if(nullptr != thePhotonEvaporation) { SetPar << 511 minEForMultiFrag = anE;
285 return thePhotonEvaporation; <<
286 } 512 }
287 513
288 G4ReactionProductVector * <<
289 G4ExcitationHandler::BreakItUp(const G4Fragmen <<
290 { <<
291 // Variables existing until end of method <<
292 G4Fragment * theInitialStatePtr = new G4Frag <<
293 if (fVerbose > 1) { <<
294 G4cout << "@@@@@@@@@@ Start G4Excitation H <<
295 G4cout << theInitialState << G4endl; <<
296 } <<
297 if (!isInitialised) { Initialise(); } <<
298 <<
299 // pointer to fragment vector which receives <<
300 G4FragmentVector * theTempResult = nullptr; <<
301 <<
302 theResults.clear(); <<
303 theEvapList.clear(); <<
304 <<
305 // Variables to describe the excited configu <<
306 G4double exEnergy = theInitialState.GetExcit <<
307 G4int A = theInitialState.GetA_asInt(); <<
308 G4int Z = theInitialState.GetZ_asInt(); <<
309 G4int nL = theInitialState.GetNumberOfLambda <<
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= " << <<
317 << " Eex/A(MeV)= " << exEnergy/A; <<
318 G4Exception("G4ExcitationHandler::BreakI <<
319 } <<
320 } <<
321 <<
322 // for hyper-nuclei subtract lambdas from th <<
323 G4double lambdaF = 0.0; <<
324 G4LorentzVector lambdaLV = theInitialStatePt <<
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 = thePartTa <<
348 if(nullptr != part) { <<
349 G4ReactionProduct* theNew = new G4Reaction <<
350 G4ThreeVector dir = G4ThreeVector( 0.0, 0. <<
351 if ( lambdaLV.vect().mag() > CLHEP:: <<
352 dir = lambdaLV.vect().unit(); <<
353 } <<
354 G4double mass = part->GetPDGMass(); <<
355 G4double etot = std::max(lambdaLV.e(), mas <<
356 dir *= std::sqrt((etot - mass)*(etot <<
357 theNew->SetMomentum(dir); <<
358 theNew->SetTotalEnergy(etot); <<
359 theNew->SetFormationTime(theInitialState.G <<
360 theNew->SetCreatorModelID(theInitialState. <<
361 G4ReactionProductVector* v = new G4Reactio <<
362 v->push_back(theNew); <<
363 return v; <<
364 } <<
365 } <<
366 } <<
367 G4double mass = theInitialStatePtr->GetGro <<
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron <<
369 <<
370 // de-excitation with neutrons instead of <<
371 theInitialStatePtr->SetZAandMomentum(lambd <<
372 514
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=" << <<
380 << " Eex/A(MeV)= " << exEnergy/A; <<
381 G4Exception("G4ExcitationHandler::BreakI <<
382 } <<
383 } <<
384 515
385 // In case A <= 1 the fragment will not perf <<
386 if (A <= 1 || !isActive) { <<
387 theResults.push_back( theInitialStatePtr ) <<
388 <<
389 // check if a fragment is stable <<
390 } else if (exEnergy < minExcitation && nist- <<
391 theResults.push_back( theInitialStatePtr ) <<
392 <<
393 // JMQ 150909: first step in de-excitation <<
394 // Fragments after the first step are stor <<
395 } else { <<
396 if ((A<maxAForFermiBreakUp && Z<maxZForFer <<
397 || exEnergy <= minEForMultiFrag*A) { <<
398 theEvapList.push_back(theInitialStatePtr <<
399 <<
400 // Statistical Multifragmentation will tak <<
401 } else { <<
402 theTempResult = theMultiFragmentation->B <<
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 secon <<
413 } else { <<
414 G4bool deletePrimary = true; <<
415 for (auto const & ptr : *theTempResult) { <<
416 if (ptr == theInitialStatePtr) { deleteP <<
417 SortSecondaryFragment(ptr); <<
418 } <<
419 if (deletePrimary) { delete theInitialStat <<
420 } <<
421 delete theTempResult; // end multifragmentat <<
422 } <<
423 } <<
424 } <<
425 if (fVerbose > 2) { <<
426 G4cout << "## After first step of handler <<
427 << " for evap; " <<
428 << theResults.size() << " results. " << G <<
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-excitatio <<
445 << " iterations \n" <<
446 << " Initial fragment: \n" << theIniti <<
447 << "\n Current fragment: \n" << *frag; <<
448 G4Exception("G4ExcitationHandler::BreakI <<
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 <<
457 << " Eex(MeV)= " << frag->GetExci <<
458 } <<
459 // Fermi Break-Up <<
460 if (theFermiModel->IsApplicable(Z, A, frag <<
461 theFermiModel->BreakFragment(&results, f <<
462 std::size_t nsec = results.size(); <<
463 if (fVerbose > 2) { G4cout << "FermiBrea <<
464 <<
465 // FBU takes care to delete input fragme <<
466 // The secondary may be excited - photo- <<
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 <<
476 // photon evaporation is possible <<
477 theEvaporation->BreakFragment(&results, fr <<
478 if (fVerbose > 3) { <<
479 G4cout << "Evaporation Nsec= " << result <<
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 << <<
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 " <<
497 << " was evap; " <<
498 << theResults.size() << " results. " << G <<
499 } <<
500 G4ReactionProductVector * theReactionProduct <<
501 new G4ReactionProductVector(); <<
502 <<
503 // MAC (24/07/08) <<
504 // To optimise the storing speed, we reserve <<
505 // in memory for the vector <<
506 theReactionProductVector->reserve( theResult <<
507 <<
508 if (fVerbose > 1) { <<
509 G4cout << "### ExcitationHandler provides <<
510 << " evaporated products:" << G4endl; <<
511 } <<
512 G4LorentzVector partOfLambdaLV; <<
513 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4d <<
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, exc <<
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 == <<
524 : std::sqrt((etot - mass)*(etot + mass))/pto <<
525 G4LorentzVector lv((frag->GetMomentum()) <<
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* theKindOfFragm <<
542 G4bool isHyperN = false; <<
543 if (fragmentA == 0) { // photon or e <<
544 theKindOfFragment = frag->GetParticleDef <<
545 } else if (fragmentA == 1 && fragmentZ == <<
546 theKindOfFragment = theNeutron; <<
547 } else if (fragmentA == 1 && fragmentZ == <<
548 theKindOfFragment = theProton; <<
549 } else if (fragmentA == 2 && fragmentZ == <<
550 theKindOfFragment = theDeuteron; <<
551 } else if (fragmentA == 3 && fragmentZ == <<
552 theKindOfFragment = theTriton; <<
553 if(0 < nL) { <<
554 const G4ParticleDefinition* p = thePar <<
555 if(nullptr != p) { <<
556 theKindOfFragment = p; <<
557 isHyperN = true; <<
558 --nL; <<
559 } <<
560 } <<
561 } else if (fragmentA == 3 && fragmentZ == <<
562 theKindOfFragment = theHe3; <<
563 } else if (fragmentA == 4 && fragmentZ == <<
564 theKindOfFragment = theAlpha; <<
565 if (0 < nL) { <<
566 const G4ParticleDefinition* p = thePar <<
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->GetFloatingLevelNumbe <<
578 if (eexc < minExcitation) { <<
579 eexc = 0.0; <<
580 idxf = 0; <<
581 } <<
582 <<
583 theKindOfFragment = theTableOfIons->GetI <<
584 <<
585 if (fVerbose > 3) { <<
586 G4cout << "### EXCH: Find ion Z= " << fragme <<
587 << " A= " << fragmentA <<
588 << " Eexc(MeV)= " << eexc/MeV << " id <<
589 << " " << theKindOfFragment->GetParti <<
590 << G4endl; <<
591 } <<
592 } <<
593 // fragment identified <<
594 if (nullptr != theKindOfFragment) { <<
595 G4ReactionProduct * theNew = new G4React <<
596 if (isHyperN) { <<
597 G4LorentzVector lv = lv0 + partOfLambd <<
598 G4ThreeVector dir = lv.vect().unit(); <<
599 G4double mass = theKindOfFragment->Get <<
600 etot = std::max(lv.e(), mass); <<
601 G4double ptot = std::sqrt((etot - mass <<
602 dir *= ptot; <<
603 theNew->SetMomentum(dir); <<
604 // remaining not compensated 4-momentum <<
605 lambdaLV += (lv0 - G4LorentzVector(dir <<
606 } else { <<
607 theNew->SetMomentum(lv0.vect()); <<
608 } <<
609 theNew->SetTotalEnergy(etot); <<
610 theNew->SetFormationTime(frag->GetCreati <<
611 theNew->SetCreatorModelID(frag->GetCreat <<
612 theReactionProductVector->push_back(theN <<
613 <<
614 // fragment not found out ground state i <<
615 } else { <<
616 theKindOfFragment = <<
617 theTableOfIons->GetIon(fragmentZ,fragmentA,0 <<
618 if (theKindOfFragment) { <<
619 G4ThreeVector mom(0.0,0.0,0.0); <<
620 G4double ionmass = theKindOfFragment->GetPDG <<
621 if (etot <= ionmass) { <<
622 etot = ionmass; <<
623 } else { <<
624 G4double ptot = std::sqrt((etot - ionmass) <<
625 mom = (frag->GetMomentum().vect().unit())* <<
626 } <<
627 G4ReactionProduct * theNew = new G4ReactionP <<
628 theNew->SetMomentum(mom); <<
629 theNew->SetTotalEnergy(etot); <<
630 theNew->SetFormationTime(frag->GetCreationTi <<
631 theNew->SetCreatorModelID(frag->GetCreatorMo <<
632 theReactionProductVector->push_back(theNew); <<
633 if (fVerbose > 3) { <<
634 G4cout << " ground state, energy <<
635 << etot << G4endl; <<
636 } <<
637 } <<
638 } <<
639 delete frag; <<
640 } <<
641 // remaining lambdas are free; conserve quan <<
642 // not 4-momentum <<
643 if (0 < nL) { <<
644 G4ThreeVector dir = G4ThreeVector(0.0, 0.0 <<
645 if (lambdaLV.vect().mag() > CLHEP::eV) { <<
646 dir = lambdaLV.vect().unit(); <<
647 } <<
648 G4double etot = std::max(lambdaLV.e()/(G4d <<
649 dir *= std::sqrt((etot - fLambdaMass)*(eto <<
650 for (G4int i=0; i<nL; ++i) { <<
651 G4ReactionProduct* theNew = new G4Reacti <<
652 theNew->SetMomentum(dir); <<
653 theNew->SetTotalEnergy(etot); <<
654 theNew->SetFormationTime(theInitialState <<
655 theNew->SetCreatorModelID(theInitialStat <<
656 theReactionProductVector->push_back(theN <<
657 } <<
658 } <<
659 if (fVerbose > 3) { <<
660 G4cout << "@@@@@@@@@@ End G4Excitation Han <<
661 } <<
662 return theReactionProductVector; <<
663 } <<
664 <<
665 void G4ExcitationHandler::ModelDescription(std <<
666 { <<
667 outFile << "G4ExcitationHandler description\ <<
668 << "This class samples de-excitation of ex <<
669 << "Fermi Break-up model for light fragmen <<
670 << "evaporation, fission, and photo-evapor <<
671 << "particle may be proton, neutron, and o <<
672 << "(Z < 13, A < 29). During photon evapor <<
673 << "or electrons due to internal conversio <<
674 } <<
675 516
676 517
677 518
678 519