- Timestamp:
- Apr 6, 2009, 12:30:29 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc
r819 r962 27 27 // Hadronic Process: Nuclear De-excitations 28 28 // by V. Lara (May 1998) 29 // 30 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 31 // cross section option 32 // 33 // Modif (24 Jul 2008) by M. A. Cortes Giraldo: 34 // -Max Z,A for Fermi Break-Up turns to 9,17 by default 35 // -BreakItUp() reorganised and bug in Evaporation loop fixed 36 // -Transform() optimised 29 37 // Modif (30 June 1998) by V. Lara: 30 38 // -Modified the Transform method for use G4ParticleTable and … … 33 41 // G4DynamicParticle 34 42 // -It uses default algorithms for: 35 // Evaporation: G4StatEvaporation 36 // MultiFragmentation: G4DummyMF (a dummy one) 37 // Fermi Breakup model: G4StatFermiBreakUp 38 43 // Evaporation: G4Evaporation 44 // MultiFragmentation: G4StatMF 45 // Fermi Breakup model: G4FermiBreakUp 46 // 47 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 48 // cross section option (default OPTxs=3) 49 // JMQ (06 September 2008) Also external choices have been added for 50 // superimposed Coulomb barrier (if useSICBis set true, by default is false) 39 51 40 52 #include "G4ExcitationHandler.hh" … … 45 57 46 58 G4ExcitationHandler::G4ExcitationHandler(): 59 // Fermi BreakUp is on and MultiFrag is off by default 60 // maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV), 47 61 maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV), 48 62 MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true), 49 MyOwnPhotonEvaporationClass(true) 63 MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false) 50 64 { 51 65 theTableOfParticles = G4ParticleTable::GetParticleTable(); … … 92 106 } 93 107 94 95 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment &theInitialState) const 96 { 97 98 G4FragmentVector * theResult = 0; 108 //////////////////////////////////////////////////////////////////////////////////////////////// 109 /// 25/07/08 16:45 Proposed by MAC //////////////////////////////////////////////////////////// 110 //////////////////////////////////////////////////////////////////////////////////////////////// 111 112 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const 113 { 114 //for inverse cross section choice 115 theEvaporation->SetOPTxs(OPTxs); 116 //for the choice of superimposed Coulomb Barrier for inverse cross sections 117 theEvaporation->UseSICB(useSICB); 118 119 // Pointer which will be used to return the final production vector 120 G4FragmentVector * theResult = 0; 121 122 // Variables existing until end of method 123 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState); 124 G4Fragment theExcitedNucleus; // object to be passed in BreakItUp methods 125 G4FragmentVector * theTempResult = 0; // pointer which receives temporal results 126 std::list<G4Fragment*> theEvapList; // list to apply Evaporation, SMF or Fermi Break-Up 127 std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation 128 std::list<G4Fragment*> theFinalStableList; // list to store final result 129 std::list<G4Fragment*>::iterator iList; 130 131 // Variables to describe the excited configuration 99 132 G4double exEnergy = theInitialState.GetExcitationEnergy(); 100 // G4cout << " first exEnergy in MeV: " << exEnergy/MeV << G4endl; 101 G4double A = theInitialState.GetA(); 102 G4int Z = static_cast<G4int>(theInitialState.GetZ()); 103 G4FragmentVector* theTempResult = 0; 104 G4Fragment theExcitedNucleus; 105 106 // Test applicability 107 if (A > 4) 133 G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 ); 134 G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 ); 135 136 // In case A <= 4 the fragment will not perform any nucleon emission 137 if (A <= 4) 108 138 { 109 // Initial State De-Excitation 110 if(A<GetMaxA()&&Z<GetMaxZ()) 111 // && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A)) { 112 { 113 theResult = theFermiModel->BreakItUp(theInitialState); 114 } 115 else if (exEnergy>GetMinE()*A) 116 { 117 theResult = theMultiFragmentation->BreakItUp(theInitialState); 118 } 119 else 120 { 121 theResult = theEvaporation->BreakItUp(theInitialState); 122 } 123 124 125 126 127 // De-Excitation loop 128 // ------------------ 129 // Check if there are excited fragments 130 std::list<G4Fragment*> theResultList; 131 G4FragmentVector::iterator j; 132 std::list<G4Fragment*>::iterator i; 133 for (j = theResult->begin(); j != theResult->end();j++) 134 { 135 theResultList.push_back(*j); 136 } 137 theResult->clear(); 138 for (i = theResultList.begin(); i != theResultList.end(); i++) 139 { 140 exEnergy = (*i)->GetExcitationEnergy(); 141 // G4cout << " exEnergy in MeV: " << exEnergy/MeV << G4endl; 142 if (exEnergy > 0.0) 143 { 144 A = (*i)->GetA(); 145 Z = static_cast<G4int>((*i)->GetZ()); 146 theExcitedNucleus = *(*i); 147 // try to de-excite this fragment 148 if( A < GetMaxA() && Z < GetMaxZ() ) 149 // && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A)) 150 { 151 // Fermi Breakup not now called for for exotic fragments for good reasons... 152 // theTempResult = theFermiModel->BreakItUp(theExcitedNucleus); 153 //if (theTempResult->size() == 1) 154 // { 155 // std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete()); 156 // delete theTempResult; 157 // } 158 theTempResult = theEvaporation->BreakItUp(theExcitedNucleus); 159 } 160 else 161 { 162 // Evaporation 163 theTempResult = theEvaporation->BreakItUp(theExcitedNucleus); 164 } 165 // The Nucleus has been fragmented? 166 if (theTempResult->size() > 1) 167 // If so : 168 { 169 // Remove excited fragment from the result 170 // delete theResult->removeAt(i--); 171 delete (*i); 172 i = theResultList.erase(i); 173 // and add theTempResult elements to theResult 174 for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin(); 175 ri != theTempResult->rend(); ++ri) 176 { 177 theResultList.push_back(*ri); 178 } 179 delete theTempResult; 180 } 181 else 182 // If not : 183 { 184 // it doesn't matter, we Follow with the next fragment but 185 // I have to clean up 186 std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete()); 187 delete theTempResult; 188 } 189 } 190 } 191 for (i = theResultList.begin(); i != theResultList.end(); i++) 192 { 193 theResult->push_back(*i); 194 } 195 theResultList.clear(); 139 // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later 140 141 theEvapStableList.push_back( theInitialStatePtr ); 196 142 } 197 else // if A > 4 143 144 else // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation 198 145 { 199 theResult = new G4FragmentVector(); 200 theResult->push_back(new G4Fragment(theInitialState)); 146 147 // Store original state in theEvapList 148 theEvapList.push_back( theInitialStatePtr ); 149 150 // ------------------------------ 151 // De-excitation loop 152 // ------------------------------ 153 154 for (iList = theEvapList.begin(); iList != theEvapList.end(); iList++) 155 { 156 exEnergy = (*iList)->GetExcitationEnergy(); 157 158 if (exEnergy > 0.0) 159 { 160 // Check conditions for each model 161 theExcitedNucleus = *(*iList); 162 A = static_cast<G4int>((*iList)->GetA()+0.5); // +0.5 to avoid bad truncation 163 Z = static_cast<G4int>((*iList)->GetZ()+0.5); 164 165 if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up 166 { 167 theTempResult = theFermiModel->BreakItUp(theExcitedNucleus); 168 } 169 else if (exEnergy > GetMinE()*A) // if satisfied apply SMF 170 { 171 theTempResult = theMultiFragmentation->BreakItUp(theExcitedNucleus); 172 } 173 else // apply Evaporation in another case 174 { 175 theTempResult = theEvaporation->BreakItUp(theExcitedNucleus); 176 } 177 178 // New configuration is stored in theTempResult, so we can free 179 // the memory where the previous configuration is 180 181 delete (*iList); 182 183 // And now the theTempResult->size() tells us if the configuration has changed 184 185 if ( theTempResult->size() > 1 ) 186 { 187 // push_back the result to the end of theEvapList (this same list) 188 189 for (G4FragmentVector::iterator j = theTempResult->begin(); 190 j != theTempResult->end(); j++) 191 { 192 theEvapList.push_back(*j); 193 } 194 } 195 else 196 { 197 // push_back the result to theEvapStableList, because 198 // is still excited, but cannot emmit more nucleons 199 200 for (G4FragmentVector::iterator j = theTempResult->begin(); 201 j != theTempResult->end(); j++) 202 { 203 theEvapStableList.push_back(*j); 204 } 205 } 206 207 // after working with theTempResult, clear and delete it 208 theTempResult->clear(); 209 delete theTempResult; 210 211 } 212 else // exEnergy = 0.0 213 { 214 // if this fragment is at ground state, 215 // store it in theFinalStableList 216 217 theFinalStableList.push_back(*iList); 218 219 } // endif (exEnergy > 0.0) 220 221 } // end of the loop over theEvapList 222 223 theEvapList.clear(); // clear all the list and do not free memory pointed by 224 // each element because this have been done before! 201 225 } 202 226 203 227 // Now we try to deexcite by means of PhotonEvaporation those fragments 204 // which are excited. 205 206 theTempResult = 0; 207 std::list<G4Fragment*> theFinalResultList; 208 //AHtest std::list<G4Fragment*> theFinalPhotonResultList; 209 std::list<G4Fragment*> theResultList; 210 std::list<G4Fragment*>::iterator j; 211 G4FragmentVector::iterator i; 212 for (i = theResult->begin(); i != theResult->end();i++) 228 // which are still excited. 229 230 231 // ----------------------- 232 // Photon-Evaporation loop 233 // ----------------------- 234 235 for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); iList++) 213 236 { 214 theResultList.push_back(*i); 215 // G4cout << " Before loop list energy in MeV: " << ((*i)->GetExcitationEnergy())/MeV << G4endl; 216 } 217 theResult->clear(); 218 219 for (j = theResultList.begin(); j != theResultList.end(); j++) { 220 // G4cout << " Test loop list: " << (*j)->GetExcitationEnergy() << " size: " << theResultList.size() << G4endl; 221 } 222 223 // for (j = theResultList.begin(); j != theResultList.end(); j++) 224 j = theResultList.begin(); //AH 225 while (j != theResultList.end()) //AH 226 { 227 if ((*j)->GetA() > 1 && (*j)->GetExcitationEnergy() > 0.1*eV) 228 { 229 theExcitedNucleus = *(*j); 230 theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus); 231 // If Gamma Evaporation has succeed then 232 if (theTempResult->size() > 1) 233 { 234 // Remove excited fragment from the result 235 // delete (*j); 236 // theResultList.erase(j--); 237 // theResultList.erase(j); don't delete as there's no push back... 238 // and add theTempResult elements to theResult 239 for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin(); 237 A = static_cast<G4int>((*iList)->GetA()+0.5); 238 exEnergy = (*iList)->GetExcitationEnergy(); 239 240 if ( A > 1 && exEnergy > 0.1*eV ) // if so, photon-evaporation is applied 241 { 242 theExcitedNucleus = *(*iList); 243 theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus); 244 245 // if there is a gamma emission then 246 if (theTempResult->size() > 1) 247 { 248 // first free the memory occupied by the previous state 249 delete (*iList); 250 251 // and now add the final state from gamma emission to the end of 252 // theEvapStableList 253 for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin(); 240 254 ri != theTempResult->rend(); ++ri) 241 { 255 // reversed is applied in order to have residual nucleus in first position 256 { 257 242 258 #ifdef PRECOMPOUND_TEST 243 if ((*ri)->GetA() == 0) 244 (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation")); 245 else 246 (*ri)->SetCreatorModel(G4String("ResidualNucleus")); 247 #endif 248 theResultList.push_back(*ri); 249 //AHtest theFinalPhotonResultList.push_back(*ri); 250 // theFinalResultList.push_back(*ri); don't add to final result as they'll go through the loop 251 } 252 delete theTempResult; 253 } 254 // In other case, just clean theTempResult and continue 255 else 256 { 257 std::for_each(theTempResult->begin(), theTempResult->end(), DeleteFragment()); 258 delete theTempResult; 259 if ((*ri)->GetA() == 0) 260 (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation")); 261 else 262 (*ri)->SetCreatorModel(G4String("ResidualNucleus")); 263 #endif 264 265 theEvapStableList.push_back(*ri); 266 } 267 268 // now we clean and remove the temporal vector 269 theTempResult->clear(); 270 delete theTempResult; 271 272 } 273 else // if theTempResult->size() = 1 274 { 275 // if there is not any gamma emission from this excited fragment 276 // we have to emmit a gamma which forces the deexcitation 277 278 // First I clean completely theTempResult 279 280 for (G4FragmentVector::iterator j = theTempResult->begin(); 281 j != theTempResult->end(); j++) 282 { 283 delete (*j); 284 } 285 286 theTempResult->clear(); 287 delete theTempResult; 288 259 289 #ifdef debugphoton 260 290 G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n" … … 263 293 << "-----------------------------------------------------------------------\n"; 264 294 #endif 265 G4double GammaEnergy = (*j)->GetExcitationEnergy(); 266 G4double cosTheta = 1. - 2. * G4UniformRand(); 267 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); 295 296 // Let's create a G4Fragment pointer representing the gamma emmited 297 G4double GammaEnergy = (*iList)->GetExcitationEnergy(); 298 G4double cosTheta = 1. - 2. * G4UniformRand(); 299 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); 268 300 G4double phi = twopi * G4UniformRand(); 269 301 G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi), … … 272 304 G4LorentzVector Gamma4P(GammaP,GammaEnergy); 273 305 G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition()); 274 275 276 277 G4double Mass = (*j)->GetGroundStateMass(); 278 G4ThreeVector ResidualP((*j)->GetMomentum().vect() - GammaP); 306 307 // And now we update momentum and energy for the nucleus 308 G4double Mass = (*iList)->GetGroundStateMass(); 309 G4ThreeVector ResidualP((*iList)->GetMomentum().vect() - GammaP); 279 310 G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass); 280 311 G4LorentzVector Residual4P(ResidualP,ResidualE); 281 (*j)->SetMomentum(Residual4P); 282 283 284 312 (*iList)->SetMomentum(Residual4P); // Now this fragment has been deexcited! 313 314 // we store the deexcited fragment in theFinalStableList 315 theFinalStableList.push_back(*iList); 316 285 317 #ifdef PRECOMPOUND_TEST 286 theHandlerPhoton->SetCreatorModel("G4ExcitationHandler"); 287 #endif 288 // theFinalPhotonResultList.push_back( theHandlerPhoton ); 289 // G4cout << " adding photon fragment " << G4endl; 290 theResultList.push_back( theHandlerPhoton ); 291 // theFinalResultList.push_back( theHandlerPhoton ); 292 theFinalResultList.push_back(*j); 318 theHandlerPhoton->SetCreatorModel("G4ExcitationHandler"); 319 #endif 320 321 // Finally, we add theHandlerPhoton to theFinalStableList 322 theFinalStableList.push_back(theHandlerPhoton); 323 293 324 #ifdef debugphoton 294 325 G4cout << "Emmited photon:\n" 295 << the ResultList.back() << '\n'326 << theFinalStableList.back() << '\n' 296 327 << "Residual nucleus after photon emission:\n" 297 << *(* j) << '\n'328 << *(*iList) << '\n' 298 329 << "-----------------------------------------------------------------------\n"; 299 330 #endif 300 //test j++; // AH only increment if not erased: 301 } 302 } else { 303 //test j++; // AH increment iterator if a proton or excitation energy small 304 theFinalResultList.push_back(*j); 331 332 } 305 333 } 306 // G4cout << " Inside loop list: " << (*j)->GetExcitationEnergy() << " size: " << theFinalResultList.size() << G4endl; 307 j++; 334 else // case of a nucleon, gamma or very small excitation energy 335 { 336 // we don't have to do anything, just store the fragment in theFinalStableList 337 338 theFinalStableList.push_back(*iList); 339 340 } // A > 1 && exEnergy > 0.1*eV 341 342 } // end of photon-evaporation loop 343 344 345 // The deexcitation from fragments inside theEvapStableList has been finished, so... 346 theEvapStableList.clear(); 347 348 // Now the final state is in theFinalStableList, and we have to send it to theResult vector 349 350 theResult = new G4FragmentVector; 351 theResult->reserve( theFinalStableList.size() * sizeof(G4Fragment*) ); 352 // We reserve enough memory to optimise the storing speed 353 354 for (iList = theFinalStableList.begin(); iList != theFinalStableList.end(); iList++) 355 { 356 theResult->push_back(*iList); 308 357 } 309 // for (j = theResultList.begin(); j != theResultList.end(); j++) 310 for (j = theFinalResultList.begin(); j != theFinalResultList.end(); j++) 311 { 312 theResult->push_back(*j); 313 } 314 315 //AHtest for (j = theFinalPhotonResultList.begin(); j != theFinalPhotonResultList.end(); j++) 316 //AHtest { 317 //AHtest theResult->push_back(*j); 318 //AHtest number_results++; 319 //AHtest } 320 321 322 theResultList.clear(); 323 theFinalResultList.clear(); 324 //AHtest theFinalPhotonResultList.clear(); 325 326 358 359 // After storing the final state , we can clear theFinalStableList 360 theFinalStableList.clear(); 361 327 362 #ifdef debug 328 363 CheckConservation(theInitialState,theResult); 329 364 #endif 330 // Change G4FragmentVector by G4DynamicParticle 365 366 // Change G4FragmentVector* to G4ReactionProductVector* 331 367 return Transform(theResult); 332 368 } 369 370 333 371 334 372 G4ReactionProductVector * … … 348 386 theNeutron->SetVerboseLevel(2); 349 387 G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector; 388 389 // MAC (24/07/08) 390 // To optimise the storing speed, we reserve space in memory for the vector 391 theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) ); 392 350 393 G4int theFragmentA, theFragmentZ; 351 394 G4LorentzVector theFragmentMomentum;
Note: See TracChangeset
for help on using the changeset viewer.