Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/models/radioactive_decay
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/radioactive_decay/History

    r819 r962  
    1414     * Reverse chronological order (last date on top), please *
    1515     ----------------------------------------------------------
     16
     1709 July 2008 Dennis Wright (radioactive_decay-V09-01-02)
     18- replace exit() with G4Exception in G4RadioactiveDecay and G4NuclearDecayChannel
     19
     2017 June 2008 Fan Lei (radioactive_decay-V09-01-01)
     21- GRIsotopeTable.cc
     22        i) change default verbosity level to 1
     23        ii) correct use G4cout instead of G4cerr 
     24
     2501 May 2008 Fan Lei (radioactive_decay-V09-01-00)
     26- G4NuclearDecayChannel.cc
     27        i) ARM is no longer applied in photon-evaporation for IT mode and
     28        is now applied at the end in DecayIt()
     29        ii) use the correct shell index in appling ARM and switch on Auger
     30        electron production
     31        iii) check the residual kinetic energy after ARM and add it to the atom
     32
    163321 June 2007 Fan Lei (radioactive_decay-V08-03-00)
    1734- Minor changes to remove compilation warnings on Windows

  • trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc

    r819 r962  
    6666#include "G4BetaFermiFunction.hh"
    6767#include "G4PhotonEvaporation.hh"
     68#include "G4AtomicTransitionManager.hh"
     69#include "G4AtomicShell.hh"
    6870#include "G4AtomicDeexcitation.hh"
    69 
    7071
    7172const G4double G4NuclearDecayChannel:: pTolerance = 0.001;
    … …  
    241242  if (decayMode == 2) theParentMass -= 2*0.511 * MeV;
    242243  //
     244#ifdef G4VERBOSE
    243245  if (GetVerboseLevel()>1) {
    244     G4cout << "G4NuclearDecayChannel::DecayIt ";
     246    G4cout << "G4NuclearDecayChannel::DecayIt "<< G4endl;
    245247    G4cout << "the decay mass = " << theParentMass << G4endl;
    246248  }
     249#endif
    247250
    248251  SetParentMass (theParentMass);
    … …  
    261264    {
    262265    case 0:
    263       if (GetVerboseLevel()>0)
    264         {
    265           G4cout << "G4NuclearDecayChannel::DecayIt ";
    266           G4cout << " daughters not defined " <<G4endl;
    267         }
     266      G4cerr << "G4NuclearDecayChannel::DecayIt ";
     267      G4cerr << " daughters not defined " <<G4endl;
    268268      break;
    269269    case 1:
    … …  
    281281      G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "G4NuclearDecayChannel::DecayIt");
    282282    }
    283   if ((products == 0) && (GetVerboseLevel()>0)) {
     283  if (products == 0) {
    284284    G4cerr << "G4NuclearDecayChannel::DecayIt ";
    285285    G4cerr << *parent_name << " can not decay " << G4endl;
    286286    DumpInfo();
    287287  }
    288 
    289   //
    290   // now we have to take care of the EC product which have go through the ARM
     288  //
     289  // If the decay is to an excited state of the daughter nuclide, we need
     290  // to apply the photo-evaporation process.
     291  //
     292  // needed to hold the shell idex after ICM
     293  G4int shellIndex = -1;
     294  //
     295  if (daughterExcitation > 0.0)
     296    {
     297      //
     298      // Pop the daughter nucleus off the product vector - we need to retain
     299      // the momentum of this particle.
     300      //
     301      dynamicDaughter = products->PopProducts();
     302      G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum();
     303      G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum));
     304      //
     305      //
     306      // Now define a G4Fragment with the correct A, Z and excitation, and declare and
     307      // initialise a G4PhotonEvaporation object.
     308      //   
     309      G4Fragment nucleus(daughterA, daughterZ, daughterMomentum);
     310      G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation;
     311      deexcitation->SetVerboseLevel(GetVerboseLevel());
     312      // switch on/off internal electron conversion
     313      deexcitation->SetICM(true);
     314      // set the maximum life-time for a level that will be treated. Level with life-time longer than this
     315      // will be outputed as meta-stable isotope
     316      //
     317      deexcitation->SetMaxHalfLife(1e-6*second);
     318      // but in IT mode, we need to force the transition
     319      if (decayMode == 0) {
     320        deexcitation->RDMForced(true);
     321      } else {
     322        deexcitation->RDMForced(false);
     323      }
     324      //
     325      // Get the gammas by deexciting the nucleus.
     326      //
     327      G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus);
     328      // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide
     329      // as its last entry.
     330      G4int nGammas=gammas->size()-1;
     331      //
     332      // Go through each gamma/e- and add it to the decay product.  The angular distribution
     333      // of the gammas is isotropic, and the residual nucleus is assumed not to have suffered
     334      // any recoil as a result of this de-excitation.
     335      //
     336      for (G4int ig=0; ig<nGammas; ig++)
     337        {
     338          G4DynamicParticle *theGammaRay = new
     339            G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(),
     340                               gammas->operator[](ig)->GetMomentum());
     341          theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime());
     342          products->PushProducts (theGammaRay);
     343        }
     344      //
     345      //      now the nucleus
     346      G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
     347      // f.lei (03/01/03) this is needed to fix the crach in test18
     348      if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ;
     349     
     350      // f.lei (07/03/05) added the delete to fix bug#711
     351      if (dynamicDaughter) delete dynamicDaughter;
     352     
     353      G4IonTable *theIonTable =  (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());     
     354      dynamicDaughter = new G4DynamicParticle
     355        (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation),
     356         daughterMomentum1);
     357      products->PushProducts (dynamicDaughter);
     358      // retrive the ICM shell index
     359      shellIndex = deexcitation->GetVacantShellNumber();
     360     
     361      //
     362      // Delete/reset variables associated with the gammas.
     363      //
     364      while (!gammas->empty()) {
     365        delete *(gammas->end()-1);
     366        gammas->pop_back();
     367      }
     368      //    gammas->clearAndDestroy();
     369      delete gammas;
     370      delete deexcitation;
     371    }
     372  //
     373  // now we have to take care of the EC product which have to go through the ARM
     374  //
     375  G4int eShell = -1;
    291376  if (decayMode == 3 || decayMode == 4 || decayMode == 5) {
    292     G4int eShell = 0;
    293377    switch (decayMode)
    294378      {
    … …  
    296380        //
    297381        {
    298           eShell = 1;
     382          eShell = 0; // --> 0 from 1 (f.lei 30/4/2008)
    299383        }
    300384        break;
    … …  
    313397      case ERROR:
    314398      default:
    315         G4cout << " There is an  error in decay mode selection! exit RDM now" << G4endl;
    316         exit(0);                     
     399        G4Exception("G4NuclearDecayChannel::DecayIt()", "601",
     400                    FatalException, "Error in decay mode selection");
    317401      }
     402  }
     403  // now deal with the IT case where ICM may have been applied
     404  //
     405  if (decayMode == 0) {
     406    eShell = shellIndex;
     407  }
     408  // now apply ARM if there is a vaccancy
     409  //
     410  if (eShell != -1) {
    318411    G4int aZ = daughterZ;
    319     if (aZ > 5 && aZ < 101) {  // only applies to 5< Z <101
    320       G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation();
    321       //no Auger electron generation
    322       // atomDeex->ActivateAugerElectronProduction(0);
    323       std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,eShell);
    324 
    325       // pop up the daughter before insertion
     412    if (aZ > 5 && aZ < 100) {  // only applies to 5< Z <100
     413      // Retrieve the corresponding identifier and binding energy of the selected shell
     414      const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
     415      const G4AtomicShell* shell = transitionManager->Shell(aZ, eShell);
     416      G4double bindingEnergy = shell->BindingEnergy();
     417      G4int shellId = shell->ShellId();
     418
     419      G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation(); 
     420      //the default is no Auger electron generation.
     421      // Switch it on/off here!
     422      atomDeex->ActivateAugerElectronProduction(true);
     423      std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,shellId);
     424
     425      // pop up the daughter before insertion;
     426      // f.lei (30/04/2008) check if the total kinetic energy is less than
     427      // the shell binding energy; if true add the difference to the daughter to conserve the energy
    326428      dynamicDaughter = products->PopProducts();
    327       for (size_t i = 0;  i < armProducts->size(); i++)
     429      G4double tARMEnergy = 0.0;
     430      for (size_t i = 0;  i < armProducts->size(); i++) {
    328431        products->PushProducts ((*armProducts)[i]);
     432        tARMEnergy += (*armProducts)[i]->GetKineticEnergy();
     433      }
     434      if ((bindingEnergy - tARMEnergy) > 0.1*keV){
     435        G4double dEnergy = dynamicDaughter->GetKineticEnergy() + (bindingEnergy - tARMEnergy);
     436        dynamicDaughter->SetKineticEnergy(dEnergy);
     437      }
     438      products->PushProducts(dynamicDaughter);
     439
     440#ifdef G4VERBOSE
     441      if (GetVerboseLevel()>0)
     442        {
     443          G4cout <<"G4NuclearDecayChannel::Selected shell number for ARM =  " <<shellId <<G4endl;
     444          G4cout <<"G4NuclearDecayChannel::ARM products =  " <<armProducts->size()<<G4endl;
     445          G4cout <<"                 The binding energy =  " << bindingEnergy << G4endl;             
     446          G4cout <<"  Total ARM particle kinetic energy =  " << tARMEnergy << G4endl;             
     447        }
     448#endif   
     449
    329450      delete armProducts;
    330451      delete atomDeex;
    331       products->PushProducts (dynamicDaughter);
    332     }
    333   }
    334 
    335   //
    336   // If the decay is to an excited state of the daughter nuclide, we need
    337   // to apply the photo-evaporation process.
    338   //
    339   if (daughterExcitation > 0.0)
    340     {
    341       //
    342       //
    343       // Pop the daughter nucleus off the product vector - we need to retain
    344       // the momentum of this particle.
    345       //
    346       dynamicDaughter = products->PopProducts();
    347       G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum();
    348       G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum));
    349       //
    350       //
    351       // Now define a G4Fragment with the correct A, Z and excitation, and declare and
    352       // initialise a G4DiscreteGammaDeexcitation object.
    353       //
    354    
    355       //  daughterMomentum.setT(daughterMomentum.t()+G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( daughterZ, daughterA )+daughterExcitation);
    356    
    357       //    daughterMomentum.setT(daughterMomentum.t()+daughterExcitation);
    358       G4Fragment nucleus(daughterA, daughterZ, daughterMomentum);
    359       //G4LorentzVector p4(0.,0.,0.,G4NucleiProperties::GetNuclearMass(daughterA,daughterZ)
    360       //               +daughterExcitation);
    361       //G4Fragment nucleus(daughterA, daughterZ, p4);
    362       //    nucleus.SetExcitationEnergy(daughterExcitation);
    363  
    364       // G4VGammaDeexcitation* deexcitation = new G4DiscreteGammaDeexcitation;
    365       G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation;
    366       deexcitation->SetVerboseLevel(GetVerboseLevel());
    367       //    deexcitation->Initialize(nucleus);
    368       deexcitation->SetICM(true);
    369       if (decayMode == 0) {
    370         deexcitation->RDMForced(true);
    371       } else {
    372         deexcitation->RDMForced(false);
    373       }
    374       // ARM in G4 is applied but no auger electrons!
    375       deexcitation->SetARM(true);
    376       // not applied
    377       //deexcitation->SetARM(false);
    378       //
    379       deexcitation->SetMaxHalfLife(1e-6*second);
    380       //
    381       // Get the gammas by deexciting the nucleus.
    382       //
    383       G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus);
    384       // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide
    385       // as its last entry.
    386       G4int nGammas=gammas->size()-1;
    387       //
    388       //
    389       // Go through each gamma/e- and add it to the decay product.  The angular distribution
    390       // of the gammas is isotropic, and the residual nucleus is assumed not to suffer
    391       // any recoil as a result of this de-excitation.
    392       //
    393       for (G4int ig=0; ig<nGammas; ig++)
    394         {
    395           //      G4double costheta = 2.0*G4UniformRand() - 1.0;
    396           //      G4double sintheta = std::sqrt((1.0 - costheta) * (1.0+costheta));
    397           //      G4double phi      = twopi * G4UniformRand();
    398           // G4ParticleMomentum gDirection
    399           //  (sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
    400           //G4double gEnergy = gammas->operator[](ig)->GetMomentum().e()
    401           //  - gammas->operator[](ig)->GetParticleDefinition()->GetPDGMass() ;
    402           G4DynamicParticle *theGammaRay = new
    403             G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(),
    404                                gammas->operator[](ig)->GetMomentum());
    405           theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime());
    406           products->PushProducts (theGammaRay);
    407         }
    408       //
    409       //      now the nucleus
    410       G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
    411       // f.lei (03/01/03) this is needed to fix the crach in test18
    412       if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ;
    413       G4IonTable *theIonTable =  (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
    414       // f.lei (07/03/05) added the delete to fix bug#711
    415       if (dynamicDaughter) delete dynamicDaughter;
    416       dynamicDaughter = new G4DynamicParticle
    417         (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation),
    418          daughterMomentum1);
    419       products->PushProducts (dynamicDaughter);
    420       //
    421       // Delete/reset variables associated with the gammas.
    422       //
    423       //    if (nGammas != 0) gammas->clearAndDestroy();
    424       while (!gammas->empty()) {
    425         delete *(gammas->end()-1);
    426         gammas->pop_back();
    427       }
    428       //    gammas->clearAndDestroy();
    429       delete gammas;
    430       delete deexcitation;
    431     }
     452    }
     453  }
     454
    432455  return products;
    433456}

  • trunk/source/processes/hadronic/models/radioactive_decay/src/G4RIsotopeTable.cc

    r819 r962  
    115115    if(fIsotopeNameList[i] == fname) j = i;}
    116116  if (j >=0) {
    117     if (GetVerboseLevel()>0) {
     117    if (GetVerboseLevel()>1) {
    118118    G4cout <<"G4RIsotopeTable::GetIsotope No. : ";
    119119    G4cout <<j<<G4endl;   
    … …  
    144144    fname = GetIsotopeName(Z, A, E);
    145145    fIsotopeNameList.push_back(fname);
    146     if (GetVerboseLevel()>0) {
     146    if (GetVerboseLevel()>1) {
    147147      G4cout <<"G4RIsotopeTable::GetIsotope create: ";
    148148      G4cout <<fname <<G4endl; 
    … …  
    160160  os <<"A"<< A << "Z" << Z <<'[' << std::setprecision(1) << E/keV << ']';
    161161  G4String name = os.str();
    162   if (GetVerboseLevel()>0) {
    163     G4cerr <<"G4RIsotopeTable::GetIsotope Name: ";
    164     G4cerr <<name <<G4endl;   
     162  if (GetVerboseLevel()>1) {
     163    G4cout <<"G4RIsotopeTable::GetIsotope Name: ";
     164    G4cout <<name <<G4endl;   
    165165  }
    166166  return name;
    … …  
    186186  if (!DecaySchemeFile )
    187187  {
    188     if (GetVerboseLevel()>0) {
     188    if (GetVerboseLevel()>1) {
    189189      G4cout <<"G4RIsotopeTable::GetMeanLife() : "
    190190             <<"cannot find ion radioactive decay file: "
    … …  
    230230    if (!found && aE )
    231231      {
    232         if (GetVerboseLevel()>0) {
     232        if (GetVerboseLevel()>1) {
    233233          G4cout <<"G4RIsotopeTable::GetMeanLife() : ";
    234234          G4cout <<"cannot find ion of required excitation E = " << aE << G4endl;
    … …  
    241241    if (!found && !aE )
    242242      {
    243         if (GetVerboseLevel()>0) {
     243        if (GetVerboseLevel()>1) {
    244244          G4cout <<"G4RIsotopeTable::GetMeanLife() : ";
    245245          G4cout <<"cannot find ion of required excitation E = " << aE << G4endl;
    … …  
    251251    DecaySchemeFile.close();
    252252  }
    253   if (GetVerboseLevel()>0) {
     253  if (GetVerboseLevel()>1) {
    254254    G4cout <<"G4RIsotopeTable::GetMeanLifeTime: ";
    255255    G4cout <<lifetime << " for " << GetIsotopeName(Z, A, aE) <<G4endl;   

  • trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc

    r819 r962  
    902902                    case ERROR:
    903903                    default:
    904                       //
    905                       //
    906                       G4cout << " There is an  error in decay mode selection! exit RDM now" << G4endl;
    907                       exit(0);
     904                      G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
     905                                  FatalException, "Error in decay mode selection");
    908906                     
    909907                    }
Note: See TracChangeset for help on using the changeset viewer.