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/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • 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.