[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // Hadronic Process: Nuclear De-excitations |
---|
| 28 | // by V. Lara (May 1998) |
---|
| 29 | // Modif (30 June 1998) by V. Lara: |
---|
| 30 | // -Modified the Transform method for use G4ParticleTable and |
---|
| 31 | // therefore G4IonTable. It makes possible to convert all kind |
---|
| 32 | // of fragments (G4Fragment) produced in deexcitation to |
---|
| 33 | // G4DynamicParticle |
---|
| 34 | // -It uses default algorithms for: |
---|
| 35 | // Evaporation: G4StatEvaporation |
---|
| 36 | // MultiFragmentation: G4DummyMF (a dummy one) |
---|
| 37 | // Fermi Breakup model: G4StatFermiBreakUp |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | #include "G4ExcitationHandler.hh" |
---|
| 41 | #include <list> |
---|
| 42 | |
---|
| 43 | //#define debugphoton |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | G4ExcitationHandler::G4ExcitationHandler(): |
---|
| 47 | maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV), |
---|
| 48 | MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true), |
---|
| 49 | MyOwnPhotonEvaporationClass(true) |
---|
| 50 | { |
---|
| 51 | theTableOfParticles = G4ParticleTable::GetParticleTable(); |
---|
| 52 | |
---|
| 53 | theEvaporation = new G4Evaporation; |
---|
| 54 | theMultiFragmentation = new G4StatMF; |
---|
| 55 | theFermiModel = new G4FermiBreakUp; |
---|
| 56 | thePhotonEvaporation = new G4PhotonEvaporation; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &) |
---|
| 60 | { |
---|
| 61 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! "); |
---|
| 62 | } |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | G4ExcitationHandler::~G4ExcitationHandler() |
---|
| 66 | { |
---|
| 67 | if (MyOwnEvaporationClass) delete theEvaporation; |
---|
| 68 | if (MyOwnMultiFragmentationClass) delete theMultiFragmentation; |
---|
| 69 | if (MyOwnFermiBreakUpClass) delete theFermiModel; |
---|
| 70 | if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation; |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &) |
---|
| 75 | { |
---|
| 76 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! "); |
---|
| 77 | |
---|
| 78 | return *this; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const |
---|
| 83 | { |
---|
| 84 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! "); |
---|
| 85 | return false; |
---|
| 86 | } |
---|
| 87 | |
---|
| 88 | G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const |
---|
| 89 | { |
---|
| 90 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! "); |
---|
| 91 | return true; |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment &theInitialState) const |
---|
| 96 | { |
---|
| 97 | |
---|
| 98 | G4FragmentVector * theResult = 0; |
---|
| 99 | 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) |
---|
| 108 | { |
---|
| 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(); |
---|
| 196 | } |
---|
| 197 | else // if A > 4 |
---|
| 198 | { |
---|
| 199 | theResult = new G4FragmentVector(); |
---|
| 200 | theResult->push_back(new G4Fragment(theInitialState)); |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | // 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++) |
---|
| 213 | { |
---|
| 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(); |
---|
| 240 | ri != theTempResult->rend(); ++ri) |
---|
| 241 | { |
---|
| 242 | #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 | #ifdef debugphoton |
---|
| 260 | G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n" |
---|
| 261 | << "-----------------------------------------------------------------------\n" |
---|
| 262 | << theExcitedNucleus << '\n' |
---|
| 263 | << "-----------------------------------------------------------------------\n"; |
---|
| 264 | #endif |
---|
| 265 | G4double GammaEnergy = (*j)->GetExcitationEnergy(); |
---|
| 266 | G4double cosTheta = 1. - 2. * G4UniformRand(); |
---|
| 267 | G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); |
---|
| 268 | G4double phi = twopi * G4UniformRand(); |
---|
| 269 | G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi), |
---|
| 270 | GammaEnergy * sinTheta * std::sin(phi), |
---|
| 271 | GammaEnergy * cosTheta ); |
---|
| 272 | G4LorentzVector Gamma4P(GammaP,GammaEnergy); |
---|
| 273 | G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition()); |
---|
| 274 | |
---|
| 275 | |
---|
| 276 | |
---|
| 277 | G4double Mass = (*j)->GetGroundStateMass(); |
---|
| 278 | G4ThreeVector ResidualP((*j)->GetMomentum().vect() - GammaP); |
---|
| 279 | G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass); |
---|
| 280 | G4LorentzVector Residual4P(ResidualP,ResidualE); |
---|
| 281 | (*j)->SetMomentum(Residual4P); |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | |
---|
| 285 | #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); |
---|
| 293 | #ifdef debugphoton |
---|
| 294 | G4cout << "Emmited photon:\n" |
---|
| 295 | << theResultList.back() << '\n' |
---|
| 296 | << "Residual nucleus after photon emission:\n" |
---|
| 297 | << *(*j) << '\n' |
---|
| 298 | << "-----------------------------------------------------------------------\n"; |
---|
| 299 | #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); |
---|
| 305 | } |
---|
| 306 | // G4cout << " Inside loop list: " << (*j)->GetExcitationEnergy() << " size: " << theFinalResultList.size() << G4endl; |
---|
| 307 | j++; |
---|
| 308 | } |
---|
| 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 | |
---|
| 327 | #ifdef debug |
---|
| 328 | CheckConservation(theInitialState,theResult); |
---|
| 329 | #endif |
---|
| 330 | // Change G4FragmentVector by G4DynamicParticle |
---|
| 331 | return Transform(theResult); |
---|
| 332 | } |
---|
| 333 | |
---|
| 334 | G4ReactionProductVector * |
---|
| 335 | G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const |
---|
| 336 | { |
---|
| 337 | if (theFragmentVector == 0) return 0; |
---|
| 338 | |
---|
| 339 | // Conversion from G4FragmentVector to G4ReactionProductVector |
---|
| 340 | G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition(); |
---|
| 341 | G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition(); |
---|
| 342 | G4ParticleDefinition *theProton = G4Proton::ProtonDefinition(); |
---|
| 343 | G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition(); |
---|
| 344 | G4ParticleDefinition *theTriton = G4Triton::TritonDefinition(); |
---|
| 345 | G4ParticleDefinition *theHelium3 = G4He3::He3Definition(); |
---|
| 346 | G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition(); |
---|
| 347 | G4ParticleDefinition *theKindOfFragment = 0; |
---|
| 348 | theNeutron->SetVerboseLevel(2); |
---|
| 349 | G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector; |
---|
| 350 | G4int theFragmentA, theFragmentZ; |
---|
| 351 | G4LorentzVector theFragmentMomentum; |
---|
| 352 | |
---|
| 353 | G4FragmentVector::iterator i; |
---|
| 354 | for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) { |
---|
| 355 | // std::cout << (*i) <<'\n'; |
---|
| 356 | theFragmentA = static_cast<G4int>((*i)->GetA()); |
---|
| 357 | theFragmentZ = static_cast<G4int>((*i)->GetZ()); |
---|
| 358 | theFragmentMomentum = (*i)->GetMomentum(); |
---|
| 359 | theKindOfFragment = 0; |
---|
| 360 | if (theFragmentA == 0 && theFragmentZ == 0) { // photon |
---|
| 361 | theKindOfFragment = theGamma; |
---|
| 362 | } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron |
---|
| 363 | theKindOfFragment = theNeutron; |
---|
| 364 | } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton |
---|
| 365 | theKindOfFragment = theProton; |
---|
| 366 | } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron |
---|
| 367 | theKindOfFragment = theDeuteron; |
---|
| 368 | } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton |
---|
| 369 | theKindOfFragment = theTriton; |
---|
| 370 | } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3 |
---|
| 371 | theKindOfFragment = theHelium3; |
---|
| 372 | } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha |
---|
| 373 | theKindOfFragment = theAlpha; |
---|
| 374 | } else { |
---|
| 375 | theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ); |
---|
| 376 | } |
---|
| 377 | if (theKindOfFragment != 0) |
---|
| 378 | { |
---|
| 379 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
---|
| 380 | theNew->SetMomentum(theFragmentMomentum.vect()); |
---|
| 381 | theNew->SetTotalEnergy(theFragmentMomentum.e()); |
---|
| 382 | theNew->SetFormationTime((*i)->GetCreationTime()); |
---|
| 383 | #ifdef PRECOMPOUND_TEST |
---|
| 384 | theNew->SetCreatorModel((*i)->GetCreatorModel()); |
---|
| 385 | #endif |
---|
| 386 | theReactionProductVector->push_back(theNew); |
---|
| 387 | } |
---|
| 388 | } |
---|
| 389 | if (theFragmentVector != 0) |
---|
| 390 | { |
---|
| 391 | std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment()); |
---|
| 392 | delete theFragmentVector; |
---|
| 393 | } |
---|
| 394 | G4ReactionProductVector::iterator debugit; |
---|
| 395 | for(debugit=theReactionProductVector->begin(); |
---|
| 396 | debugit!=theReactionProductVector->end(); debugit++) |
---|
| 397 | { |
---|
| 398 | if((*debugit)->GetTotalEnergy()<1.*eV) |
---|
| 399 | { |
---|
| 400 | if(getenv("G4DebugPhotonevaporationData")) |
---|
| 401 | { |
---|
| 402 | G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl; |
---|
| 403 | G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = " |
---|
| 404 | << (*debugit)->GetTotalEnergy()/MeV << "MeV" |
---|
| 405 | << G4endl; |
---|
| 406 | } |
---|
| 407 | delete (*debugit); |
---|
| 408 | *debugit = 0; |
---|
| 409 | } |
---|
| 410 | } |
---|
| 411 | G4ReactionProduct* tmpPtr=0; |
---|
| 412 | theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(), |
---|
| 413 | theReactionProductVector->end(), |
---|
| 414 | std::bind2nd(std::equal_to<G4ReactionProduct*>(), |
---|
| 415 | tmpPtr)), |
---|
| 416 | theReactionProductVector->end()); |
---|
| 417 | return theReactionProductVector; |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | #ifdef debug |
---|
| 422 | void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState, |
---|
| 423 | G4FragmentVector * Result) const |
---|
| 424 | { |
---|
| 425 | G4double ProductsEnergy =0; |
---|
| 426 | G4ThreeVector ProductsMomentum; |
---|
| 427 | G4int ProductsA = 0; |
---|
| 428 | G4int ProductsZ = 0; |
---|
| 429 | G4FragmentVector::iterator h; |
---|
| 430 | for (h = Result->begin(); h != Result->end(); h++) { |
---|
| 431 | G4LorentzVector tmp = (*h)->GetMomentum(); |
---|
| 432 | ProductsEnergy += tmp.e(); |
---|
| 433 | ProductsMomentum += tmp.vect(); |
---|
| 434 | ProductsA += static_cast<G4int>((*h)->GetA()); |
---|
| 435 | ProductsZ += static_cast<G4int>((*h)->GetZ()); |
---|
| 436 | } |
---|
| 437 | |
---|
| 438 | if (ProductsA != theInitialState.GetA()) { |
---|
| 439 | G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 440 | G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments" |
---|
| 441 | << G4endl; |
---|
| 442 | G4cout << "Initial A = " << theInitialState.GetA() |
---|
| 443 | << " Fragments A = " << ProductsA << " Diference --> " |
---|
| 444 | << theInitialState.GetA() - ProductsA << G4endl; |
---|
| 445 | } |
---|
| 446 | if (ProductsZ != theInitialState.GetZ()) { |
---|
| 447 | G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 448 | G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments" |
---|
| 449 | << G4endl; |
---|
| 450 | G4cout << "Initial Z = " << theInitialState.GetZ() |
---|
| 451 | << " Fragments Z = " << ProductsZ << " Diference --> " |
---|
| 452 | << theInitialState.GetZ() - ProductsZ << G4endl; |
---|
| 453 | } |
---|
| 454 | if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) { |
---|
| 455 | G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 456 | G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments" |
---|
| 457 | << G4endl; |
---|
| 458 | G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV" |
---|
| 459 | << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> " |
---|
| 460 | << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl; |
---|
| 461 | } |
---|
| 462 | if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || |
---|
| 463 | std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV || |
---|
| 464 | std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) { |
---|
| 465 | G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 466 | G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments" |
---|
| 467 | << G4endl; |
---|
| 468 | G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV" |
---|
| 469 | << " Fragments P = " << ProductsMomentum << " MeV Diference --> " |
---|
| 470 | << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl; |
---|
| 471 | } |
---|
| 472 | return; |
---|
| 473 | } |
---|
| 474 | #endif |
---|
| 475 | |
---|
| 476 | |
---|
| 477 | |
---|
| 478 | |
---|