[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) |
---|
[962] | 29 | // |
---|
| 30 | // |
---|
[1196] | 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 | // |
---|
[962] | 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 |
---|
[819] | 42 | // Modif (30 June 1998) by V. Lara: |
---|
| 43 | // -Modified the Transform method for use G4ParticleTable and |
---|
| 44 | // therefore G4IonTable. It makes possible to convert all kind |
---|
| 45 | // of fragments (G4Fragment) produced in deexcitation to |
---|
| 46 | // G4DynamicParticle |
---|
| 47 | // -It uses default algorithms for: |
---|
[962] | 48 | // Evaporation: G4Evaporation |
---|
| 49 | // MultiFragmentation: G4StatMF |
---|
| 50 | // Fermi Breakup model: G4FermiBreakUp |
---|
| 51 | // |
---|
[819] | 52 | |
---|
| 53 | #include "G4ExcitationHandler.hh" |
---|
| 54 | #include <list> |
---|
| 55 | |
---|
| 56 | //#define debugphoton |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | G4ExcitationHandler::G4ExcitationHandler(): |
---|
[1196] | 60 | // JMQ 160909 Fermi BreakUp & MultiFrag are on by default |
---|
| 61 | // This is needed for activation of such models when G4BinaryLightIonReaction is used |
---|
| 62 | // since no interface (for external activation via macro input file) is still available in this case. |
---|
| 63 | //maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV), |
---|
[819] | 64 | maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV), |
---|
| 65 | MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true), |
---|
[962] | 66 | MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false) |
---|
[819] | 67 | { |
---|
| 68 | theTableOfParticles = G4ParticleTable::GetParticleTable(); |
---|
| 69 | |
---|
| 70 | theEvaporation = new G4Evaporation; |
---|
| 71 | theMultiFragmentation = new G4StatMF; |
---|
| 72 | theFermiModel = new G4FermiBreakUp; |
---|
| 73 | thePhotonEvaporation = new G4PhotonEvaporation; |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &) |
---|
| 77 | { |
---|
| 78 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! "); |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | G4ExcitationHandler::~G4ExcitationHandler() |
---|
| 83 | { |
---|
| 84 | if (MyOwnEvaporationClass) delete theEvaporation; |
---|
| 85 | if (MyOwnMultiFragmentationClass) delete theMultiFragmentation; |
---|
| 86 | if (MyOwnFermiBreakUpClass) delete theFermiModel; |
---|
| 87 | if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation; |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &) |
---|
| 92 | { |
---|
| 93 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! "); |
---|
[1196] | 94 | |
---|
[819] | 95 | return *this; |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const |
---|
| 100 | { |
---|
| 101 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! "); |
---|
| 102 | return false; |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const |
---|
| 106 | { |
---|
| 107 | throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! "); |
---|
| 108 | return true; |
---|
| 109 | } |
---|
| 110 | |
---|
[962] | 111 | //////////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 112 | /// 25/07/08 16:45 Proposed by MAC //////////////////////////////////////////////////////////// |
---|
| 113 | //////////////////////////////////////////////////////////////////////////////////////////////// |
---|
[819] | 114 | |
---|
[962] | 115 | G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const |
---|
| 116 | { |
---|
| 117 | //for inverse cross section choice |
---|
| 118 | theEvaporation->SetOPTxs(OPTxs); |
---|
| 119 | //for the choice of superimposed Coulomb Barrier for inverse cross sections |
---|
| 120 | theEvaporation->UseSICB(useSICB); |
---|
[1196] | 121 | |
---|
[962] | 122 | // Pointer which will be used to return the final production vector |
---|
| 123 | G4FragmentVector * theResult = 0; |
---|
[1196] | 124 | |
---|
[962] | 125 | // Variables existing until end of method |
---|
[1196] | 126 | G4Fragment * theInitialStatePtr = const_cast<G4Fragment*>(&theInitialState); |
---|
| 127 | // G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState); |
---|
[962] | 128 | G4Fragment theExcitedNucleus; // object to be passed in BreakItUp methods |
---|
| 129 | G4FragmentVector * theTempResult = 0; // pointer which receives temporal results |
---|
| 130 | std::list<G4Fragment*> theEvapList; // list to apply Evaporation, SMF or Fermi Break-Up |
---|
| 131 | std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation |
---|
| 132 | std::list<G4Fragment*> theFinalStableList; // list to store final result |
---|
| 133 | std::list<G4Fragment*>::iterator iList; |
---|
[1196] | 134 | // |
---|
| 135 | |
---|
| 136 | |
---|
[962] | 137 | // Variables to describe the excited configuration |
---|
[819] | 138 | G4double exEnergy = theInitialState.GetExcitationEnergy(); |
---|
[962] | 139 | G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 ); |
---|
| 140 | G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 ); |
---|
[1196] | 141 | |
---|
| 142 | // JMQ 150909: first step in de-excitation chain (SMM will be used only here) |
---|
| 143 | |
---|
[962] | 144 | // In case A <= 4 the fragment will not perform any nucleon emission |
---|
| 145 | if (A <= 4) |
---|
[819] | 146 | { |
---|
[962] | 147 | // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later |
---|
[1196] | 148 | |
---|
[962] | 149 | theEvapStableList.push_back( theInitialStatePtr ); |
---|
| 150 | } |
---|
[1196] | 151 | |
---|
[962] | 152 | else // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation |
---|
| 153 | { |
---|
[1196] | 154 | |
---|
| 155 | // JMQ 150909: first step in de-excitation is treated separately |
---|
| 156 | // Fragments after the first step are stored in theEvapList |
---|
| 157 | // Statistical Multifragmentation will take place (just in case) only here |
---|
| 158 | // |
---|
| 159 | // Test applicability |
---|
| 160 | // Initial State De-Excitation |
---|
| 161 | if(A<GetMaxA()&&Z<GetMaxZ()) |
---|
| 162 | { |
---|
| 163 | theTempResult = theFermiModel->BreakItUp(theInitialState); |
---|
| 164 | } |
---|
| 165 | else if (exEnergy>GetMinE()*A) |
---|
| 166 | { |
---|
| 167 | theTempResult = theMultiFragmentation->BreakItUp(theInitialState); |
---|
| 168 | } |
---|
| 169 | else |
---|
| 170 | { |
---|
| 171 | theTempResult = theEvaporation->BreakItUp(theInitialState); |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | |
---|
[962] | 175 | // Store original state in theEvapList |
---|
[1196] | 176 | G4FragmentVector::iterator j; |
---|
| 177 | for (j = theTempResult->begin(); j != theTempResult->end(); j++) |
---|
| 178 | { |
---|
| 179 | theEvapList.push_back(*j); |
---|
| 180 | // This memory release works with a list but not with a G4FragmentVector |
---|
| 181 | // delete (*j); |
---|
| 182 | } |
---|
| 183 | delete theTempResult; |
---|
| 184 | //theTempResult->clear(); |
---|
| 185 | // |
---|
| 186 | // JMQ 150909: Further steps in de-excitation chain follow .. |
---|
| 187 | |
---|
[962] | 188 | // ------------------------------ |
---|
| 189 | // De-excitation loop |
---|
| 190 | // ------------------------------ |
---|
[1196] | 191 | |
---|
[962] | 192 | for (iList = theEvapList.begin(); iList != theEvapList.end(); iList++) |
---|
| 193 | { |
---|
[1196] | 194 | A = static_cast<G4int>((*iList)->GetA()+0.5); // +0.5 to avoid bad truncation |
---|
| 195 | Z = static_cast<G4int>((*iList)->GetZ()+0.5); |
---|
| 196 | |
---|
| 197 | // In case A <= 4 the fragment will not perform any nucleon emission |
---|
| 198 | if (A <= 4) |
---|
[962] | 199 | { |
---|
[1196] | 200 | // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later |
---|
| 201 | |
---|
| 202 | theEvapStableList.push_back(*iList ); |
---|
| 203 | } |
---|
| 204 | |
---|
| 205 | else // If A > 4 we try to apply theFermiModel or theEvaporation |
---|
| 206 | { |
---|
| 207 | exEnergy = (*iList)->GetExcitationEnergy(); |
---|
| 208 | |
---|
| 209 | if (exEnergy > 0.0) |
---|
[962] | 210 | { |
---|
[1196] | 211 | // Check conditions for each model |
---|
| 212 | theExcitedNucleus = *(*iList); |
---|
| 213 | |
---|
| 214 | if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up |
---|
[962] | 215 | { |
---|
[1196] | 216 | theTempResult = theFermiModel->BreakItUp(theExcitedNucleus); |
---|
[962] | 217 | } |
---|
[1196] | 218 | else // apply Evaporation in another case |
---|
[962] | 219 | { |
---|
[1196] | 220 | theTempResult = theEvaporation->BreakItUp(theExcitedNucleus); |
---|
[962] | 221 | } |
---|
[1196] | 222 | |
---|
| 223 | // New configuration is stored in theTempResult, so we can free |
---|
| 224 | // the memory where the previous configuration is |
---|
| 225 | |
---|
| 226 | delete (*iList); |
---|
| 227 | |
---|
| 228 | // And now the theTempResult->size() tells us if the configuration has changed |
---|
| 229 | |
---|
| 230 | if ( theTempResult->size() > 1 ) |
---|
| 231 | { |
---|
| 232 | // push_back the result to the end of theEvapList (this same list) |
---|
| 233 | |
---|
| 234 | for (G4FragmentVector::iterator j = theTempResult->begin(); |
---|
| 235 | j != theTempResult->end(); j++) |
---|
| 236 | { |
---|
| 237 | theEvapList.push_back(*j); |
---|
| 238 | } |
---|
| 239 | } |
---|
| 240 | else |
---|
| 241 | { |
---|
| 242 | // push_back the result to theEvapStableList, because |
---|
| 243 | // is still excited, but cannot emmit more nucleons |
---|
| 244 | |
---|
| 245 | for (G4FragmentVector::iterator j = theTempResult->begin(); |
---|
| 246 | j != theTempResult->end(); j++) |
---|
| 247 | { |
---|
| 248 | theEvapStableList.push_back(*j); |
---|
| 249 | } |
---|
| 250 | } |
---|
| 251 | |
---|
| 252 | // after working with theTempResult, clear and delete it |
---|
| 253 | theTempResult->clear(); |
---|
| 254 | delete theTempResult; |
---|
| 255 | |
---|
[962] | 256 | } |
---|
[1196] | 257 | else // exEnergy = 0.0 |
---|
| 258 | { |
---|
| 259 | // if this fragment is at ground state, |
---|
| 260 | // store it in theFinalStableList |
---|
| 261 | |
---|
| 262 | theFinalStableList.push_back(*iList); |
---|
| 263 | |
---|
| 264 | } // endif (exEnergy > 0.0) |
---|
| 265 | } // endif (A <=4) |
---|
[962] | 266 | } // end of the loop over theEvapList |
---|
[1196] | 267 | |
---|
[962] | 268 | theEvapList.clear(); // clear all the list and do not free memory pointed by |
---|
| 269 | // each element because this have been done before! |
---|
[819] | 270 | } |
---|
[1196] | 271 | |
---|
[819] | 272 | // Now we try to deexcite by means of PhotonEvaporation those fragments |
---|
[962] | 273 | // which are still excited. |
---|
[1196] | 274 | |
---|
| 275 | |
---|
[962] | 276 | // ----------------------- |
---|
| 277 | // Photon-Evaporation loop |
---|
| 278 | // ----------------------- |
---|
[1196] | 279 | |
---|
[962] | 280 | for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); iList++) |
---|
[819] | 281 | { |
---|
[962] | 282 | A = static_cast<G4int>((*iList)->GetA()+0.5); |
---|
| 283 | exEnergy = (*iList)->GetExcitationEnergy(); |
---|
[1196] | 284 | |
---|
[962] | 285 | if ( A > 1 && exEnergy > 0.1*eV ) // if so, photon-evaporation is applied |
---|
| 286 | { |
---|
| 287 | theExcitedNucleus = *(*iList); |
---|
| 288 | theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus); |
---|
[1196] | 289 | |
---|
| 290 | |
---|
[962] | 291 | // if there is a gamma emission then |
---|
| 292 | if (theTempResult->size() > 1) |
---|
| 293 | { |
---|
| 294 | // first free the memory occupied by the previous state |
---|
| 295 | delete (*iList); |
---|
[1196] | 296 | |
---|
[962] | 297 | // and now add the final state from gamma emission to the end of |
---|
| 298 | // theEvapStableList |
---|
| 299 | for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin(); |
---|
[819] | 300 | ri != theTempResult->rend(); ++ri) |
---|
[962] | 301 | // reversed is applied in order to have residual nucleus in first position |
---|
| 302 | { |
---|
[1196] | 303 | |
---|
[819] | 304 | #ifdef PRECOMPOUND_TEST |
---|
[962] | 305 | if ((*ri)->GetA() == 0) |
---|
| 306 | (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation")); |
---|
| 307 | else |
---|
| 308 | (*ri)->SetCreatorModel(G4String("ResidualNucleus")); |
---|
[819] | 309 | #endif |
---|
[1196] | 310 | |
---|
[962] | 311 | theEvapStableList.push_back(*ri); |
---|
| 312 | } |
---|
[1196] | 313 | |
---|
[962] | 314 | // now we clean and remove the temporal vector |
---|
| 315 | theTempResult->clear(); |
---|
| 316 | delete theTempResult; |
---|
[1196] | 317 | |
---|
[962] | 318 | } |
---|
| 319 | else // if theTempResult->size() = 1 |
---|
| 320 | { |
---|
| 321 | // if there is not any gamma emission from this excited fragment |
---|
| 322 | // we have to emmit a gamma which forces the deexcitation |
---|
[1196] | 323 | |
---|
[962] | 324 | // First I clean completely theTempResult |
---|
[1196] | 325 | |
---|
[962] | 326 | for (G4FragmentVector::iterator j = theTempResult->begin(); |
---|
| 327 | j != theTempResult->end(); j++) |
---|
| 328 | { |
---|
| 329 | delete (*j); |
---|
| 330 | } |
---|
[1196] | 331 | |
---|
[962] | 332 | theTempResult->clear(); |
---|
| 333 | delete theTempResult; |
---|
[1196] | 334 | |
---|
[819] | 335 | #ifdef debugphoton |
---|
| 336 | G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n" |
---|
| 337 | << "-----------------------------------------------------------------------\n" |
---|
| 338 | << theExcitedNucleus << '\n' |
---|
| 339 | << "-----------------------------------------------------------------------\n"; |
---|
| 340 | #endif |
---|
[1196] | 341 | |
---|
[962] | 342 | // Let's create a G4Fragment pointer representing the gamma emmited |
---|
| 343 | G4double GammaEnergy = (*iList)->GetExcitationEnergy(); |
---|
| 344 | G4double cosTheta = 1. - 2. * G4UniformRand(); |
---|
| 345 | G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); |
---|
[819] | 346 | G4double phi = twopi * G4UniformRand(); |
---|
| 347 | G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi), |
---|
| 348 | GammaEnergy * sinTheta * std::sin(phi), |
---|
| 349 | GammaEnergy * cosTheta ); |
---|
| 350 | G4LorentzVector Gamma4P(GammaP,GammaEnergy); |
---|
| 351 | G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition()); |
---|
[1196] | 352 | |
---|
[962] | 353 | // And now we update momentum and energy for the nucleus |
---|
| 354 | G4double Mass = (*iList)->GetGroundStateMass(); |
---|
| 355 | G4ThreeVector ResidualP((*iList)->GetMomentum().vect() - GammaP); |
---|
[819] | 356 | G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass); |
---|
| 357 | G4LorentzVector Residual4P(ResidualP,ResidualE); |
---|
[962] | 358 | (*iList)->SetMomentum(Residual4P); // Now this fragment has been deexcited! |
---|
[1196] | 359 | |
---|
[962] | 360 | // we store the deexcited fragment in theFinalStableList |
---|
| 361 | theFinalStableList.push_back(*iList); |
---|
[1196] | 362 | |
---|
[819] | 363 | #ifdef PRECOMPOUND_TEST |
---|
[962] | 364 | theHandlerPhoton->SetCreatorModel("G4ExcitationHandler"); |
---|
[819] | 365 | #endif |
---|
[1196] | 366 | |
---|
[962] | 367 | // Finally, we add theHandlerPhoton to theFinalStableList |
---|
| 368 | theFinalStableList.push_back(theHandlerPhoton); |
---|
[1196] | 369 | |
---|
[819] | 370 | #ifdef debugphoton |
---|
| 371 | G4cout << "Emmited photon:\n" |
---|
[962] | 372 | << theFinalStableList.back() << '\n' |
---|
[819] | 373 | << "Residual nucleus after photon emission:\n" |
---|
[962] | 374 | << *(*iList) << '\n' |
---|
[819] | 375 | << "-----------------------------------------------------------------------\n"; |
---|
| 376 | #endif |
---|
[1196] | 377 | |
---|
[962] | 378 | } |
---|
[819] | 379 | } |
---|
[962] | 380 | else // case of a nucleon, gamma or very small excitation energy |
---|
| 381 | { |
---|
| 382 | // we don't have to do anything, just store the fragment in theFinalStableList |
---|
[1196] | 383 | |
---|
[962] | 384 | theFinalStableList.push_back(*iList); |
---|
[1196] | 385 | |
---|
[962] | 386 | } // A > 1 && exEnergy > 0.1*eV |
---|
[1196] | 387 | |
---|
[962] | 388 | } // end of photon-evaporation loop |
---|
[1196] | 389 | |
---|
| 390 | |
---|
[962] | 391 | // The deexcitation from fragments inside theEvapStableList has been finished, so... |
---|
| 392 | theEvapStableList.clear(); |
---|
[1196] | 393 | |
---|
[962] | 394 | // Now the final state is in theFinalStableList, and we have to send it to theResult vector |
---|
[1196] | 395 | |
---|
[962] | 396 | theResult = new G4FragmentVector; |
---|
| 397 | theResult->reserve( theFinalStableList.size() * sizeof(G4Fragment*) ); |
---|
| 398 | // We reserve enough memory to optimise the storing speed |
---|
[1196] | 399 | |
---|
[962] | 400 | for (iList = theFinalStableList.begin(); iList != theFinalStableList.end(); iList++) |
---|
[819] | 401 | { |
---|
[962] | 402 | theResult->push_back(*iList); |
---|
[819] | 403 | } |
---|
[1196] | 404 | |
---|
[962] | 405 | // After storing the final state , we can clear theFinalStableList |
---|
| 406 | theFinalStableList.clear(); |
---|
[1196] | 407 | |
---|
[819] | 408 | #ifdef debug |
---|
| 409 | CheckConservation(theInitialState,theResult); |
---|
| 410 | #endif |
---|
[1196] | 411 | |
---|
| 412 | |
---|
[962] | 413 | // Change G4FragmentVector* to G4ReactionProductVector* |
---|
[819] | 414 | return Transform(theResult); |
---|
| 415 | } |
---|
| 416 | |
---|
[962] | 417 | |
---|
| 418 | |
---|
[819] | 419 | G4ReactionProductVector * |
---|
| 420 | G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const |
---|
| 421 | { |
---|
| 422 | if (theFragmentVector == 0) return 0; |
---|
| 423 | |
---|
| 424 | // Conversion from G4FragmentVector to G4ReactionProductVector |
---|
| 425 | G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition(); |
---|
| 426 | G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition(); |
---|
| 427 | G4ParticleDefinition *theProton = G4Proton::ProtonDefinition(); |
---|
| 428 | G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition(); |
---|
| 429 | G4ParticleDefinition *theTriton = G4Triton::TritonDefinition(); |
---|
| 430 | G4ParticleDefinition *theHelium3 = G4He3::He3Definition(); |
---|
| 431 | G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition(); |
---|
| 432 | G4ParticleDefinition *theKindOfFragment = 0; |
---|
| 433 | theNeutron->SetVerboseLevel(2); |
---|
| 434 | G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector; |
---|
[962] | 435 | |
---|
| 436 | // MAC (24/07/08) |
---|
| 437 | // To optimise the storing speed, we reserve space in memory for the vector |
---|
| 438 | theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) ); |
---|
| 439 | |
---|
[819] | 440 | G4int theFragmentA, theFragmentZ; |
---|
| 441 | G4LorentzVector theFragmentMomentum; |
---|
| 442 | |
---|
| 443 | G4FragmentVector::iterator i; |
---|
| 444 | for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) { |
---|
| 445 | // std::cout << (*i) <<'\n'; |
---|
| 446 | theFragmentA = static_cast<G4int>((*i)->GetA()); |
---|
| 447 | theFragmentZ = static_cast<G4int>((*i)->GetZ()); |
---|
| 448 | theFragmentMomentum = (*i)->GetMomentum(); |
---|
| 449 | theKindOfFragment = 0; |
---|
| 450 | if (theFragmentA == 0 && theFragmentZ == 0) { // photon |
---|
| 451 | theKindOfFragment = theGamma; |
---|
| 452 | } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron |
---|
| 453 | theKindOfFragment = theNeutron; |
---|
| 454 | } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton |
---|
| 455 | theKindOfFragment = theProton; |
---|
| 456 | } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron |
---|
| 457 | theKindOfFragment = theDeuteron; |
---|
| 458 | } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton |
---|
| 459 | theKindOfFragment = theTriton; |
---|
| 460 | } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3 |
---|
| 461 | theKindOfFragment = theHelium3; |
---|
| 462 | } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha |
---|
| 463 | theKindOfFragment = theAlpha; |
---|
| 464 | } else { |
---|
| 465 | theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ); |
---|
| 466 | } |
---|
| 467 | if (theKindOfFragment != 0) |
---|
| 468 | { |
---|
| 469 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
---|
| 470 | theNew->SetMomentum(theFragmentMomentum.vect()); |
---|
| 471 | theNew->SetTotalEnergy(theFragmentMomentum.e()); |
---|
| 472 | theNew->SetFormationTime((*i)->GetCreationTime()); |
---|
| 473 | #ifdef PRECOMPOUND_TEST |
---|
| 474 | theNew->SetCreatorModel((*i)->GetCreatorModel()); |
---|
| 475 | #endif |
---|
| 476 | theReactionProductVector->push_back(theNew); |
---|
| 477 | } |
---|
| 478 | } |
---|
| 479 | if (theFragmentVector != 0) |
---|
| 480 | { |
---|
| 481 | std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment()); |
---|
| 482 | delete theFragmentVector; |
---|
| 483 | } |
---|
| 484 | G4ReactionProductVector::iterator debugit; |
---|
| 485 | for(debugit=theReactionProductVector->begin(); |
---|
| 486 | debugit!=theReactionProductVector->end(); debugit++) |
---|
| 487 | { |
---|
| 488 | if((*debugit)->GetTotalEnergy()<1.*eV) |
---|
| 489 | { |
---|
| 490 | if(getenv("G4DebugPhotonevaporationData")) |
---|
| 491 | { |
---|
| 492 | G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl; |
---|
| 493 | G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = " |
---|
| 494 | << (*debugit)->GetTotalEnergy()/MeV << "MeV" |
---|
| 495 | << G4endl; |
---|
| 496 | } |
---|
| 497 | delete (*debugit); |
---|
| 498 | *debugit = 0; |
---|
| 499 | } |
---|
| 500 | } |
---|
| 501 | G4ReactionProduct* tmpPtr=0; |
---|
| 502 | theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(), |
---|
| 503 | theReactionProductVector->end(), |
---|
| 504 | std::bind2nd(std::equal_to<G4ReactionProduct*>(), |
---|
| 505 | tmpPtr)), |
---|
| 506 | theReactionProductVector->end()); |
---|
| 507 | return theReactionProductVector; |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | |
---|
| 511 | #ifdef debug |
---|
| 512 | void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState, |
---|
| 513 | G4FragmentVector * Result) const |
---|
| 514 | { |
---|
| 515 | G4double ProductsEnergy =0; |
---|
| 516 | G4ThreeVector ProductsMomentum; |
---|
| 517 | G4int ProductsA = 0; |
---|
| 518 | G4int ProductsZ = 0; |
---|
| 519 | G4FragmentVector::iterator h; |
---|
| 520 | for (h = Result->begin(); h != Result->end(); h++) { |
---|
| 521 | G4LorentzVector tmp = (*h)->GetMomentum(); |
---|
| 522 | ProductsEnergy += tmp.e(); |
---|
| 523 | ProductsMomentum += tmp.vect(); |
---|
| 524 | ProductsA += static_cast<G4int>((*h)->GetA()); |
---|
| 525 | ProductsZ += static_cast<G4int>((*h)->GetZ()); |
---|
| 526 | } |
---|
| 527 | |
---|
| 528 | if (ProductsA != theInitialState.GetA()) { |
---|
| 529 | G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 530 | G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments" |
---|
| 531 | << G4endl; |
---|
| 532 | G4cout << "Initial A = " << theInitialState.GetA() |
---|
| 533 | << " Fragments A = " << ProductsA << " Diference --> " |
---|
| 534 | << theInitialState.GetA() - ProductsA << G4endl; |
---|
| 535 | } |
---|
| 536 | if (ProductsZ != theInitialState.GetZ()) { |
---|
| 537 | G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 538 | G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments" |
---|
| 539 | << G4endl; |
---|
| 540 | G4cout << "Initial Z = " << theInitialState.GetZ() |
---|
| 541 | << " Fragments Z = " << ProductsZ << " Diference --> " |
---|
| 542 | << theInitialState.GetZ() - ProductsZ << G4endl; |
---|
| 543 | } |
---|
| 544 | if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) { |
---|
| 545 | G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 546 | G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments" |
---|
| 547 | << G4endl; |
---|
| 548 | G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV" |
---|
| 549 | << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> " |
---|
| 550 | << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl; |
---|
| 551 | } |
---|
| 552 | if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || |
---|
| 553 | std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV || |
---|
| 554 | std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) { |
---|
| 555 | G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 556 | G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments" |
---|
| 557 | << G4endl; |
---|
| 558 | G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV" |
---|
| 559 | << " Fragments P = " << ProductsMomentum << " MeV Diference --> " |
---|
| 560 | << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl; |
---|
| 561 | } |
---|
| 562 | return; |
---|
| 563 | } |
---|
| 564 | #endif |
---|
| 565 | |
---|
| 566 | |
---|
| 567 | |
---|
| 568 | |
---|