Changeset 1315 for trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ExcitationHandler.cc,v 1.39 2010/05/18 15:46:11 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 26 28 // 27 29 // Hadronic Process: Nuclear De-excitations … … 29 31 // 30 32 // 31 // Modif (September 2009) by J. M. Quesada: 32 // according to Igor Pshenichnov, SMM will be applied (just in case) only once . 33 // 34 // Modif (September 2008) by J. M. Quesada. External choices have been added for : 35 // -inverse cross section option (default OPTxs=3) 36 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false) 37 // 38 // Modif (24 Jul 2008) by M. A. Cortes Giraldo: 39 // -Max Z,A for Fermi Break-Up turns to 9,17 by default 40 // -BreakItUp() reorganised and bug in Evaporation loop fixed 41 // -Transform() optimised 42 // Modif (30 June 1998) by V. Lara: 33 // Modified: 34 // (30 June 1998) by V. Lara: 43 35 // -Modified the Transform method for use G4ParticleTable and 44 36 // therefore G4IonTable. It makes possible to convert all kind … … 49 41 // MultiFragmentation: G4StatMF 50 42 // Fermi Breakup model: G4FermiBreakUp 43 // (24 Jul 2008) by M. A. Cortes Giraldo: 44 // -Max Z,A for Fermi Break-Up turns to 9,17 by default 45 // -BreakItUp() reorganised and bug in Evaporation loop fixed 46 // -Transform() optimised 47 // (September 2008) by J. M. Quesada. External choices have been added for : 48 // -inverse cross section option (default OPTxs=3) 49 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false) 50 // (September 2009) by J. M. Quesada: 51 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once. 52 // (27 Nov 2009) by V.Ivanchenko: 53 // -cleanup the logic, reduce number internal vectors, fixed memory leak. 54 // (11 May 2010) by V.Ivanchenko: 55 // -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for 56 // final photon deexcitation; used check on adundance of a fragment, decay 57 // unstable fragments with A <5 51 58 // 52 59 … … 54 61 #include "globals.hh" 55 62 #include "G4LorentzVector.hh" 63 #include "G4NistManager.hh" 56 64 #include <list> 57 58 //#define debugphoton59 60 65 61 66 G4ExcitationHandler::G4ExcitationHandler(): 62 67 // JMQ 160909 Fermi BreakUp & MultiFrag are on by default 63 68 // This is needed for activation of such models when G4BinaryLightIonReaction is used 64 // since no interface (for external activation via macro input file) is still available in this case. 65 //maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV), 66 maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV), 69 // since no interface (for external activation via macro input file) is still available. 70 maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV), 71 // maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV), 72 //maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV), 73 minExcitation(CLHEP::keV), 67 74 MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true), 68 75 MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false) … … 74 81 theFermiModel = new G4FermiBreakUp; 75 82 thePhotonEvaporation = new G4PhotonEvaporation; 83 SetParameters(); 76 84 } 77 78 G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)79 {80 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! ");81 }82 83 85 84 86 G4ExcitationHandler::~G4ExcitationHandler() … … 90 92 } 91 93 92 93 const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &) 94 void G4ExcitationHandler::SetParameters() 94 95 { 95 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");96 97 return *this;98 }99 100 101 G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const102 {103 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");104 return false;105 }106 107 G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const108 {109 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");110 return true;111 }112 113 ////////////////////////////////////////////////////////////////////////////////////////////////114 /// 25/07/08 16:45 Proposed by MAC ////////////////////////////////////////////////////////////115 ////////////////////////////////////////////////////////////////////////////////////////////////116 117 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const118 {119 96 //for inverse cross section choice 120 97 theEvaporation->SetOPTxs(OPTxs); 121 98 //for the choice of superimposed Coulomb Barrier for inverse cross sections 122 99 theEvaporation->UseSICB(useSICB); 123 124 // Pointer which will be used to return the final production vector 125 //G4FragmentVector * theResult = new G4FragmentVector; 100 theEvaporation->Initialise(); 101 } 102 103 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const 104 { 126 105 127 106 // Variables existing until end of method 128 //G4Fragment * theInitialStatePtr = const_cast<G4Fragment*>(&theInitialState);129 107 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState); 108 130 109 G4FragmentVector * theTempResult = 0; // pointer which receives temporal results 131 std::list<G4Fragment*> theEvapList; // list to apply Evaporation , SMFor Fermi Break-Up110 std::list<G4Fragment*> theEvapList; // list to apply Evaporation or Fermi Break-Up 132 111 std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation 133 112 std::list<G4Fragment*> theResults; // list to store final result 134 113 std::list<G4Fragment*>::iterator iList; 135 114 // 136 //G4cout << "@@@@@@@@@@ Start G4Ex itation Handler @@@@@@@@@@@@@" << G4endl;115 //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl; 137 116 //G4cout << theInitialState << G4endl; 138 117 139 118 // Variables to describe the excited configuration 140 119 G4double exEnergy = theInitialState.GetExcitationEnergy(); 141 G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 ); 142 G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 ); 120 G4int A = theInitialState.GetA_asInt(); 121 G4int Z = theInitialState.GetZ_asInt(); 122 123 G4NistManager* nist = G4NistManager::Instance(); 143 124 144 125 // JMQ 150909: first step in de-excitation chain (SMM will be used only here) 145 126 146 // In case A <= 4the fragment will not perform any nucleon emission147 if (A <= 4)127 // In case A <= 1 the fragment will not perform any nucleon emission 128 if (A <= 1) 148 129 { 149 // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later 150 theEvapStableList.push_back( theInitialStatePtr ); 130 theResults.push_back( theInitialStatePtr ); 151 131 } 152 153 else // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation132 // check if a fragment is stable 133 else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) 154 134 { 155 135 theResults.push_back( theInitialStatePtr ); 136 } 137 else // In all cases apply once theFermiModel, theMultiFragmentation or theEvaporation 138 { 156 139 // JMQ 150909: first step in de-excitation is treated separately 157 140 // Fragments after the first step are stored in theEvapList 158 // Statistical Multifragmentation will take place (just in case) only here141 // Statistical Multifragmentation will take place only once 159 142 // 160 // Test applicability 161 // Initial State De-Excitation 162 if(A<GetMaxA()&&Z<GetMaxZ()) 143 // Fermi Break-Up always first 144 if((A < 5) || (A<GetMaxA() && Z<GetMaxZ())) 163 145 { 164 146 theTempResult = theFermiModel->BreakItUp(theInitialState); 147 // fragment was not decaied try to evaporate 148 if(1 == theTempResult->size()) { 149 if((*theTempResult)[0] != theInitialStatePtr) { 150 delete (*theTempResult)[0]; 151 } 152 delete theTempResult; 153 theTempResult = theEvaporation->BreakItUp(theInitialState); 154 } 165 155 } 166 156 else if (exEnergy>GetMinE()*A) 167 157 { 168 158 theTempResult = theMultiFragmentation->BreakItUp(theInitialState); 159 // fragment was not decaied try to evaporate 160 if(1 == theTempResult->size()) { 161 if((*theTempResult)[0] != theInitialStatePtr) { 162 delete (*theTempResult)[0]; 163 } 164 delete theTempResult; 165 theTempResult = theEvaporation->BreakItUp(theInitialState); 166 } 169 167 } 170 168 else … … 174 172 175 173 G4bool deletePrimary = true; 176 if(theTempResult->size() > 0) 174 G4int nsec=theTempResult->size(); 175 if(nsec > 0) 177 176 { 178 // S tore original state in theEvapList177 // Sort out secondary fragments 179 178 G4FragmentVector::iterator j; 180 179 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) 181 180 { 182 181 if((*j) == theInitialStatePtr) { deletePrimary = false; } 183 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation 184 185 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n 186 else if(A <= 4) { theEvapStableList.push_back(*j); } // evaporation is not possible 187 else { theEvapList.push_back(*j); } // evaporation is possible 182 A = (*j)->GetA_asInt(); 183 184 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n 185 else if(1 == nsec) { theEvapStableList.push_back(*j); } // evaporation is not possible 186 else // Analyse fragment 187 { 188 G4double exEnergy = (*j)->GetExcitationEnergy(); 189 if(exEnergy < minExcitation) { 190 Z = (*j)->GetZ_asInt(); 191 if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 192 theResults.push_back(*j); // stable fragment 193 } else { 194 theEvapList.push_back(*j); // unstable cold fragment 195 } 196 } else { theEvapList.push_back(*j); } // hot fragment 197 } 188 198 } 189 199 } … … 204 214 for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList) 205 215 { 206 A = static_cast<G4int>((*iList)->GetA()+0.5); // +0.5 to avoid bad truncation207 Z = static_cast<G4int>((*iList)->GetZ()+0.5);216 A = (*iList)->GetA_asInt(); 217 Z = (*iList)->GetZ_asInt(); 208 218 209 // In case A <= 4 the fragment will not perform any nucleon emission210 if ( A <= 4)219 // Fermi Break-Up 220 if ((A < 5) || (A < GetMaxA() && Z < GetMaxZ())) 211 221 { 212 // storing G4Fragment* in theEvapStableList to apply thePhotonEvaporation later 213 theEvapStableList.push_back(*iList ); 222 theTempResult = theFermiModel->BreakItUp(*(*iList)); 223 // fragment was not decaied try to evaporate 224 if(1 == theTempResult->size()) { 225 if((*theTempResult)[0] != theInitialStatePtr) { delete (*theTempResult)[0]; } 226 delete theTempResult; 227 theTempResult = theEvaporation->BreakItUp(*(*iList)); 228 } 214 229 } 215 216 else // If A > 4 we try to apply theFermiModel or theEvaporation 217 { 218 // stable fragment 219 if ((*iList)->GetExcitationEnergy() <= 0.1*eV) 220 { 221 theResults.push_back(*iList); 230 else // apply Evaporation in another case 231 { 232 theTempResult = theEvaporation->BreakItUp(*(*iList)); 233 } 234 235 G4bool deletePrimary = true; 236 G4int nsec = theTempResult->size(); 237 238 // Sort out secondary fragments 239 if ( nsec > 0 ) 240 { 241 G4FragmentVector::iterator j; 242 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) 243 { 244 if((*j) == (*iList)) { deletePrimary = false; } 245 A = (*j)->GetA_asInt(); 246 247 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n 248 else if(1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation 249 else // Analyse fragment 250 { 251 G4double exEnergy = (*j)->GetExcitationEnergy(); 252 if(exEnergy < minExcitation) { 253 Z = (*j)->GetZ_asInt(); 254 if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 255 theResults.push_back(*j); // stable fragment 256 } else { 257 theEvapList.push_back(*j); // unstable cold fragment 258 } 259 } else { theEvapList.push_back(*j); } // hot fragment 260 } 222 261 } 223 else 224 { 225 if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up 226 { 227 theTempResult = theFermiModel->BreakItUp(*(*iList)); 228 } 229 else // apply Evaporation in another case 230 { 231 theTempResult = theEvaporation->BreakItUp(*(*iList)); 232 } 233 234 // New configuration is stored in theTempResult, so we can free 235 // the memory where the previous configuration is 236 237 G4bool deletePrimary = true; 238 G4int nsec = theTempResult->size(); 239 240 // The number of secondaries tells us if the configuration has changed 241 if ( nsec > 0 ) 242 { 243 G4FragmentVector::iterator j; 244 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) 245 { 246 if((*j) == (*iList)) { deletePrimary = false; } 247 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation 248 249 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n 250 else if(A <= 4 || 1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation 251 else { theEvapList.push_back(*j); } 252 } 253 } 254 if( deletePrimary ) { delete (*iList); } 255 delete theTempResult; 256 257 } 258 } // endif (A <=4) 262 } 263 if( deletePrimary ) { delete (*iList); } 264 delete theTempResult; 259 265 } // end of the loop over theEvapList 260 266 … … 267 273 // ----------------------- 268 274 275 // normally should not reach this point - it is kind of work around 269 276 for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList) 270 277 { 271 // take out stable particles and fragments272 A = static_cast<G4int>((*iList)->GetA()+0.5);273 if ( A <= 1 ) { theResults.push_back(*iList); }274 else if ((*iList)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*iList); } 275 276 else278 // photon-evaporation is applied 279 theTempResult = thePhotonEvaporation->BreakUpFragment(*iList); 280 G4int nsec = theTempResult->size(); 281 282 // if there is a gamma emission then 283 if (nsec > 0) 277 284 { 278 // photon-evaporation is applied 279 theTempResult = thePhotonEvaporation->BreakItUp(*(*iList)); 280 281 G4bool deletePrimary = true; 282 G4int nsec = theTempResult->size(); 283 284 // if there is a gamma emission then 285 if (nsec > 1) 285 G4FragmentVector::iterator j; 286 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) 286 287 { 287 G4FragmentVector::iterator j; 288 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) 289 { 290 if((*j) == (*iList)) { deletePrimary = false; } 291 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation 292 293 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n 294 else if((*j)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*j); } // stable fragment 295 else { theEvapStableList.push_back(*j); } 296 } 288 theResults.push_back(*j); 297 289 } 298 299 else if(1 == nsec) 300 { 301 G4FragmentVector::iterator j = theTempResult->begin(); 302 if((*j) == (*iList)) { deletePrimary = false; } 303 // Let's create a G4Fragment pointer representing the gamma emmited 304 G4LorentzVector lv = (*j)->GetMomentum(); 305 G4double Mass = (*j)->GetGroundStateMass(); 306 G4double Ecm = lv.m(); 307 if(Ecm - Mass > 0.1*eV) 308 { 309 G4ThreeVector bst = lv.boostVector(); 310 G4double GammaEnergy = 0.5*(Ecm - Mass)*(Ecm + Mass)/Ecm; 311 G4double cosTheta = 1. - 2. * G4UniformRand(); 312 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); 313 G4double phi = twopi * G4UniformRand(); 314 G4LorentzVector Gamma4P(GammaEnergy * sinTheta * std::cos(phi), 315 GammaEnergy * sinTheta * std::sin(phi), 316 GammaEnergy * cosTheta, 317 GammaEnergy); 318 Gamma4P.boost(bst); 319 G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition()); 320 theResults.push_back(theHandlerPhoton); 321 322 // And now we update momentum and energy for the nucleus 323 lv -= Gamma4P; 324 (*j)->SetMomentum(lv); // Now this fragment has been deexcited! 325 } 326 // we store the deexcited fragment 327 theResults.push_back(*j); 328 } 329 if( deletePrimary ) { delete (*iList); } 330 delete theTempResult; 331 } 290 } 291 delete theTempResult; 292 293 // priamry fragment is kept 294 theResults.push_back(*iList); 295 332 296 } // end of photon-evaporation loop 333 297 … … 356 320 theFragmentMomentum = (*i)->GetMomentum(); 357 321 G4ParticleDefinition* theKindOfFragment = 0; 358 if (theFragmentA == 0 && theFragmentZ == 0) { // photon359 theKindOfFragment = G4Gamma::GammaDefinition();322 if (theFragmentA == 0) { // photon or e- 323 theKindOfFragment = (*i)->GetParticleDefinition(); 360 324 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron 361 325 theKindOfFragment = G4Neutron::NeutronDefinition();
Note: See TracChangeset
for help on using the changeset viewer.
