[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 | // |
---|
[1315] | 26 | // $Id: G4ExcitationHandler.cc,v 1.39 2010/05/18 15:46:11 vnivanch Exp $ |
---|
[1340] | 27 | // GEANT4 tag $Name: geant4-09-03-ref-09 $ |
---|
[819] | 28 | // |
---|
| 29 | // Hadronic Process: Nuclear De-excitations |
---|
| 30 | // by V. Lara (May 1998) |
---|
[962] | 31 | // |
---|
| 32 | // |
---|
[1315] | 33 | // Modified: |
---|
| 34 | // (30 June 1998) by V. Lara: |
---|
[819] | 35 | // -Modified the Transform method for use G4ParticleTable and |
---|
| 36 | // therefore G4IonTable. It makes possible to convert all kind |
---|
| 37 | // of fragments (G4Fragment) produced in deexcitation to |
---|
| 38 | // G4DynamicParticle |
---|
| 39 | // -It uses default algorithms for: |
---|
[962] | 40 | // Evaporation: G4Evaporation |
---|
| 41 | // MultiFragmentation: G4StatMF |
---|
| 42 | // Fermi Breakup model: G4FermiBreakUp |
---|
[1315] | 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 |
---|
[962] | 58 | // |
---|
[819] | 59 | |
---|
| 60 | #include "G4ExcitationHandler.hh" |
---|
[1228] | 61 | #include "globals.hh" |
---|
| 62 | #include "G4LorentzVector.hh" |
---|
[1315] | 63 | #include "G4NistManager.hh" |
---|
[819] | 64 | #include <list> |
---|
| 65 | |
---|
| 66 | G4ExcitationHandler::G4ExcitationHandler(): |
---|
[1196] | 67 | // JMQ 160909 Fermi BreakUp & MultiFrag are on by default |
---|
| 68 | // This is needed for activation of such models when G4BinaryLightIonReaction is used |
---|
[1315] | 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), |
---|
[819] | 74 | MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true), |
---|
[962] | 75 | MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false) |
---|
[819] | 76 | { |
---|
| 77 | theTableOfParticles = G4ParticleTable::GetParticleTable(); |
---|
| 78 | |
---|
| 79 | theEvaporation = new G4Evaporation; |
---|
| 80 | theMultiFragmentation = new G4StatMF; |
---|
| 81 | theFermiModel = new G4FermiBreakUp; |
---|
| 82 | thePhotonEvaporation = new G4PhotonEvaporation; |
---|
[1315] | 83 | SetParameters(); |
---|
[819] | 84 | } |
---|
| 85 | |
---|
| 86 | G4ExcitationHandler::~G4ExcitationHandler() |
---|
| 87 | { |
---|
| 88 | if (MyOwnEvaporationClass) delete theEvaporation; |
---|
| 89 | if (MyOwnMultiFragmentationClass) delete theMultiFragmentation; |
---|
| 90 | if (MyOwnFermiBreakUpClass) delete theFermiModel; |
---|
| 91 | if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation; |
---|
| 92 | } |
---|
| 93 | |
---|
[1315] | 94 | void G4ExcitationHandler::SetParameters() |
---|
[819] | 95 | { |
---|
[1315] | 96 | //for inverse cross section choice |
---|
| 97 | theEvaporation->SetOPTxs(OPTxs); |
---|
| 98 | //for the choice of superimposed Coulomb Barrier for inverse cross sections |
---|
| 99 | theEvaporation->UseSICB(useSICB); |
---|
| 100 | theEvaporation->Initialise(); |
---|
[819] | 101 | } |
---|
| 102 | |
---|
[962] | 103 | G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const |
---|
| 104 | { |
---|
[1196] | 105 | |
---|
[962] | 106 | // Variables existing until end of method |
---|
[1228] | 107 | G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState); |
---|
[1315] | 108 | |
---|
[962] | 109 | G4FragmentVector * theTempResult = 0; // pointer which receives temporal results |
---|
[1315] | 110 | std::list<G4Fragment*> theEvapList; // list to apply Evaporation or Fermi Break-Up |
---|
[962] | 111 | std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation |
---|
[1228] | 112 | std::list<G4Fragment*> theResults; // list to store final result |
---|
[962] | 113 | std::list<G4Fragment*>::iterator iList; |
---|
[1196] | 114 | // |
---|
[1315] | 115 | //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl; |
---|
[1228] | 116 | //G4cout << theInitialState << G4endl; |
---|
[1196] | 117 | |
---|
[962] | 118 | // Variables to describe the excited configuration |
---|
[819] | 119 | G4double exEnergy = theInitialState.GetExcitationEnergy(); |
---|
[1315] | 120 | G4int A = theInitialState.GetA_asInt(); |
---|
| 121 | G4int Z = theInitialState.GetZ_asInt(); |
---|
| 122 | |
---|
| 123 | G4NistManager* nist = G4NistManager::Instance(); |
---|
[1196] | 124 | |
---|
| 125 | // JMQ 150909: first step in de-excitation chain (SMM will be used only here) |
---|
| 126 | |
---|
[1315] | 127 | // In case A <= 1 the fragment will not perform any nucleon emission |
---|
| 128 | if (A <= 1) |
---|
[819] | 129 | { |
---|
[1315] | 130 | theResults.push_back( theInitialStatePtr ); |
---|
[962] | 131 | } |
---|
[1315] | 132 | // check if a fragment is stable |
---|
| 133 | else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) |
---|
[962] | 134 | { |
---|
[1315] | 135 | theResults.push_back( theInitialStatePtr ); |
---|
| 136 | } |
---|
| 137 | else // In all cases apply once theFermiModel, theMultiFragmentation or theEvaporation |
---|
| 138 | { |
---|
[1196] | 139 | // JMQ 150909: first step in de-excitation is treated separately |
---|
| 140 | // Fragments after the first step are stored in theEvapList |
---|
[1315] | 141 | // Statistical Multifragmentation will take place only once |
---|
[1196] | 142 | // |
---|
[1315] | 143 | // Fermi Break-Up always first |
---|
| 144 | if((A < 5) || (A<GetMaxA() && Z<GetMaxZ())) |
---|
[1196] | 145 | { |
---|
| 146 | theTempResult = theFermiModel->BreakItUp(theInitialState); |
---|
[1315] | 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 | } |
---|
[1196] | 155 | } |
---|
| 156 | else if (exEnergy>GetMinE()*A) |
---|
| 157 | { |
---|
| 158 | theTempResult = theMultiFragmentation->BreakItUp(theInitialState); |
---|
[1315] | 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 | } |
---|
[1196] | 167 | } |
---|
| 168 | else |
---|
| 169 | { |
---|
| 170 | theTempResult = theEvaporation->BreakItUp(theInitialState); |
---|
| 171 | } |
---|
| 172 | |
---|
[1228] | 173 | G4bool deletePrimary = true; |
---|
[1315] | 174 | G4int nsec=theTempResult->size(); |
---|
| 175 | if(nsec > 0) |
---|
[1228] | 176 | { |
---|
[1315] | 177 | // Sort out secondary fragments |
---|
[1228] | 178 | G4FragmentVector::iterator j; |
---|
| 179 | for (j = theTempResult->begin(); j != theTempResult->end(); ++j) |
---|
| 180 | { |
---|
| 181 | if((*j) == theInitialStatePtr) { deletePrimary = false; } |
---|
[1315] | 182 | A = (*j)->GetA_asInt(); |
---|
[1228] | 183 | |
---|
[1315] | 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 | } |
---|
[1228] | 198 | } |
---|
[1196] | 199 | } |
---|
[1228] | 200 | if( deletePrimary ) { delete theInitialStatePtr; } |
---|
[1196] | 201 | delete theTempResult; |
---|
[1228] | 202 | } |
---|
| 203 | // |
---|
| 204 | // JMQ 150909: Further steps in de-excitation chain follow .. |
---|
[1196] | 205 | |
---|
[1228] | 206 | //G4cout << "## After first step " << theEvapList.size() << " for evap; " |
---|
| 207 | // << theEvapStableList.size() << " for photo-evap; " |
---|
| 208 | // << theResults.size() << " results. " << G4endl; |
---|
| 209 | |
---|
| 210 | // ------------------------------ |
---|
| 211 | // De-excitation loop |
---|
| 212 | // ------------------------------ |
---|
[1196] | 213 | |
---|
[1228] | 214 | for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList) |
---|
| 215 | { |
---|
[1315] | 216 | A = (*iList)->GetA_asInt(); |
---|
| 217 | Z = (*iList)->GetZ_asInt(); |
---|
[1228] | 218 | |
---|
[1315] | 219 | // Fermi Break-Up |
---|
| 220 | if ((A < 5) || (A < GetMaxA() && Z < GetMaxZ())) |
---|
[962] | 221 | { |
---|
[1315] | 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 | } |
---|
[1228] | 229 | } |
---|
[1315] | 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) |
---|
[1196] | 243 | { |
---|
[1315] | 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 |
---|
[962] | 250 | { |
---|
[1315] | 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 |
---|
[1196] | 258 | } |
---|
[1315] | 259 | } else { theEvapList.push_back(*j); } // hot fragment |
---|
[962] | 260 | } |
---|
[1228] | 261 | } |
---|
[1315] | 262 | } |
---|
| 263 | if( deletePrimary ) { delete (*iList); } |
---|
| 264 | delete theTempResult; |
---|
[1228] | 265 | } // end of the loop over theEvapList |
---|
| 266 | |
---|
| 267 | //G4cout << "## After 2nd step " << theEvapList.size() << " was evap; " |
---|
| 268 | // << theEvapStableList.size() << " for photo-evap; " |
---|
| 269 | // << theResults.size() << " results. " << G4endl; |
---|
[1196] | 270 | |
---|
[962] | 271 | // ----------------------- |
---|
| 272 | // Photon-Evaporation loop |
---|
| 273 | // ----------------------- |
---|
[1196] | 274 | |
---|
[1315] | 275 | // normally should not reach this point - it is kind of work around |
---|
[1228] | 276 | for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList) |
---|
[819] | 277 | { |
---|
[1315] | 278 | // 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) |
---|
[962] | 284 | { |
---|
[1315] | 285 | G4FragmentVector::iterator j; |
---|
| 286 | for (j = theTempResult->begin(); j != theTempResult->end(); ++j) |
---|
[962] | 287 | { |
---|
[1315] | 288 | theResults.push_back(*j); |
---|
| 289 | } |
---|
| 290 | } |
---|
| 291 | delete theTempResult; |
---|
[1228] | 292 | |
---|
[1315] | 293 | // priamry fragment is kept |
---|
| 294 | theResults.push_back(*iList); |
---|
| 295 | |
---|
[962] | 296 | } // end of photon-evaporation loop |
---|
[1228] | 297 | |
---|
| 298 | //G4cout << "## After 3d step " << theEvapList.size() << " was evap; " |
---|
| 299 | // << theEvapStableList.size() << " was photo-evap; " |
---|
| 300 | // << theResults.size() << " results. " << G4endl; |
---|
| 301 | |
---|
| 302 | #ifdef debug |
---|
| 303 | CheckConservation(theInitialState,*theResults); |
---|
| 304 | #endif |
---|
| 305 | |
---|
| 306 | G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector; |
---|
| 307 | |
---|
| 308 | // MAC (24/07/08) |
---|
| 309 | // To optimise the storing speed, we reserve space in memory for the vector |
---|
| 310 | theReactionProductVector->reserve( theResults.size() ); |
---|
| 311 | |
---|
| 312 | G4int theFragmentA, theFragmentZ; |
---|
| 313 | G4LorentzVector theFragmentMomentum; |
---|
| 314 | |
---|
| 315 | std::list<G4Fragment*>::iterator i; |
---|
| 316 | for (i = theResults.begin(); i != theResults.end(); ++i) |
---|
[819] | 317 | { |
---|
[1228] | 318 | theFragmentA = static_cast<G4int>((*i)->GetA()); |
---|
| 319 | theFragmentZ = static_cast<G4int>((*i)->GetZ()); |
---|
| 320 | theFragmentMomentum = (*i)->GetMomentum(); |
---|
| 321 | G4ParticleDefinition* theKindOfFragment = 0; |
---|
[1315] | 322 | if (theFragmentA == 0) { // photon or e- |
---|
| 323 | theKindOfFragment = (*i)->GetParticleDefinition(); |
---|
[1228] | 324 | } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron |
---|
| 325 | theKindOfFragment = G4Neutron::NeutronDefinition(); |
---|
| 326 | } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton |
---|
| 327 | theKindOfFragment = G4Proton::ProtonDefinition(); |
---|
| 328 | } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron |
---|
| 329 | theKindOfFragment = G4Deuteron::DeuteronDefinition(); |
---|
| 330 | } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton |
---|
| 331 | theKindOfFragment = G4Triton::TritonDefinition(); |
---|
| 332 | } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3 |
---|
| 333 | theKindOfFragment = G4He3::He3Definition(); |
---|
| 334 | } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha |
---|
| 335 | theKindOfFragment = G4Alpha::AlphaDefinition();; |
---|
| 336 | } else { |
---|
| 337 | theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ); |
---|
| 338 | } |
---|
| 339 | if (theKindOfFragment != 0) |
---|
| 340 | { |
---|
| 341 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
---|
| 342 | theNew->SetMomentum(theFragmentMomentum.vect()); |
---|
| 343 | theNew->SetTotalEnergy(theFragmentMomentum.e()); |
---|
| 344 | theNew->SetFormationTime((*i)->GetCreationTime()); |
---|
| 345 | theReactionProductVector->push_back(theNew); |
---|
| 346 | } |
---|
| 347 | delete (*i); |
---|
[819] | 348 | } |
---|
[1228] | 349 | |
---|
| 350 | return theReactionProductVector; |
---|
[819] | 351 | } |
---|
| 352 | |
---|
[962] | 353 | |
---|
| 354 | |
---|
[819] | 355 | G4ReactionProductVector * |
---|
| 356 | G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const |
---|
| 357 | { |
---|
| 358 | if (theFragmentVector == 0) return 0; |
---|
| 359 | |
---|
| 360 | // Conversion from G4FragmentVector to G4ReactionProductVector |
---|
| 361 | G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition(); |
---|
| 362 | G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition(); |
---|
| 363 | G4ParticleDefinition *theProton = G4Proton::ProtonDefinition(); |
---|
| 364 | G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition(); |
---|
| 365 | G4ParticleDefinition *theTriton = G4Triton::TritonDefinition(); |
---|
| 366 | G4ParticleDefinition *theHelium3 = G4He3::He3Definition(); |
---|
| 367 | G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition(); |
---|
| 368 | G4ParticleDefinition *theKindOfFragment = 0; |
---|
| 369 | theNeutron->SetVerboseLevel(2); |
---|
| 370 | G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector; |
---|
[962] | 371 | |
---|
| 372 | // MAC (24/07/08) |
---|
| 373 | // To optimise the storing speed, we reserve space in memory for the vector |
---|
| 374 | theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) ); |
---|
| 375 | |
---|
[819] | 376 | G4int theFragmentA, theFragmentZ; |
---|
| 377 | G4LorentzVector theFragmentMomentum; |
---|
| 378 | |
---|
| 379 | G4FragmentVector::iterator i; |
---|
| 380 | for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) { |
---|
| 381 | // std::cout << (*i) <<'\n'; |
---|
| 382 | theFragmentA = static_cast<G4int>((*i)->GetA()); |
---|
| 383 | theFragmentZ = static_cast<G4int>((*i)->GetZ()); |
---|
| 384 | theFragmentMomentum = (*i)->GetMomentum(); |
---|
| 385 | theKindOfFragment = 0; |
---|
| 386 | if (theFragmentA == 0 && theFragmentZ == 0) { // photon |
---|
| 387 | theKindOfFragment = theGamma; |
---|
| 388 | } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron |
---|
| 389 | theKindOfFragment = theNeutron; |
---|
| 390 | } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton |
---|
| 391 | theKindOfFragment = theProton; |
---|
| 392 | } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron |
---|
| 393 | theKindOfFragment = theDeuteron; |
---|
| 394 | } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton |
---|
| 395 | theKindOfFragment = theTriton; |
---|
| 396 | } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3 |
---|
| 397 | theKindOfFragment = theHelium3; |
---|
| 398 | } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha |
---|
| 399 | theKindOfFragment = theAlpha; |
---|
| 400 | } else { |
---|
| 401 | theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ); |
---|
| 402 | } |
---|
| 403 | if (theKindOfFragment != 0) |
---|
| 404 | { |
---|
| 405 | G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment); |
---|
| 406 | theNew->SetMomentum(theFragmentMomentum.vect()); |
---|
| 407 | theNew->SetTotalEnergy(theFragmentMomentum.e()); |
---|
| 408 | theNew->SetFormationTime((*i)->GetCreationTime()); |
---|
| 409 | #ifdef PRECOMPOUND_TEST |
---|
| 410 | theNew->SetCreatorModel((*i)->GetCreatorModel()); |
---|
| 411 | #endif |
---|
| 412 | theReactionProductVector->push_back(theNew); |
---|
| 413 | } |
---|
| 414 | } |
---|
| 415 | if (theFragmentVector != 0) |
---|
| 416 | { |
---|
| 417 | std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment()); |
---|
| 418 | delete theFragmentVector; |
---|
| 419 | } |
---|
| 420 | G4ReactionProductVector::iterator debugit; |
---|
| 421 | for(debugit=theReactionProductVector->begin(); |
---|
| 422 | debugit!=theReactionProductVector->end(); debugit++) |
---|
| 423 | { |
---|
| 424 | if((*debugit)->GetTotalEnergy()<1.*eV) |
---|
| 425 | { |
---|
| 426 | if(getenv("G4DebugPhotonevaporationData")) |
---|
| 427 | { |
---|
| 428 | G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl; |
---|
| 429 | G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = " |
---|
| 430 | << (*debugit)->GetTotalEnergy()/MeV << "MeV" |
---|
| 431 | << G4endl; |
---|
| 432 | } |
---|
| 433 | delete (*debugit); |
---|
| 434 | *debugit = 0; |
---|
| 435 | } |
---|
| 436 | } |
---|
| 437 | G4ReactionProduct* tmpPtr=0; |
---|
| 438 | theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(), |
---|
| 439 | theReactionProductVector->end(), |
---|
| 440 | std::bind2nd(std::equal_to<G4ReactionProduct*>(), |
---|
| 441 | tmpPtr)), |
---|
| 442 | theReactionProductVector->end()); |
---|
| 443 | return theReactionProductVector; |
---|
| 444 | } |
---|
| 445 | |
---|
| 446 | |
---|
| 447 | #ifdef debug |
---|
| 448 | void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState, |
---|
| 449 | G4FragmentVector * Result) const |
---|
| 450 | { |
---|
| 451 | G4double ProductsEnergy =0; |
---|
| 452 | G4ThreeVector ProductsMomentum; |
---|
| 453 | G4int ProductsA = 0; |
---|
| 454 | G4int ProductsZ = 0; |
---|
| 455 | G4FragmentVector::iterator h; |
---|
| 456 | for (h = Result->begin(); h != Result->end(); h++) { |
---|
| 457 | G4LorentzVector tmp = (*h)->GetMomentum(); |
---|
| 458 | ProductsEnergy += tmp.e(); |
---|
| 459 | ProductsMomentum += tmp.vect(); |
---|
| 460 | ProductsA += static_cast<G4int>((*h)->GetA()); |
---|
| 461 | ProductsZ += static_cast<G4int>((*h)->GetZ()); |
---|
| 462 | } |
---|
| 463 | |
---|
| 464 | if (ProductsA != theInitialState.GetA()) { |
---|
| 465 | G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 466 | G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments" |
---|
| 467 | << G4endl; |
---|
| 468 | G4cout << "Initial A = " << theInitialState.GetA() |
---|
| 469 | << " Fragments A = " << ProductsA << " Diference --> " |
---|
| 470 | << theInitialState.GetA() - ProductsA << G4endl; |
---|
| 471 | } |
---|
| 472 | if (ProductsZ != theInitialState.GetZ()) { |
---|
| 473 | G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 474 | G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments" |
---|
| 475 | << G4endl; |
---|
| 476 | G4cout << "Initial Z = " << theInitialState.GetZ() |
---|
| 477 | << " Fragments Z = " << ProductsZ << " Diference --> " |
---|
| 478 | << theInitialState.GetZ() - ProductsZ << G4endl; |
---|
| 479 | } |
---|
| 480 | if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) { |
---|
| 481 | G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 482 | G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments" |
---|
| 483 | << G4endl; |
---|
| 484 | G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV" |
---|
| 485 | << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> " |
---|
| 486 | << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl; |
---|
| 487 | } |
---|
| 488 | if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || |
---|
| 489 | std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV || |
---|
| 490 | std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) { |
---|
| 491 | G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 492 | G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments" |
---|
| 493 | << G4endl; |
---|
| 494 | G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV" |
---|
| 495 | << " Fragments P = " << ProductsMomentum << " MeV Diference --> " |
---|
| 496 | << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl; |
---|
| 497 | } |
---|
| 498 | return; |
---|
| 499 | } |
---|
| 500 | #endif |
---|
| 501 | |
---|
| 502 | |
---|
| 503 | |
---|
| 504 | |
---|