Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/de_excitation/photon_evaporation
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4ContinuumGammaDeexcitation.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4ContinuumGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    6466
    6567  // Destructor
    66   ~G4ContinuumGammaDeexcitation();
     68  virtual ~G4ContinuumGammaDeexcitation();
    6769
    6870  // Functions
    … …  
    7274  virtual G4VGammaTransition* CreateTransition();
    7375
    74   virtual G4bool CanDoTransition() const;
     76  virtual G4bool CanDoTransition();
    7577
    7678private:

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaDeexcitation.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4DiscreteGammaDeexcitation.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    5153// -------------------------------------------------------------------
    5254//
     55
    5356#ifndef G4DiscreteGammaDeexcitation_hh
    5457#define G4DiscreteGammaDeexcitation_hh
    … …  
    7174
    7275  // Destructor
    73   ~G4DiscreteGammaDeexcitation();
     76  virtual ~G4DiscreteGammaDeexcitation();
    7477
    7578  // Functions
    … …  
    7982  virtual G4VGammaTransition * CreateTransition();
    8083
    81   virtual G4bool CanDoTransition() const;
     84  virtual G4bool CanDoTransition();
    8285
    8386  void SetICM(G4bool hl) { _icm = hl; };

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4DiscreteGammaTransition.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4DiscreteGammaTransition.hh,v 1.4 2010/04/25 18:43:21 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    5254//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
    5355//              Added creation time evaluation for products of evaporation
     56//
    5457//     
     58//        19 April 2010, J. M. Quesada.
     59//              Corrections added for taking into account mismatch between tabulated
     60//              gamma energies and level energy differences (fake photons eliminated)
     61//   
    5562// -------------------------------------------------------------------
    5663
    … …  
    6269#include "G4NuclearLevel.hh"
    6370
     71//JMQ 180410
     72class G4NuclearLevelManager;
     73
    6474class G4DiscreteGammaTransition : public G4VGammaTransition
    6575{
    … …  
    6878  // Constructor
    6979  G4DiscreteGammaTransition(const G4NuclearLevel& level);
    70   G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z);
     80//JMQ 180410
     81//  G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z);
     82  G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z, G4int A);
    7183
    7284  // Destructor
    … …  
    99111  G4double _excitation;
    100112  G4double _gammaCreationTime;
     113  //JMQ 180410
     114  G4int _A;
     115  G4int _Z;
     116  G4NuclearLevelManager * _levelManager;
     117  //JMQ 190410
     118  G4double _tolerance;
    101119};
    102120

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4E1Probability.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4E1Probability.hh,v 1.5 2010/05/19 10:21:44 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
     28//
     29//---------------------------------------------------------------------
     30//
     31// Geant4 header G4E1Probability
     32//
     33// by V. Lara (May 2003)
     34//
     35// Modifications:
     36// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
     37//            by usage of G4Pow, integer A and introduction of const members
    2638//
    2739//
    … …  
    3345#include "G4VEmissionProbability.hh"
    3446#include "G4Fragment.hh"
    35 #include "G4VLevelDensityParameter.hh"
     47
     48class G4Pow;
    3649
    3750class G4E1Probability : public G4VEmissionProbability
    … …  
    4053public:
    4154
    42   G4E1Probability() {};
     55  G4E1Probability();
    4356
    44   ~G4E1Probability();
     57  virtual ~G4E1Probability();
    4558
    4659  G4double EmissionProbability(const G4Fragment& frag, const G4double excite);
    … …  
    4861
    4962private:
    50 
    51   // G4E1Probability() {};
    5263
    5364  G4E1Probability(const G4E1Probability& right);
    … …  
    5869
    5970  // Integrator (simple Gaussian quadrature)
    60 
    6171  G4double EmissionIntegration(const G4Fragment& frag, const G4double excite,
    6272                               const G4double lowLim, const G4double upLim,
    6373                               const G4int numIters);
    6474
    65   // G4VLevelDensityParameter* _levelDensity; // Don't need this
     75  // members
     76  G4Pow*   fG4pow;
     77  G4double theLevelDensityParameter;
     78  G4double normC;
    6679
    6780};

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4PhotonEvaporation.hh

    r962 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4PhotonEvaporation.hh,v 1.7 2010/05/11 11:22:14 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    3638//      Creation date: 23 October 1998
    3739//
    38 //      Modifications:
     40//Modifications:
    3941//
    40 //        18 October 2002, Fan Lei (flei@space.qinetiq.com)
     42// 18 October 2002, Fan Lei (flei@space.qinetiq.com)
    4143//   
    4244//        Implementation of Internal Convertion process in discrete deexcitation
    … …  
    5456//            G4ElectronOccupancy _eOccupancy;
    5557//            G4int _vShellNumber;
     58//
     59// 11 May 2010, V.Ivanchenko added EmittedFragment and BreakUpFragment
     60//                           methods
    5661//
    5762// -------------------------------------------------------------------
    … …  
    6166
    6267#include "globals.hh"
    63 #include "G4VPhotonEvaporation.hh"
    6468#include "G4VEvaporationChannel.hh"
    6569#include "G4VEmissionProbability.hh"
    … …  
    6771#include "G4ElectronOccupancy.hh"
    6872
    69 //#define debug
    70 
    7173class G4Fragment;
    7274
    73 class G4PhotonEvaporation : public G4VPhotonEvaporation, public G4VEvaporationChannel {
     75class G4PhotonEvaporation : public G4VEvaporationChannel {
    7476
    7577public:
    … …  
    7981    virtual ~G4PhotonEvaporation();
    8082
     83    virtual void Initialize(const G4Fragment & fragment);
     84
     85    virtual G4Fragment* EmittedFragment(G4Fragment* theNucleus);
     86
     87    virtual G4FragmentVector* BreakUpFragment(G4Fragment* theNucleus);
     88
    8189    virtual G4FragmentVector * BreakItUp(const G4Fragment & nucleus);
    82 
    83     virtual void Initialize(const G4Fragment & fragment);
    8490
    8591    virtual G4FragmentVector * BreakUp(const G4Fragment & nucleus);
    … …  
    99105    void SetEOccupancy( G4ElectronOccupancy  eOccupancy) ;
    100106
    101 
    102107    G4ElectronOccupancy GetEOccupancy () { return _eOccupancy;} ;
    103108   
    … …  
    111116    G4VGammaDeexcitation * _discrDeexcitation;
    112117    G4VGammaDeexcitation * _contDeexcitation;
    113   //    G4VGammaDeexcitation * _cdDeexcitation;
    114118
    115119    G4ElectronOccupancy _eOccupancy;
    116120    G4int _vShellNumber;
    117121
    118     G4Fragment _nucleus;
     122    G4Fragment* _nucleus;
    119123    G4double _gammaE;
    120124
    121125    G4PhotonEvaporation(const G4PhotonEvaporation & right);
    122 
    123126    const G4PhotonEvaporation & operator = (const G4PhotonEvaporation & right);
    124127
    125     // MGP - Check == and != multiple inheritance... must be a mess!
    126128    G4bool operator == (const G4PhotonEvaporation & right) const;
    127129    G4bool operator != (const G4PhotonEvaporation & right) const;
    128130
    129 #ifdef debug
    130     void CheckConservation(const G4Fragment & theInitialState, G4FragmentVector * Result) const;
    131 #endif
     131  //#ifdef debug
     132  //  void CheckConservation(const G4Fragment & theInitialState, G4FragmentVector * Result) const;
     133  //#endif
    132134
    133135

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VGammaDeexcitation.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4VGammaDeexcitation.hh,v 1.8 2010/04/28 14:22:40 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    3739//
    3840//      Modifications:
     41//         
     42//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
     43//              Added creation time evaluation for products of evaporation
     44//
     45//        21 Nov 2001, Fan Lei (flei@space.qinetiq.com)
     46//           Modified GenerateGamma() and UpdateUncleus() for implementation
     47//           of Internal Conversion processs
     48//
     49//        8 March 2002, Fan Lei (flei@space.qinetiq.com)
     50//          Added  SetEO () , GetEO(), UpdateElectrons() to allow the assignment
     51//          and modification of electron configuration.
    3952//
    4053//        18 October 2002, F. Lei
    … …  
    4356//          Added SetVaccantSN(). It is need to to re-set _vSN after each
    4457//          IC happened.
    45 //
    46 //        8 March 2002, Fan Lei (flei@space.qinetiq.com)
    47 //          Added  SetEO () , GetEO(), UpdateElectrons() to allow the assignment
    48 //          and modification of electron configuration.
    49 //         
    50 //        21 Nov 2001, Fan Lei (flei@space.qinetiq.com)
    51 //           Modified GenerateGamma() and UpdateUncleus() for implementation
    52 //           of Internal Conversion processs
    5358//
    54 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
    55 //              Added creation time evaluation for products of evaporation
     59//        28 April 2010, V.Ivanchenko cleanup methods
    5660//
    5761// -------------------------------------------------------------------
     62//
    5863
    5964#ifndef G4VGAMMADEEXCITATION_HH
    … …  
    7075public:
    7176
    72     G4VGammaDeexcitation();
     77  G4VGammaDeexcitation();
    7378
    74     virtual ~G4VGammaDeexcitation();
     79  virtual ~G4VGammaDeexcitation();
    7580
    76     virtual G4VGammaTransition * CreateTransition() = 0;
    77     virtual G4bool CanDoTransition() const = 0;
     81  virtual G4VGammaTransition * CreateTransition() = 0;
     82  virtual G4bool CanDoTransition() = 0;
    7883
    79     // Single gamma transition
    80     virtual G4FragmentVector * DoTransition();
     84  // Single gamma transition
     85  G4FragmentVector * DoTransition();
    8186
    82     // Chain of gamma transitions
    83     virtual G4FragmentVector * DoChain();
     87  // Chain of gamma transitions
     88  G4FragmentVector * DoChain();
    8489
    85     virtual G4Fragment * GenerateGamma();
     90  G4Fragment * GenerateGamma();
    8691
    87     virtual const G4Fragment & GetNucleus() const;
     92  inline G4Fragment* GetNucleus();
    8893
    89     virtual void SetNucleus(const G4Fragment & nucleus);
     94  inline void SetNucleus(G4Fragment* nucleus);
    9095
    91     virtual void SetVerboseLevel(G4int verbose);
     96  inline void SetVerboseLevel(G4int verbose);
    9297
    93     void SetEO(G4ElectronOccupancy eo) { _electronO = eo; };
    94     void SetVaccantSN( G4int val ) { _vSN = val;};
     98  inline void Initialize();
     99
     100  void SetEO(G4ElectronOccupancy eo) { _electronO = eo; };
     101  void SetVaccantSN( G4int val ) { _vSN = val;};
    95102 
    96     G4ElectronOccupancy GetEO() { return _electronO; };   
    97     G4int GetVacantSN() {return _vSN;};
     103  G4ElectronOccupancy GetEO() { return _electronO; };   
     104  G4int GetVacantSN() {return _vSN;};
     105
    98106protected:
    99107
    100     void Initialize();
    101     void UpdateNucleus(const G4Fragment * gamma);
    102     void UpdateElectrons ();
    103     void Update();
     108  void Update();
    104109
    105     G4VGammaTransition * _transition; // Owned pointer
    106     G4int _verbose;
     110  G4VGammaTransition* _transition; // Owned pointer
     111  G4int _verbose;
    107112
    108113private:
    109114
    110     G4Fragment _nucleus;
    111     G4ElectronOccupancy _electronO;
    112     G4int _vSN;
     115  G4Fragment* _nucleus;
     116  G4ElectronOccupancy _electronO;
     117  G4int _vSN;
    113118
    114     G4VGammaDeexcitation(const G4VGammaDeexcitation & right);
    115 
    116     const G4VGammaDeexcitation & operator = (const G4VGammaDeexcitation & right);
    117     G4bool operator == (const G4VGammaDeexcitation & right) const;
    118     G4bool operator != (const G4VGammaDeexcitation & right) const;
     119  G4VGammaDeexcitation(const G4VGammaDeexcitation & right);
     120  const G4VGammaDeexcitation & operator = (const G4VGammaDeexcitation & right);
     121  G4bool operator == (const G4VGammaDeexcitation & right) const;
     122  G4bool operator != (const G4VGammaDeexcitation & right) const;
    119123
    120124};
     125
     126inline G4Fragment* G4VGammaDeexcitation::GetNucleus()
     127{
     128  return _nucleus;
     129}
     130
     131inline void G4VGammaDeexcitation::SetNucleus(G4Fragment* nucleus)
     132{
     133  _nucleus = nucleus;
     134}
     135
     136inline void G4VGammaDeexcitation::SetVerboseLevel(G4int verbose)
     137{
     138  _verbose = verbose;
     139}
     140
     141inline void G4VGammaDeexcitation::Initialize()
     142{
     143  if (_transition != 0) { delete _transition; }
     144  _transition = CreateTransition();
     145  if (_transition != 0) {
     146    _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
     147  }
     148}
    121149
    122150#endif

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/include/G4VPhotonEvaporation.hh

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4VPhotonEvaporation.hh,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaDeexcitation.cc

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4ContinuumGammaDeexcitation.cc,v 1.7 2010/04/30 16:08:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    4042//
    4143//      02 May 2003,   Vladimir Ivanchenko change interface to G4NuclearlevelManager
     44//
     45//      19 April 2010 J. M. Quesada: smaller value of tolerance parameter
    4246//     
    4347// -------------------------------------------------------------------
    … …  
    6165//
    6266
    63 G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation(): _nucleusZ(0), _nucleusA(0), _levelManager(0)
     67G4ContinuumGammaDeexcitation::G4ContinuumGammaDeexcitation()
     68  : _nucleusZ(0), _nucleusA(0), _levelManager(0)
    6469{ }
    6570
    … …  
    7176G4VGammaTransition* G4ContinuumGammaDeexcitation::CreateTransition()
    7277{
    73   G4Fragment nucleus = GetNucleus();
    74   G4int Z = static_cast<G4int>(nucleus.GetZ());
    75   G4int A = static_cast<G4int>(nucleus.GetA());
    76   G4double excitation = nucleus.GetExcitationEnergy();
     78  G4Fragment* nucleus = GetNucleus();
     79  G4int Z = static_cast<G4int>(nucleus->GetZ());
     80  G4int A = static_cast<G4int>(nucleus->GetA());
     81  G4double excitation = nucleus->GetExcitationEnergy();
    7782
    7883  if (_nucleusA != A || _nucleusZ != Z)
    … …  
    8388    }
    8489
    85   if (_verbose > 1)
     90  if (_verbose > 1) {
    8691    G4cout << "G4ContinuumGammaDeexcitation::CreateTransition - Created" << G4endl;
    87 
     92  }
    8893  G4VGammaTransition* gt =  new G4ContinuumGammaTransition(_levelManager,Z,A,excitation,_verbose );
    8994
    … …  
    9297   
    9398
    94 G4bool G4ContinuumGammaDeexcitation::CanDoTransition() const
     99G4bool G4ContinuumGammaDeexcitation::CanDoTransition()
    95100{
    96  G4bool canDo = true;
     101  //JMQ: far too small, creating sometimes continuum gammas instead of the right discrete ones
     102  // (when excitation energy is slightly over maximum discrete  energy): changed
     103  //  G4double tolerance = 10*eV;
     104  const G4double tolerance = CLHEP::keV;
    97105
    98106  if (_transition == 0)
    99107    {
    100       canDo = false;
    101 
    102       if (_verbose > 0)
    103         G4cout
    104           << "G4ContinuumGammaDeexcitation::CanDoTransition - Null transition "
    105           << G4endl;
     108      if (_verbose > 0) {
     109        G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition - Null transition "
     110               << G4endl;
     111      }
     112      return false;
    106113    }
    107114
    108   G4Fragment nucleus = GetNucleus();
    109   G4double excitation = nucleus.GetExcitationEnergy();
     115  G4Fragment* nucleus = GetNucleus();
     116  G4double excitation = nucleus->GetExcitationEnergy();
    110117
    111   G4double A = nucleus.GetA();
    112   G4double Z = nucleus.GetZ();
    113   if (A <2 || Z<3)
     118  //  G4int A = (G4int)nucleus->GetA();
     119  // G4int Z = (G4int)nucleus->GetZ();
     120  if (_nucleusA<2 || _nucleusZ<3)
    114121    {
    115       canDo = false;
    116       if (_verbose > 0)
    117         G4cout
    118           << "G4ContinuumGammaDeexcitation::CanDoTransition - n/p/H"
    119           << G4endl;
     122      if (_verbose > 1) {
     123        G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition - n/p/H"
     124               << G4endl;
     125      }
     126      return false;
    120127    }
    121128
    122 
    123 
    124   if (excitation <= 0.)
     129  if (excitation <= tolerance)
    125130    {
    126       canDo = false;
    127       if (_verbose > 0)
    128         G4cout
    129           << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation <= 0"
    130           << G4endl;
     131      if (_verbose > 1) {
     132        G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation "
     133               << excitation/CLHEP::keV << " keV is too small"
     134               << G4endl;
     135      }
     136      return false;
    131137    }
    132   G4double tolerance = 10*eV;
    133   if (excitation <= (_levelManager->MaxLevelEnergy()+ tolerance))
    134     {
    135       canDo = false;
    136       if (_verbose > 0)
    137       G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation "
    138              << excitation << " below max discrete level "
    139              << _levelManager->MaxLevelEnergy() << G4endl;
     138  if (excitation <= (_levelManager->MaxLevelEnergy() + tolerance))
     139    {
     140      if (_verbose > 0) {
     141        G4cout << "G4ContinuumGammaDeexcitation::CanDoTransition -  Excitation "
     142               << excitation << " below max discrete level "
     143               << _levelManager->MaxLevelEnergy() << G4endl;
     144      }
     145      return false;
    140146    }
    141147 
    142   if (canDo)
    143     { if (_verbose > 1)
    144       G4cout <<"G4ContinuumGammaDeexcitation::CanDoTransition - CanDo"
    145              << G4endl;
    146     }
    147 
    148   return canDo;
    149 
     148  if (_verbose > 1) {
     149    G4cout <<"G4ContinuumGammaDeexcitation::CanDoTransition - CanDo"
     150           << " Eex(keV)= " << excitation/CLHEP::keV
     151           << " Emax(keV)= " << _levelManager->MaxLevelEnergy()/CLHEP::keV
     152           << " Z= " << _nucleusZ << " A= " << _nucleusA
     153           << G4endl;
     154  }
     155  return true;
    150156}
    151157

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaDeexcitation.cc

    r962 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4DiscreteGammaDeexcitation.cc,v 1.15 2010/05/10 07:20:40 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    6567  _rdm(false), _levelManager(0)
    6668{
    67   _tolerance = 0.1 * MeV;
     69  _tolerance = CLHEP::keV;
    6870}
    6971
    70 
    71 G4DiscreteGammaDeexcitation::~G4DiscreteGammaDeexcitation() {}
    72 
     72G4DiscreteGammaDeexcitation::~G4DiscreteGammaDeexcitation()
     73{}
    7374
    7475G4VGammaTransition* G4DiscreteGammaDeexcitation::CreateTransition()
    7576{
    76   G4Fragment nucleus = GetNucleus();
    77   G4int A = static_cast<G4int>(nucleus.GetA());
    78   G4int Z = static_cast<G4int>(nucleus.GetZ());
    79 
     77  G4Fragment* nucleus = GetNucleus();
     78  G4int A = static_cast<G4int>(nucleus->GetA());
     79  G4int Z = static_cast<G4int>(nucleus->GetZ());
     80  //  _verbose =2;
     81  //  G4cout << "G4DiscreteGammaDeexcitation::CreateTransition: " << nucleus << G4endl;
    8082  if (_nucleusA != A || _nucleusZ != Z)
    8183    {
    … …  
    8486      _levelManager = G4NuclearLevelStore::GetInstance()->GetManager(Z,A);
    8587    }
    86 
    87 
    8888
    8989  if (_levelManager->IsValid())
    … …  
    9696        }
    9797       
    98       G4double excitation = nucleus.GetExcitationEnergy();
    99       //      const G4NuclearLevel* level =_levelManager.NearestLevel(excitation, _tolerance);
     98      G4double excitation = nucleus->GetExcitationEnergy();
    10099      const G4NuclearLevel* level =_levelManager->NearestLevel(excitation);
    101100       
    102101      if (level != 0) 
    103102        {
    104           if (_verbose > 0)
     103          if (_verbose > 0) {
    105104            G4cout
    106105              << "G4DiscreteGammaDeexcitation::CreateTransition - Created from level energy "
    107106              << level->Energy() << ", excitation is "
    108107              << excitation << G4endl;
    109           G4DiscreteGammaTransition* dtransit = new G4DiscreteGammaTransition(*level,Z);
     108          }
     109          G4DiscreteGammaTransition* dtransit = new G4DiscreteGammaTransition(*level,Z,A);
    110110          dtransit->SetICM(_icm); 
    111111          return dtransit;
    … …  
    113113      else
    114114        {
    115           if (_verbose > 0)
     115          if (_verbose > 0) {
    116116            G4cout
    117117              << "G4DiscreteGammaDeexcitation::CreateTransition - No transition created from "
    118118              << excitation << " within tolerance " << _tolerance << G4endl;
    119            
     119          }
    120120          return 0;
    121121        }
    122122    }
    123   else return 0;
     123  return 0;
    124124}
    125125
    126126
    127 G4bool G4DiscreteGammaDeexcitation::CanDoTransition() const
     127G4bool G4DiscreteGammaDeexcitation::CanDoTransition()
    128128{
    129129
    … …  
    138138        << G4endl;
    139139  }
    140   G4Fragment nucleus = GetNucleus();
    141140  if (canDo)  {
    142     G4double A = nucleus.GetA();
    143     G4double Z = nucleus.GetZ();
    144     if (Z<2 || A<3 || Z>98)
     141    if (_nucleusZ<2 || _nucleusA<3 || _nucleusZ>98)
    145142      {
    146143        canDo = false;
    … …  
    152149  }
    153150
    154   G4double excitation = nucleus.GetExcitationEnergy();
     151  G4Fragment* nucleus = GetNucleus();
     152  G4double excitation = nucleus->GetExcitationEnergy();
     153  //G4cout << "G4DiscreteGammaDeexcitation::CanDoTransition: " << nucleus << G4endl;
    155154
    156155  if (canDo) {
    157     if (excitation <= 0.) {
     156    if (excitation <= _tolerance) {
    158157      canDo = false;
    159       if (_verbose > 0)
     158      if (_verbose > 0) {
    160159        G4cout
    161160          << "G4DiscreteGammaDeexcitation::CanDoTransition -  Excitation <= 0"
     161          << excitation << "  " << excitation - _tolerance
    162162          << G4endl;
     163      }
    163164    } else {
    164165      if (excitation > _levelManager->MaxLevelEnergy() + _tolerance) canDo = false;
    165       if (excitation < _levelManager->MinLevelEnergy() - _tolerance) canDo = false; 
     166      //if (excitation < _levelManager->MinLevelEnergy() - _tolerance) canDo = false; 
    166167      // The following is a protection to avoid looping in case of elements with very low
    167168      // ensdf levels
    168       if (excitation < _levelManager->MinLevelEnergy() * 0.9) canDo = false; 
     169      //if (excitation < _levelManager->MinLevelEnergy() * 0.9) canDo = false; 
    169170 
    170171      if (_verbose > 0) {

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4DiscreteGammaTransition.cc

    r819 r1315  
    3737//
    3838//      Modifications:
     39//        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
     40//              Added creation time evaluation for products of evaporation
     41//
     42//        21 Nov. 2001, Fan Lei (flei@space.qinetiq.com)
     43//              i) added G4int _nucleusZ initialise it through the constructor
     44//              ii) modified SelectGamma() to allow the generation of conversion electrons   
     45//              iii) added #include G4AtomicShells.hh
     46//     
    3947//        09 Sep. 2002, Fan Lei  (flei@space.qinetiq.com)
    4048//              Added renormalization to determine whether transition leads to
    4149//              electron or gamma in SelectGamma()
    4250//
    43 //        21 Nov. 2001, Fan Lei (flei@space.qinetiq.com)
    44 //              i) added G4int _nucleusZ initialise it through the constructor
    45 //              ii) modified SelectGamma() to allow the generation of conversion electrons     
    46 //              iii) added #include G4AtomicShells.hh
    47 //     
    48 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
    49 //              Added creation time evaluation for products of evaporation
     51//        19 April 2010, J. M. Quesada.
     52//              Corrections added for taking into account mismatch between tabulated
     53//              gamma energies and level energy differences (fake photons eliminated)
     54//
     55//        9 May 2010, V.Ivanchenko
     56//              Removed unphysical corretions of gamma energy; fixed default particle
     57//              as gamma; do not subtract bounding energy in case of electron emmision
    5058//     
    5159// -------------------------------------------------------------------
    … …  
    5563#include "G4RandGeneralTmp.hh"
    5664#include "G4AtomicShells.hh"
     65//JMQ:
     66#include "G4NuclearLevelStore.hh"
     67#include "G4Pow.hh"
    5768
    5869G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level):
    … …  
    6071{ }
    6172
    62 G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z):
     73//JMQ: now A is also needed in the constructor
     74//G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z):
     75G4DiscreteGammaTransition::G4DiscreteGammaTransition(const G4NuclearLevel& level, G4int Z, G4int A):
    6376  _nucleusZ(Z), _orbitE(-1), _bondE(0.), _aGamma(true), _icm(false), _gammaEnergy(0.),
    64   _level(level), _excitation(0.),  _gammaCreationTime(0.)
     77  _level(level), _excitation(0.),  _gammaCreationTime(0.),_A(A),_Z(Z)
    6578{
    6679  _verbose = 0;
     80  //JMQ: added tolerence in the mismatch
     81  _tolerance = CLHEP::keV;
    6782}
    6883
    … …  
    7489void G4DiscreteGammaTransition::SelectGamma()
    7590{
    76 
     91  // default gamma
     92  _aGamma = true;   
    7793  _gammaEnergy = 0.;
    78 
     94 
    7995  G4int nGammas = _level.NumberOfGammas();
    8096  if (nGammas > 0)
    8197    {
    8298      G4double random = G4UniformRand();
    83 
    84       G4int iGamma = 0;
    85       for(iGamma=0;iGamma < nGammas;iGamma++)
     99     
     100      G4int iGamma;
     101      for(iGamma=0; iGamma<nGammas; ++iGamma)
    86102        {
    87103          if(random <= (_level.GammaCumulativeProbabilities())[iGamma])
    88             break;
     104            { break; }
    89105        }
    90 
    91 
     106     
    92107      // Small correction due to the fact that there are mismatches between
    93108      // nominal level energies and emitted gamma energies
    94 
    95       G4double eCorrection = _level.Energy() - _excitation;
    96 
     109     
     110      // 09.05.2010 VI : it is an error ?
     111      G4double eCorrection = _level.Energy() - _excitation;     
    97112      _gammaEnergy = (_level.GammaEnergies())[iGamma] - eCorrection;
    98 
     113           
     114      //JMQ:
     115      //1)If chosen gamma energy is close enough to excitation energy, the later
     116      //  is used instead for gamma dacey to gs (it guarantees energy conservation)
     117      //2)For energy conservation, level energy differences instead of  tabulated
     118      //  gamma energies must be used (origin of final fake photons)
     119     
     120      if(_excitation - _gammaEnergy < _tolerance)     
     121        {
     122          _gammaEnergy =_excitation;
     123        }
     124      else
     125        {               
     126          _levelManager = G4NuclearLevelStore::GetInstance()->GetManager(_Z,_A);
     127          _gammaEnergy = _excitation -
     128            _levelManager->NearestLevel(_excitation - _gammaEnergy)->Energy();
     129        }
     130     
    99131      //  Warning: the following check is needed to avoid loops:
    100132      //  Due essentially to missing nuclear levels in data files, it is
    … …  
    106138      //          but this change needs a more complex revision of actual design.
    107139      //          I leave this for a later revision.
    108 
    109       if (_gammaEnergy < _level.Energy()*10e-5) _gammaEnergy = _excitation;
     140     
     141      if (_gammaEnergy < _level.Energy()*10.e-5) _gammaEnergy = _excitation;
     142
     143      //G4cout << "G4DiscreteGammaTransition::SelectGamma: " << _gammaEnergy
     144      //             << " _icm: " << _icm << G4endl;
     145
    110146      // now decide whether Internal Coversion electron should be emitted instead
    111147      if (_icm) {
    … …  
    144180              }
    145181            }
    146             if (_verbose > 0)
     182            _bondE = G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell);
     183            if (_verbose > 0) {
    147184              G4cout << "G4DiscreteGammaTransition: _nucleusZ = " <<_nucleusZ
    148185                     << " , iShell = " << iShell 
    149                      << " , Shell binding energy = " << G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell) / keV
     186                     << " , Shell binding energy = " << _bondE/keV
    150187                     << " keV " << G4endl;
    151             _bondE = G4AtomicShells::GetBindingEnergy(_nucleusZ, iShell);
    152             _gammaEnergy = _gammaEnergy - _bondE;
     188            }
     189
     190            // 09.05.2010 VI : it is an error - cannot subtract bond energy from
     191            //                 transition energy here
     192            //_gammaEnergy = _gammaEnergy - _bondE;
     193            //G4cout << "_gammaEnergy = " << _gammaEnergy << G4endl;
     194
    153195            _orbitE = iShell;     
    154196            _aGamma = false ;   // emitted is not a gamma now
    155197          }
    156198      }
    157    
    158       G4double tau = _level.HalfLife() / std::log(2.0);
    159 
    160       G4double tMin = 0;
    161       G4double tMax = 10.0 * tau;
     199     
     200      G4double tau = _level.HalfLife() / G4Pow::GetInstance()->logZ(2);
     201
     202      //09.05.2010 VI rewrite samling of decay time
     203      //              assuming ordinary exponential low
     204      _gammaCreationTime = 0.;     
     205      if(tau > 0.0) {  _gammaCreationTime = -tau*log(G4UniformRand()); }
     206
     207      //G4double tMin = 0;
     208      //G4double tMax = 10.0 * tau;
    162209      //  Original code, not very efficent
    163210      //  G4int nBins = 200;
    … …  
    172219      //  G4RandGeneralTmp randGeneral(sampleArray, nBins);
    173220      //G4double random = randGeneral.shoot();
    174  
     221     
    175222      //_gammaCreationTime = tMin + (tMax - tMin) * random;
    176223
    177224      // new code by Fan Lei
    178225      //
     226      /*
    179227      if (tau != 0 )
    180228      {
    … …  
    186234          //       << _gammaCreationTime/second << G4endl;
    187235       } else { _gammaCreationTime=0.; }
     236      */
    188237    }
    189238  return;
    190239}
    191240
    192 
    193 //G4bool G4DiscreteGammaTransition::IsAGamma()
    194 //{
    195 //  return _aGamma;
    196 //}
    197 
    198 
    199241G4double G4DiscreteGammaTransition::GetGammaEnergy()
    200242{
    … …  
    210252{
    211253  _excitation = energy;
    212   return;
    213 }
    214 
    215 
    216 
    217 
    218 
    219 
     254}
     255
     256
     257
     258
     259
     260

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc

    r819 r1315  
    2525//
    2626//
    27 //  Class G4E1Probability.cc
     27// $Id: G4E1Probability.cc,v 1.9 2010/05/25 10:47:24 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
     29//
     30//---------------------------------------------------------------------
     31//
     32// Geant4 class G4E1Probability
     33//
     34// by V. Lara (May 2003)
     35//
     36// Modifications:
     37// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
     38//            by usage of G4Pow, integer A and introduction of const members
     39//
    2840//
    2941
    3042#include "G4E1Probability.hh"
    31 #include "G4ConstantLevelDensityParameter.hh"
     43//#include "G4ConstantLevelDensityParameter.hh"
    3244#include "Randomize.hh"
     45#include "G4Pow.hh"
    3346
    3447// Constructors and operators
    3548//
    3649
    37 G4E1Probability::G4E1Probability(const G4E1Probability& ) : G4VEmissionProbability()
    38 {
    39 
    40   throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::copy_constructor meant to not be accessible");
    41 
    42 }
    43 
    44 const G4E1Probability& G4E1Probability::
    45 operator=(const G4E1Probability& )
    46 {
    47 
    48   throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::operator= meant to not be accessible");
    49   return *this;
    50 }
    51 
    52 G4bool G4E1Probability::operator==(const G4E1Probability& ) const
    53 {
    54 
    55   return false;
    56 
    57 }
    58 
    59 G4bool G4E1Probability::operator!=(const G4E1Probability& ) const
    60 {
    61 
    62   return true;
    63 
    64 }
     50G4E1Probability::G4E1Probability():G4VEmissionProbability()
     51{
     52  G4double x = CLHEP::pi*CLHEP::hbarc;
     53  normC = 1.0 / (x*x);
     54  theLevelDensityParameter = 0.125/MeV;
     55  fG4pow = G4Pow::GetInstance();
     56}
     57
     58G4E1Probability::~G4E1Probability()
     59{}
    6560
    6661// Calculate the emission probability
    … …  
    6863
    6964G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag,
    70                                                  const G4double gammaE)
     65                                              const G4double gammaE)
    7166{
    7267
    … …  
    7671  // the probability density for photon evaporation from U to U - gammaE
    7772  // (U = nucleus excitation energy, gammaE = total evaporated photon
    78   // energy).
    79   // fragment = nuclear fragment BEFORE de-excitation
     73  // energy). Fragment = nuclear fragment BEFORE de-excitation
    8074
    8175  G4double theProb = 0.0;
    8276
    83   const G4double Afrag = frag.GetA();
    84   const G4double Zfrag = frag.GetZ();
    85   const G4double Uexcite = frag.GetExcitationEnergy();
    86 
    87   if( (Uexcite-gammaE) < 0.0 || gammaE < 0 || Uexcite <= 0) return theProb;
     77  G4int Afrag = frag.GetA_asInt();
     78  G4double Uexcite = frag.GetExcitationEnergy();
     79
     80  if( (Uexcite-gammaE) < 0.0 || gammaE < 0) { return theProb; }
    8881
    8982  // Need a level density parameter.
    … …  
    9184  // nuclei).
    9285
    93   static G4ConstantLevelDensityParameter a;
    94 
    95   G4double aLevelDensityParam = a.LevelDensityParameter(static_cast<G4int>(Afrag),
    96                                                         static_cast<G4int>(Zfrag),
    97                                                         Uexcite);
    98 
    99   G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*Uexcite));
    100   G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
    101 
     86  G4double aLevelDensityParam = Afrag*theLevelDensityParameter;
     87
     88  //  G4double levelDensBef = std::exp(2*std::sqrt(aLevelDensityParam*Uexcite));
     89  //  G4double levelDensAft = std::exp(2*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
     90  // VI reduce number of calls to exp
     91  G4double levelDens = 1.0;
     92  if( aLevelDensityParam > 0.0 ) {
     93    levelDens = std::exp(2*(std::sqrt(aLevelDensityParam*(Uexcite-gammaE))
     94      - std::sqrt(aLevelDensityParam*Uexcite)));
     95  }
    10296  // Now form the probability density
    10397
    … …  
    107101  G4double sigma0 = 2.5 * Afrag * millibarn;  // millibarns
    108102
    109   G4double Egdp = (40.3 / std::pow(Afrag,0.2) )*MeV;
     103  G4double Egdp   = (40.3 / fG4pow->powZ(Afrag,0.2) )*MeV;
    110104  G4double GammaR = 0.30 * Egdp;
    111105 
    112   static G4double normC = 1.0 / ((pi * hbarc)*(pi * hbarc));
    113 
    114106  // CD
    115107  //cout<<"  PROB TESTS "<<G4endl;
    … …  
    124116  //cout<<" normC = "<<normC<<G4endl;
    125117
    126   G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
    127   G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
    128            (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
    129 
    130   G4double sigmaAbs = numerator/denominator;
    131 
    132   theProb = normC * sigmaAbs * gammaE*gammaE *
    133             levelDensAft/levelDensBef;
     118  // VI implementation 18.05.2010
     119  G4double gammaE2 = gammaE*gammaE;
     120  G4double gammaR2 = gammaE2*GammaR*GammaR;
     121  G4double egdp2   = gammaE2 - Egdp*Egdp;
     122  G4double sigmaAbs = sigma0*gammaR2/(egdp2*egdp2 + gammaR2);
     123  theProb = normC * sigmaAbs * gammaE2 * levelDens;
     124
     125  // old implementation
     126  //  G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
     127  // G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
     128  //         (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
     129
     130  //G4double sigmaAbs = numerator/denominator;
     131  //theProb = normC * sigmaAbs * gammaE2 * levelDensAft/levelDensBef;
    134132
    135133  // CD
    … …  
    211209}
    212210
    213 G4E1Probability::~G4E1Probability() {}
    214 
    215 
     211

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4PhotonEvaporation.cc

    r962 r1315  
    2424// ********************************************************************
    2525//
    26 
     26// $Id: G4PhotonEvaporation.cc,v 1.13 2010/05/17 11:47:43 flei Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2728//
    2829// -------------------------------------------------------------------
    … …  
    3738//      Creation date: 23 October 1998
    3839//
    39 //      Modifications:
     40// Modifications:
    4041//     
    41 //        8 March 2002, Fan Lei (flei@space.qinetiq.com)
     42// 8 March 2002, Fan Lei (flei@space.qinetiq.com)
    4243//   
    4344//        Implementation of Internal Convertion process in discrete deexcitation
    … …  
    5051//       G4ElectronOccupancy GetEOccupancy () ;
    5152//
     53// 11 May 2010, V.Ivanchenko add implementation of EmittedFragment and
     54//                           BreakUpFragment methods; cleanup logic
     55//
    5256// -------------------------------------------------------------------
     57//
    5358
    5459#include "G4PhotonEvaporation.hh"
    … …  
    8691void G4PhotonEvaporation::Initialize(const G4Fragment& fragment)
    8792{
    88   _nucleus = fragment;
    89   return;
    90 }
    91 
     93  _nucleus = const_cast<G4Fragment*>(&fragment);
     94}
     95
     96G4Fragment* G4PhotonEvaporation::EmittedFragment(G4Fragment* nucleus)
     97{
     98  //G4cout << "G4PhotonEvaporation::EmittedFragment" << G4endl;
     99  _nucleus = nucleus;
     100 
     101  // Do one photon emission by the continues deexcitation 
     102  _contDeexcitation->SetNucleus(_nucleus);
     103  _contDeexcitation->Initialize();
     104
     105  if(_contDeexcitation->CanDoTransition()) { 
     106    G4Fragment* gamma = _contDeexcitation->GenerateGamma();
     107    if(gamma) {
     108      if (_verbose > 0) {
     109        G4cout << "G4PhotonEvaporation::EmittedFragment continium deex: "   
     110               << gamma << G4endl;
     111        G4cout << "   Residual: " << nucleus << G4endl;
     112      }
     113      return gamma;
     114    }
     115  }
     116
     117  // Do one photon emission by the discrete deexcitation
     118  _discrDeexcitation->SetNucleus(_nucleus);
     119  _discrDeexcitation->Initialize();
     120
     121  if(_discrDeexcitation->CanDoTransition()) { 
     122    G4Fragment* gamma = _discrDeexcitation->GenerateGamma();
     123    if(gamma) {
     124      if (_verbose > 0) {
     125        G4cout << "G4PhotonEvaporation::EmittedFragment discrete deex: "   
     126               << gamma << G4endl;
     127        G4cout << "   Residual: " << nucleus << G4endl;
     128      }
     129      return gamma;
     130    }
     131  }
     132
     133  if (_verbose > 0) {
     134    G4cout << "G4PhotonEvaporation unable emit gamma: "
     135           << nucleus << G4endl;
     136  }
     137  return 0;
     138}
     139
     140G4FragmentVector* G4PhotonEvaporation::BreakUpFragment(G4Fragment* nucleus)
     141{
     142  //G4cout << "G4PhotonEvaporation::BreakUpFragment" << G4endl;
     143  // The same pointer of primary nucleus
     144  _nucleus = nucleus;
     145  _contDeexcitation->SetNucleus(_nucleus);
     146  _discrDeexcitation->SetNucleus(_nucleus);
     147
     148  // Do the whole gamma chain
     149  G4FragmentVector* products = _contDeexcitation->DoChain(); 
     150  if( !products ) { products = new G4FragmentVector(); }
     151
     152  if (_verbose > 0) {
     153    G4cout << "G4PhotonEvaporation::BreakUpFragment " << products->size()
     154           << " gammas from ContinuumDeexcitation " << G4endl;
     155    G4cout << "   Residual: " << nucleus << G4endl;
     156  }
     157  // Products from discrete gamma transitions
     158  G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
     159  if(discrProducts) {
     160    _eOccupancy = _discrDeexcitation->GetEO();
     161    _vShellNumber = _discrDeexcitation->GetVacantSN();
     162
     163    // not sure if the following line is needed!
     164    _discrDeexcitation->SetVaccantSN(-1);
     165
     166    if (_verbose > 0) {
     167      G4cout << "G4PhotonEvaporation::BreakUpFragment " << discrProducts->size()
     168             << " gammas from DiscreteDeexcitation " << G4endl;
     169      G4cout << "   Residual: " << nucleus << G4endl;
     170    }
     171    G4FragmentVector::iterator i;
     172    for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
     173      {
     174        products->push_back(*i);
     175      }
     176    delete discrProducts;
     177  }
     178
     179  if (_verbose > 0) {
     180    G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
     181  }
     182  return products;
     183}
    92184
    93185G4FragmentVector* G4PhotonEvaporation::BreakUp(const G4Fragment& nucleus)
    94186{
    95   _nucleus = nucleus;
    96 
    97   G4FragmentVector* products = new G4FragmentVector;
    98  
    99   _contDeexcitation->SetNucleus(nucleus);
    100   _discrDeexcitation->SetNucleus(nucleus);
     187  //G4cout << "G4PhotonEvaporation::BreakUp" << G4endl;
     188  _nucleus = new G4Fragment(nucleus);
     189
     190  _contDeexcitation->SetNucleus(_nucleus);
     191  _discrDeexcitation->SetNucleus(_nucleus);
    101192 
    102193  // Do one photon emission
    … …  
    104195  // Products from continuum gamma transitions
    105196
    106   G4FragmentVector* contProducts = _contDeexcitation->DoTransition(); 
    107 
    108   G4int nCont = 0;
    109   if (contProducts != 0) nCont = contProducts->size();
    110 
    111   G4FragmentVector::iterator i;
    112   if (nCont > 0)
    113     {
    114       G4Fragment modifiedNucleus = _contDeexcitation->GetNucleus();
    115       _discrDeexcitation->SetNucleus(modifiedNucleus);
    116       for (i = contProducts->begin(); i != contProducts->end(); i++)
    117         {
    118           products->push_back(*i);
    119         }
    120     }
    121   else
     197  G4FragmentVector* products = _contDeexcitation->DoTransition(); 
     198  if( !products ) { products = new G4FragmentVector(); }
     199  else if(_verbose > 0) {
     200    G4cout << "G4PhotonEvaporation::BreakUp " << products->size()
     201           << " gammas from ContinuesDeexcitation " << G4endl;
     202    G4cout << "   Residual: " << nucleus << G4endl;
     203  }
     204
     205  if (0 == products->size())
    122206    {
    123207      // Products from discrete gamma transitions
    124208      G4FragmentVector* discrProducts = _discrDeexcitation->DoTransition();
    125      
    126       G4int nDiscr = 0;
    127       if (discrProducts != 0) nDiscr = discrProducts->size();
    128      
    129       if (_verbose > 0)
    130         G4cout << " = BreakUp = " << nDiscr
    131                << " gammas from DiscreteDeexcitation "
    132                << G4endl;
    133      
    134       for (i = discrProducts->begin(); i != discrProducts->end(); i++)
    135         {
    136           products->push_back(*i);
     209
     210      if (discrProducts) {
     211        _eOccupancy = _discrDeexcitation->GetEO();
     212        _vShellNumber = _discrDeexcitation->GetVacantSN();
     213
     214        // not sure if the following line is needed!
     215        _discrDeexcitation->SetVaccantSN(-1);
     216        //
     217        if (_verbose > 0) {
     218          G4cout << " = BreakUp = " << discrProducts->size()
     219                 << " gammas from DiscreteDeexcitation "
     220                 << G4endl;
     221          G4cout << "   Residual: " << nucleus << G4endl;
    137222        }
    138       discrProducts->clear();
    139       delete discrProducts;
     223        G4FragmentVector::iterator i;
     224        for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
     225          {
     226            products->push_back(*i);
     227          }
     228        delete discrProducts;
     229      }
    140230    }
    141 
    142 
    143   _gammaE = 0.;
    144   if (products->size() > 0)
    145     {
    146       _gammaE = (*(products->begin()))->GetMomentum().e();
    147     }
    148 
    149   contProducts->clear();
    150   delete contProducts;  // delete vector, not fragments
    151 
     231 
    152232  // Add deexcited nucleus to products
    153   G4Fragment* finalNucleus = new G4Fragment(_discrDeexcitation->GetNucleus());
    154   products->push_back(finalNucleus);
    155 
    156 
    157   if (_verbose > 0)
     233  products->push_back(_nucleus);
     234
     235  if (_verbose > 0) {
    158236    G4cout << "*-*-*-* Photon evaporation: " << products->size() << G4endl;
     237  }
    159238
    160239  return products;
    161240}
    162241
    163 
    164242G4FragmentVector* G4PhotonEvaporation::BreakItUp(const G4Fragment& nucleus)
    165243{
    166   _nucleus = nucleus;
    167 
    168   G4FragmentVector* products = new G4FragmentVector;
    169 
    170   _contDeexcitation->SetNucleus(nucleus);
    171   _discrDeexcitation->SetNucleus(nucleus);
     244  // The same pointer of primary nucleus
     245  _nucleus = new G4Fragment(nucleus);
     246  _contDeexcitation->SetNucleus(_nucleus);
     247  _discrDeexcitation->SetNucleus(_nucleus);
     248
     249  //G4cout << "G4PhotonEvaporation::BreakItUp:  " << nucleus << G4endl;
    172250
    173251  // Do the whole gamma chain
    174 
    175   G4FragmentVector* contProducts = _contDeexcitation->DoChain(); 
     252  G4FragmentVector* products = _contDeexcitation->DoChain(); 
     253  if( !products ) { products = new G4FragmentVector; }
    176254
    177255  // Products from continuum gamma transitions
    178   G4int nCont = 0;
    179   if (contProducts != 0) nCont = contProducts->size();
    180 
    181   if (_verbose > 0)
    182     G4cout << " = BreakItUp = " << nCont
     256  if (_verbose > 0) {
     257    G4cout << " = BreakItUp = " << products->size()
    183258           << " gammas from ContinuumDeexcitation " << G4endl;
    184 
    185   G4FragmentVector::iterator i;
    186   if (nCont > 0)
    187     {
    188       G4Fragment modifiedNucleus = _contDeexcitation->GetNucleus();
    189       _discrDeexcitation->SetNucleus(modifiedNucleus);
    190       for (i = contProducts->begin(); i != contProducts->end(); i++)
    191         {
    192           products->push_back(*i);
    193         }
    194     }
     259  }
    195260
    196261  // Products from discrete gamma transitions
    197262  G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
    198   _eOccupancy = _discrDeexcitation->GetEO();
    199   _vShellNumber = _discrDeexcitation->GetVacantSN();
    200 
    201   // not sure if the following line is needed!
    202   _discrDeexcitation->SetVaccantSN(-1);
    203 
    204   G4int nDiscr = 0;
    205   if (discrProducts != 0) nDiscr = discrProducts->size();
    206 
    207   if (_verbose > 0)
    208     G4cout << " = BreakItUp = " << nDiscr
    209            << " gammas from DiscreteDeexcitation " << G4endl;
    210 
    211   for (i = discrProducts->begin(); i != discrProducts->end(); i++)
    212     {
    213       products->push_back(*i);
     263  if(discrProducts) {
     264    _eOccupancy = _discrDeexcitation->GetEO();
     265    _vShellNumber = _discrDeexcitation->GetVacantSN();
     266
     267    // not sure if the following line is needed!
     268    _discrDeexcitation->SetVaccantSN(-1);
     269
     270    if (_verbose > 0) {
     271      G4cout << " = BreakItUp = " << discrProducts->size()
     272             << " gammas from DiscreteDeexcitation " << G4endl;
    214273    }
    215 
     274    G4FragmentVector::iterator i;
     275    for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
     276      {
     277        products->push_back(*i);
     278      }
     279    delete discrProducts;
     280  }
    216281  // Add deexcited nucleus to products
    217   G4Fragment* finalNucleus = new G4Fragment(_discrDeexcitation->GetNucleus());
    218   products->push_back(finalNucleus);
    219   if (_verbose > 0)
    220     G4cout << " = BreakItUp = Nucleus added to products" << G4endl;
    221 
    222   if (_verbose > 0)
     282  products->push_back(_nucleus);
     283
     284  if (_verbose > 0) {
    223285    G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
    224 
    225 #ifdef debug
    226   CheckConservation(nucleus,products);
    227 #endif
    228   contProducts->clear();
    229   discrProducts->clear();
    230   delete contProducts;  // delete vector, not fragments
    231   delete discrProducts;
     286  }
    232287  return products;
    233288}
    … …  
    236291{
    237292  G4double prob = 0.;
    238   if (_probAlgorithm != 0) prob = _probAlgorithm->EmissionProbability(_nucleus,_gammaE);
     293  if (_probAlgorithm != 0) {
     294    prob = _probAlgorithm->EmissionProbability(*_nucleus,_gammaE);
     295  }
    239296  return prob;
    240297}

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VGammaDeexcitation.cc

    r819 r1315  
    2424// ********************************************************************
    2525//
     26// $Id: G4VGammaDeexcitation.cc,v 1.17 2010/05/09 17:31:23 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    4547//              Added creation time evaluation for products of evaporation
    4648//
     49//        19 April 2010, J. M. Quesada calculations in CM system
     50//              pending final boost to lab system  (not critical)
     51//
     52//        23 April 2010, V.Ivanchenko rewite kinematic part using PDG formula
     53//                                    for 2-body decay
    4754//
    4855// -------------------------------------------------------------------
    … …  
    6471#include "G4DiscreteGammaTransition.hh"
    6572
     73
    6674G4VGammaDeexcitation::G4VGammaDeexcitation(): _transition(0), _verbose(0),
    6775                                              _electronO (0), _vSN(-1)
    6876{ }
    6977
    70 
    7178G4VGammaDeexcitation::~G4VGammaDeexcitation()
    7279{
    73   //  if (_transition != 0) delete _transition;
    74 }
    75 
     80  if (_transition != 0) { delete _transition; }
     81}
    7682
    7783G4FragmentVector* G4VGammaDeexcitation::DoTransition()
    7884{
    79   // Template method
    80 
    81  Initialize();
    82  G4FragmentVector* products = new G4FragmentVector;
    83  
     85  Initialize();
     86  G4FragmentVector* products = new G4FragmentVector();
     87 
    8488  if (CanDoTransition())
    8589    {
    86      G4Fragment* gamma = GenerateGamma();
    87      if (gamma != 0)
    88        {
    89          products->push_back(gamma);
    90          UpdateNucleus(gamma);
    91          Update();
    92        }
    93     }
    94 
    95   if (_verbose > 1)
     90      G4Fragment* gamma = GenerateGamma();
     91      if (gamma != 0) { products->push_back(gamma); }
     92    }
     93 
     94  if (_verbose > 1) {
    9695    G4cout << "G4VGammaDeexcitation::DoTransition - Transition deleted " << G4endl;
    97 
    98   if (_transition != 0) delete _transition;
    99 
     96  }
     97 
    10098  return products;
    10199}
    … …  
    103101G4FragmentVector* G4VGammaDeexcitation::DoChain()
    104102{
     103  if (_verbose > 1) { G4cout << "G4VGammaDeexcitation::DoChain" << G4endl; }
     104  const G4double tolerance = CLHEP::keV;
     105
    105106  Initialize();
    106   G4FragmentVector* products = new G4FragmentVector;
    107 
     107  G4FragmentVector* products = new G4FragmentVector();
     108 
    108109  while (CanDoTransition())
    109     {
    110       if (_verbose > 5) G4cout << "G4VGammaDeexcitation::DoChain -  Looping" << G4endl;
    111 
     110    {     
     111      _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
    112112      G4Fragment* gamma = GenerateGamma();
    113113      if (gamma != 0)
    114114        {
    115115          products->push_back(gamma);
    116           UpdateNucleus(gamma);
    117           UpdateElectrons ();
     116          //G4cout << "Eex(keV)= " << _nucleus->GetExcitationEnergy()/keV << G4endl;
     117          if(_nucleus->GetExcitationEnergy() <= tolerance) { break; }
     118          Update();
    118119        }
    119      Update();
    120120    }
    121 
    122   if (_verbose > 1)
    123       G4cout << "G4VGammaDeexcitation::DoChain - Transition deleted, end of chain " << G4endl;
    124 
    125   if (_transition != 0) delete _transition;
     121 
     122  if (_verbose > 1) {
     123    G4cout << "G4VGammaDeexcitation::DoChain - Transition deleted, end of chain " << G4endl;
     124  }
    126125 
    127126  return products;
    128127}
    129128
    130 
    131 const G4Fragment& G4VGammaDeexcitation::GetNucleus() const
    132 {
    133   return _nucleus;
    134 }
    135 
    136 
    137 void G4VGammaDeexcitation::SetNucleus(const G4Fragment& nucleus)
    138 {
    139   _nucleus = G4Fragment(nucleus);
    140 }
    141 
    142 
    143129G4Fragment* G4VGammaDeexcitation::GenerateGamma()
    144130{
     131  // 23/04/10 V.Ivanchenko rewrite complitely
    145132  G4double eGamma = 0.;
    146 
     133 
    147134  if (_transition != 0) {
    148135    _transition->SelectGamma();  // it can be conversion electron too
    149136    eGamma = _transition->GetGammaEnergy();
    150   }
     137    if(eGamma <= 0.0) { return 0; }
     138  }
     139  G4double excitation = _nucleus->GetExcitationEnergy() - eGamma;
     140  if(excitation < 0.0) { excitation = 0.0; }
    151141  if (_verbose > 1 && _transition != 0 )
    152142    {
    153       G4cout << "G4VGammaDeexcitation::GenerateGamma - Gamma energy " << eGamma
    154              << " ** New excitation " << _nucleus.GetExcitationEnergy() - eGamma
     143      G4cout << "G4VGammaDeexcitation::GenerateGamma - Edeexc(MeV)= " << eGamma
     144             << " ** left Eexc(MeV)= " << excitation
    155145             << G4endl;
    156146    }
    157147 
    158   // Photon momentum isotropically generated
    159   // the same for electron
    160  
    161   if (eGamma > 0.)
    162     {
    163       G4double cosTheta = 1. - 2. * G4UniformRand();
    164       G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
    165       G4double phi = twopi * G4UniformRand();
    166       G4double pM = eGamma;
    167       G4DiscreteGammaTransition* dtransition = 0;
    168       dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
    169       if ( dtransition && !( dtransition->IsAGamma()) ) {
    170         G4double eMass = G4Electron::ElectronDefinition()->GetPDGMass();
    171         pM =std::sqrt(eGamma*(eGamma + 2.0*eMass));
    172         eGamma = eGamma + eMass;
    173       }
    174       G4ThreeVector pGamma( pM * sinTheta * std::cos(phi),
    175                             pM * sinTheta * std::sin(phi),
    176                             pM * cosTheta );
    177       G4LorentzVector gamma(pGamma, eGamma);
    178       //        gamma.boost(_nucleus.GetMomentum().boostVector() );
    179       G4Fragment* gammaFragment ;
    180       if ( dtransition && !(dtransition->IsAGamma()) ){
    181         gammaFragment = new G4Fragment(gamma,G4Electron::ElectronDefinition());
    182       } else {
    183         gammaFragment = new G4Fragment(gamma,G4Gamma::GammaDefinition());
    184       }
    185       G4double gammaTime = _transition->GetGammaCreationTime();
    186       gammaTime += _nucleus.GetCreationTime();
    187       gammaFragment->SetCreationTime(gammaTime);
    188      
    189       if (_verbose > 1 && dtransition )
    190         G4cout << "G4VGammaDeexcitation::GenerateGamma -  Gamma fragment generated:  "
    191                << (dtransition->IsAGamma() ? " Gamma" : " Electron" ) << G4endl;
    192       return gammaFragment;
    193     }
    194   else
    195     {
    196       return 0;
    197     }
    198 }
    199 
    200 void G4VGammaDeexcitation::UpdateNucleus(const G4Fragment*  gamma)
    201 {
    202   G4LorentzVector p4Gamma = gamma->GetMomentum();
    203   G4ThreeVector pGamma(p4Gamma.vect());
    204  
    205   G4double eGamma = 0.;
    206   if (_transition != 0)
    207     eGamma = _transition->GetGammaEnergy();
     148  // Do complete Lorentz computation
     149
     150  G4LorentzVector lv = _nucleus->GetMomentum();
     151  G4double Mass = _nucleus->GetGroundStateMass() + excitation;
     152
     153  // select secondary
     154  G4ParticleDefinition* gamma = G4Gamma::Gamma();
     155
    208156  G4DiscreteGammaTransition* dtransition = 0;
    209157  dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
    210   if (dtransition && !(dtransition->IsAGamma()) )     
    211     eGamma += dtransition->GetBondEnergy();
    212 
    213   G4LorentzVector p4Nucleus(_nucleus.GetMomentum() );
    214 
    215 // New tetravector calculation:
    216 
    217 //  G4double Mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(_nucleus.GetZ(),_nucleus.GetA());
    218   G4double m1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(static_cast<G4int>(_nucleus.GetZ()),
    219                                                                                static_cast<G4int>(_nucleus.GetA()));
    220   G4double m2 = _nucleus.GetZ() *  G4Proton::Proton()->GetPDGMass() +
    221     (_nucleus.GetA()- _nucleus.GetZ())*G4Neutron::Neutron()->GetPDGMass();
    222 
    223   G4double Mass = std::min(m1,m2);
    224 
    225 
    226   G4double newExcitation = p4Nucleus.mag() - Mass - eGamma;
    227   if(newExcitation < 0)
    228     newExcitation = 0;
    229  
    230   G4ThreeVector p3Residual(p4Nucleus.vect() - pGamma);
    231   G4double newEnergy = std::sqrt(p3Residual * p3Residual +
    232                             (Mass + newExcitation) * (Mass + newExcitation));
    233   G4LorentzVector p4Residual(p3Residual, newEnergy);
    234  
    235   //  G4LorentzVector p4Residual(-pGamma, p4Nucleus.e() - eGamma);
    236   //  p4Residual.boost( _nucleus.GetMomentum().boostVector() );
    237  
    238   // Update excited nucleus parameters
    239 
    240   _nucleus.SetMomentum(p4Residual);
    241   _nucleus.SetCreationTime(gamma->GetCreationTime());
    242 
    243 //  if (_transition != 0)
    244 //  {
    245 //    G4double excitation =_transition->GetEnergyTo();
    246 //    if (excitation < 0.) excitation = 0.0;
    247 //    _nucleus.SetExcitationEnergy(excitation);
    248 //  }
    249 
    250   return;
    251 }
    252 
    253 void G4VGammaDeexcitation::UpdateElectrons( )
    254 {
    255   G4DiscreteGammaTransition* dtransition = 0;
    256   dtransition = dynamic_cast <G4DiscreteGammaTransition*> (_transition);
    257   if (dtransition && !(dtransition->IsAGamma()) ) {   
    258    
     158  if ( dtransition && !( dtransition->IsAGamma()) ) {
     159    gamma = G4Electron::Electron();
    259160    _vSN = dtransition->GetOrbitNumber();   
    260161    _electronO.RemoveElectron(_vSN);
    261   }
    262   return;
     162    lv += G4LorentzVector(0.0,0.0,0.0,CLHEP::electron_mass_c2 - dtransition->GetBondEnergy());
     163  }
     164
     165  // check consistency 
     166  G4double eMass = gamma->GetPDGMass();
     167
     168  G4double Ecm       = lv.mag();
     169  G4ThreeVector bst  = lv.boostVector();
     170
     171  G4double GammaEnergy = 0.5*((Ecm - Mass)*(Ecm + Mass) + eMass*eMass)/Ecm;
     172  if(GammaEnergy <= eMass) { return 0; }
     173
     174  G4double cosTheta = 1. - 2. * G4UniformRand();
     175  G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
     176  G4double phi = twopi * G4UniformRand();
     177  G4double mom = sqrt((GammaEnergy - eMass)*(GammaEnergy + eMass));
     178  G4LorentzVector Gamma4P(mom * sinTheta * std::cos(phi),
     179                          mom * sinTheta * std::sin(phi),
     180                          mom * cosTheta,
     181                          GammaEnergy);
     182  Gamma4P.boost(bst); 
     183  G4Fragment * thePhoton = new G4Fragment(Gamma4P,gamma);
     184
     185  G4double gammaTime = _nucleus->GetCreationTime() + _transition->GetGammaCreationTime();
     186  thePhoton->SetCreationTime(gammaTime);
     187
     188  lv -= Gamma4P;
     189  _nucleus->SetMomentum(lv);
     190  _nucleus->SetCreationTime(gammaTime);
     191
     192  //G4cout << "G4VGammaDeexcitation::GenerateGamma left nucleus: " << _nucleus << G4endl;
     193  return thePhoton;
    263194}
    264195
    … …  
    269200      delete _transition;
    270201      _transition = 0;
    271       if (_verbose > 1)
     202      if (_verbose > 1) {
    272203        G4cout << "G4VGammaDeexcitation::Update - Transition deleted " << G4endl;
    273     }
    274 
     204      }
     205    }
     206 
    275207  _transition = CreateTransition();
    276208  if (_transition != 0)
    277209    {
    278       _transition->SetEnergyFrom(_nucleus.GetExcitationEnergy());
    279       //      if ( _vSN != -1) (dynamic_cast <G4DiscreteGammaTransition*> (_transition))->SetICM(false);
     210      _transition->SetEnergyFrom(_nucleus->GetExcitationEnergy());
     211      // if ( _vSN != -1) (dynamic_cast <G4DiscreteGammaTransition*> (_transition))->SetICM(false);
    280212      // the above line is commented out for bug fix #952. It was intruduced for reason that
    281213      // the k-shell electron is most likely one to be kicked out and there is no time for
    … …  
    283215      // reported in #952
    284216    }
    285 
     217 
    286218  return;
    287219}
    288 
    289 
    290 void G4VGammaDeexcitation::Initialize()
    291 {
    292   _transition = CreateTransition();
    293   if (_transition != 0) {
    294     _transition->SetEnergyFrom(_nucleus.GetExcitationEnergy());
    295   }
    296   return;
    297 }
    298 
    299 
    300 void G4VGammaDeexcitation::SetVerboseLevel(G4int verbose)
    301 {
    302   _verbose = verbose;
    303   return;
    304 }
    305 
    306 
    307 
    308 
    309 
    310 

  • trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4VPhotonEvaporation.cc

    r819 r1315  
    3939//     
    4040// -------------------------------------------------------------------
     41//
     42// $Id: G4VPhotonEvaporation.cc,v 1.3 2010/04/25 18:43:21 vnivanch Exp $
     43// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    4144
    4245
Note: See TracChangeset for help on using the changeset viewer.