Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (14 years ago)
Author:
garnier
Message:

geant4 tag 9.4

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

Legend:

Unmodified
Added
Removed

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

    r1340 r1347  
    1414     * Reverse chronological order (last date on top), please *
    1515     ----------------------------------------------------------
     1618 November 2010  F.Lei (radioactive_decay-V09-03-04)
     17- src/G4RadioactiveDecay.cc:1593 get rid of the compalition warning with gcc4.5.1
     18
     1917 November 2010  F.Lei (radioactive_decay-V09-03-03)
     20- Completed the implementation of generating the activity table in VR mode
     21- Set the default h-l threshold to 1 micros in VR mode.
     22- General improvement in VR mode implementation. 
     23- G4RadioactiveDecaymessenger.cc: icmCMD,armCMD & hlThreshold are available at all states.
     24
     2511 November 2010  F.Lei
     26- further updates to G4RadioactivityTable.hh .cc and G4RadioactiveDecay.hh .cc
     27
     2810 November 2010  Dennis Wright (radioactive_decay-V09-03-02)
     29-------------------------------------------------------------
     30- G4RadioactiveDecay.cc : replaced incorrect use of "->" with "." for G4Track
     31  in DecayIt(const G4Track&, )
     32
     33- G4RadioactivityTable.cc : add include file for <map> and replace map
     34  with std::map in method AddIsotope
     35
     3629 Oct 2010 F. Lei
     37------------------
     38- Added G4RadioactivityTable.hh .cc files for tallying the accumulated radioactivitties in VR mode
     39- Added in G4RadioactiveDecay.hh:
     40        std::vector<G4RadioactivityTable*> GetTheRadioactivityTables() {return theRadioactivityTables;}
     41        // this is how the radioactivity tables can be retrieved by the user
     42        std::vector<G4RadioactivityTable*>      theRadioactivityTables;
     43        std::map<G4int,G4int>             decayWindows;
     44- In G4RadioactiveDecay.cc:
     45        starting line 1291: changes need to setup the radioactivity tables
     46        line 1624: include the track weight in the weight calculation
     47        line 1627: add the rate to the radioactivity tables
     48- improved formatting of all the class files
     49       
    1650
    175111 Oct 2010 F. Lei (radioactive_decay-V09-03-01)

  • trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh

    r1315 r1347  
    6161#include "G4RadioactiveDecayRateVector.hh"
    6262#include "G4RIsotopeTable.hh"
     63#include "G4RadioactivityTable.hh"
     64
    6365#include <vector>
     66
    6467
    6568//class G4UserlimitsForRD;
    … …  
    7477class G4RadioactiveDecay : public G4VRestDiscreteProcess
    7578{
    76   // class description
    77  
    78   // Implementation of nuclear radioactive decay process.
    79   // It simulates the decays of radioactive nuclei, which is submitted to RDM in the form
    80   // of G4Ions. Half-liffe and decay schemes are hold in the Radioactivity database
    81   // All decay products are submitted back to the particle tracking process through
    82   // the G4ParticleChangeForRadDecay object.
    83 
    84   // class description - end
    85  
     79        // class description
     80
     81        // Implementation of nuclear radioactive decay process.
     82        // It simulates the decays of radioactive nuclei, which is submitted to RDM in the form
     83        // of G4Ions. Half-liffe and decay schemes are hold in the Radioactivity database
     84        // All decay products are submitted back to the particle tracking process through
     85        // the G4ParticleChangeForRadDecay object.
     86
     87        // class description - end
     88
    8689public: // with description
    8790
    88   G4RadioactiveDecay (const G4String& processName="RadioactiveDecay");
    89   ~G4RadioactiveDecay();
    90   // constructor and destructor
    91   //
    92 
    93   G4bool IsApplicable(const G4ParticleDefinition &);
    94   // Returns true if the specified isotope is
    95   //  1) defined as "nucleus" and
    96   //  2) it is within theNucleusLimit
    97   //
    98   G4bool IsLoaded (const G4ParticleDefinition &);
    99   // Returns true if the decay table of the specified nuclei is ready.
    100   //
    101   void SelectAVolume(const G4String aVolume);
    102   // Select a logical volume in which RDM applies.
    103   //
    104   void DeselectAVolume(const G4String aVolume);
    105   // remove a logical volume from the RDM applied list
    106   //
    107   void SelectAllVolumes();
    108   // Select all logical volumes for the application of RDM.
    109   //
    110   void DeselectAllVolumes();
    111   // Remove all logical volumes from RDM applications.
    112   //
    113   void SetDecayBias (G4String filename);
    114   //   Sets the decay biasing scheme using the data in "filename"
    115   //
    116   void SetHLThreshold (G4double hl) {halflifethreshold = hl;}
    117   // Set the half-life threshold for isomer production
    118   //
    119   void SetICM (G4bool icm) {applyICM = icm;}
    120   // Enable/disable ICM
    121   //
    122   void SetARM (G4bool arm) {applyARM = arm;}
    123   // Enable/disable ARM
    124   //
    125   void SetSourceTimeProfile (G4String filename) ;
    126   //  Sets source exposure function using histograms in "filename"
    127   //
    128   G4bool IsRateTableReady(const G4ParticleDefinition &);
    129   // Returns true if the coefficient and decay time table for all the
    130   // descendants of the specified isotope are ready.
    131   //
    132   // used in VR decay mode only
    133   //
    134   void AddDecayRateTable(const G4ParticleDefinition &);
    135   // Calculates the coefficient and decay time table for all the descendants
    136   // of the specified isotope.  Adds the calculated table to the private data
    137   // member "theDecayRateTableVector".
    138   //
    139   //
    140   // used in VR decay mode only
    141   //
    142   void GetDecayRateTable(const G4ParticleDefinition &);
    143   // Used to retrieve the coefficient and decay time table for all the
    144   // descendants of the specified isotope from "theDecayRateTableVector"
    145   // and place it in "theDecayRateTable".
    146   //
    147   //
    148   // used in VR decay mode only
    149   //
    150   void SetDecayRate(G4int,G4int,G4double, G4int, std::vector<G4double>,
    151                     std::vector<G4double>);
    152   // Sets "theDecayRate" with data supplied in the arguements.
    153   //
    154   //  //
    155   // used in VR decay mode only
    156   //
    157 
    158   G4DecayTable *LoadDecayTable (G4ParticleDefinition & theParentNucleus);
    159   // Load the decay data of isotope theParentNucleus.
    160   //
    161   inline void  SetVerboseLevel(G4int value) {verboseLevel = value;}
    162   // Sets the VerboseLevel which controls duggering display.
    163   //
    164   inline G4int GetVerboseLevel() const {return verboseLevel;}
    165   // Returns the VerboseLevel which controls level of debugging output.
    166   //
    167   inline void SetNucleusLimits(G4NucleusLimits theNucleusLimits1)
    168   {theNucleusLimits = theNucleusLimits1 ;}
    169   //  Sets theNucleusLimits which specifies the range of isotopes
    170   //  the G4RadioactiveDecay applies.
    171   //
    172   inline G4NucleusLimits GetNucleusLimits() const
    173   {return theNucleusLimits;}
    174   //  Returns theNucleusLimits which specifies the range of isotopes
    175   //  the G4RadioactiveDecay applies.
    176   //
    177   inline void SetAnalogueMonteCarlo (G4bool r ) { AnalogueMC  = r; }
    178   //   Controls whether G4RadioactiveDecay runs in analogue mode or
    179   //   variance reduction mode.
    180   inline void SetFBeta (G4bool r ) { FBeta  = r; }
    181   //   Controls whether G4RadioactiveDecay uses fast beta simulation mode
    182   //
    183   inline G4bool IsAnalogueMonteCarlo () {return AnalogueMC;}
    184   //   Returns true if the simulation is an analogue Monte Carlo, and false if
    185   //   any of the biassing schemes have been selected.
    186   //
    187   inline void SetBRBias (G4bool r ) { BRBias  = r;
    188   SetAnalogueMonteCarlo(0);}
    189   //   Sets whether branching ration bias scheme applies.
    190   //
    191   inline void SetSplitNuclei (G4int r) { NSplit = r; SetAnalogueMonteCarlo(0); }
    192   //  Sets the N number for the Nuclei spliting bias scheme
    193   //
    194   inline G4int GetSplitNuclei () {return NSplit;}
    195   //  Returns the N number used for the Nuclei spliting bias scheme
    196   //
     91        G4RadioactiveDecay (const G4String& processName="RadioactiveDecay");
     92        ~G4RadioactiveDecay();
     93        // constructor and destructor
     94        //
     95
     96        G4bool IsApplicable(const G4ParticleDefinition &);
     97        // Returns true if the specified isotope is
     98        //  1) defined as "nucleus" and
     99        //  2) it is within theNucleusLimit
     100        //
     101        G4bool IsLoaded (const G4ParticleDefinition &);
     102        // Returns true if the decay table of the specified nuclei is ready.
     103        //
     104        void SelectAVolume(const G4String aVolume);
     105        // Select a logical volume in which RDM applies.
     106        //
     107        void DeselectAVolume(const G4String aVolume);
     108        // remove a logical volume from the RDM applied list
     109        //
     110        void SelectAllVolumes();
     111        // Select all logical volumes for the application of RDM.
     112        //
     113        void DeselectAllVolumes();
     114        // Remove all logical volumes from RDM applications.
     115        //
     116        void SetDecayBias (G4String filename);
     117        //   Sets the decay biasing scheme using the data in "filename"
     118        //
     119        void SetHLThreshold (G4double hl) {halflifethreshold = hl;}
     120        // Set the half-life threshold for isomer production
     121        //
     122        void SetICM (G4bool icm) {applyICM = icm;}
     123        // Enable/disable ICM
     124        //
     125        void SetARM (G4bool arm) {applyARM = arm;}
     126        // Enable/disable ARM
     127        //
     128        void SetSourceTimeProfile (G4String filename) ;
     129        //  Sets source exposure function using histograms in "filename"
     130        //
     131        G4bool IsRateTableReady(const G4ParticleDefinition &);
     132        // Returns true if the coefficient and decay time table for all the
     133        // descendants of the specified isotope are ready.
     134        //
     135        // used in VR decay mode only
     136        //
     137        void AddDecayRateTable(const G4ParticleDefinition &);
     138        // Calculates the coefficient and decay time table for all the descendants
     139        // of the specified isotope.  Adds the calculated table to the private data
     140        // member "theDecayRateTableVector".
     141        //
     142        //
     143        // used in VR decay mode only
     144        //
     145        void GetDecayRateTable(const G4ParticleDefinition &);
     146        // Used to retrieve the coefficient and decay time table for all the
     147        // descendants of the specified isotope from "theDecayRateTableVector"
     148        // and place it in "theDecayRateTable".
     149        //
     150        //
     151        // used in VR decay mode only
     152        //
     153        void SetDecayRate(G4int,G4int,G4double, G4int, std::vector<G4double>,
     154                std::vector<G4double>);
     155        // Sets "theDecayRate" with data supplied in the arguements.
     156        //
     157        //  //
     158        // used in VR decay mode only
     159        //
     160
     161        std::vector<G4RadioactivityTable*> GetTheRadioactivityTables() {return theRadioactivityTables;}
     162        // return the vector of G4Radioactivity map
     163        // should be used in VR mode only
     164
     165        G4DecayTable *LoadDecayTable (G4ParticleDefinition & theParentNucleus);
     166        // Load the decay data of isotope theParentNucleus.
     167        //
     168        inline void  SetVerboseLevel(G4int value) {verboseLevel = value;}
     169        // Sets the VerboseLevel which controls duggering display.
     170        //
     171        inline G4int GetVerboseLevel() const {return verboseLevel;}
     172        // Returns the VerboseLevel which controls level of debugging output.
     173        //
     174        inline void SetNucleusLimits(G4NucleusLimits theNucleusLimits1)
     175        {theNucleusLimits = theNucleusLimits1 ;}
     176        //  Sets theNucleusLimits which specifies the range of isotopes
     177        //  the G4RadioactiveDecay applies.
     178        //
     179        inline G4NucleusLimits GetNucleusLimits() const
     180        {return theNucleusLimits;}
     181        //  Returns theNucleusLimits which specifies the range of isotopes
     182        //  the G4RadioactiveDecay applies.
     183        //
     184        inline void SetAnalogueMonteCarlo (G4bool r ) {
     185          AnalogueMC  = r;
     186          if (!AnalogueMC) halflifethreshold = 1e-6*s;
     187        }
     188        //   Controls whether G4RadioactiveDecay runs in analogue mode or
     189        //   variance reduction mode.
     190        inline void SetFBeta (G4bool r ) { FBeta  = r; }
     191        //   Controls whether G4RadioactiveDecay uses fast beta simulation mode
     192        //
     193        inline G4bool IsAnalogueMonteCarlo () {return AnalogueMC;}
     194        //   Returns true if the simulation is an analogue Monte Carlo, and false if
     195        //   any of the biassing schemes have been selected.
     196        //
     197        inline void SetBRBias (G4bool r ) { BRBias  = r;
     198        SetAnalogueMonteCarlo(0);}
     199        //   Sets whether branching ration bias scheme applies.
     200        //
     201        inline void SetSplitNuclei (G4int r) { NSplit = r; SetAnalogueMonteCarlo(0); }
     202        //  Sets the N number for the Nuclei spliting bias scheme
     203        //
     204        inline G4int GetSplitNuclei () {return NSplit;}
     205        //  Returns the N number used for the Nuclei spliting bias scheme
     206        //
    197207public:
    198  
    199   void BuildPhysicsTable(const G4ParticleDefinition &);
     208
     209        void BuildPhysicsTable(const G4ParticleDefinition &);
    200210
    201211protected:
    202212
    203   G4VParticleChange* DecayIt( const G4Track& theTrack,
    204                               const G4Step&  theStep);
    205 
    206   G4DecayProducts* DoDecay(G4ParticleDefinition& theParticleDef);
    207 
    208   G4double GetMeanFreePath(const G4Track& theTrack,
    209                            G4double previousStepSize,
    210                            G4ForceCondition* condition);
    211 
    212   G4double GetMeanLifeTime(const G4Track& theTrack,
    213                            G4ForceCondition* condition);
    214 
    215   G4double GetTaoTime(G4double,G4double);
    216 
    217   G4double GetDecayTime();
    218 
    219   G4int GetDecayTimeBin(const G4double aDecayTime);
     213        G4VParticleChange* DecayIt( const G4Track& theTrack,
     214                const G4Step&  theStep);
     215
     216        G4DecayProducts* DoDecay(G4ParticleDefinition& theParticleDef);
     217
     218        G4double GetMeanFreePath(const G4Track& theTrack,
     219                G4double previousStepSize,
     220                G4ForceCondition* condition);
     221
     222        G4double GetMeanLifeTime(const G4Track& theTrack,
     223                G4ForceCondition* condition);
     224
     225        G4double GetTaoTime(G4double,G4double);
     226
     227        G4double GetDecayTime();
     228
     229        G4int GetDecayTimeBin(const G4double aDecayTime);
    220230
    221231private:
    222232
    223   G4RadioactiveDecay(const G4RadioactiveDecay &right);
    224   G4RadioactiveDecay & operator=(const G4RadioactiveDecay &right);
     233        G4RadioactiveDecay(const G4RadioactiveDecay &right);
     234        G4RadioactiveDecay & operator=(const G4RadioactiveDecay &right);
    225235
    226236private:
    227237
    228   G4RadioactiveDecaymessenger  *theRadioactiveDecaymessenger;
    229 
    230   G4PhysicsTable               *aPhysicsTable;
    231 
    232   G4RIsotopeTable              *theIsotopeTable;
    233 
    234   G4NucleusLimits               theNucleusLimits;
    235 
    236   const G4double                HighestBinValue;
    237   const G4double                LowestBinValue;
    238 
    239   const G4int                   TotBin;
    240 
    241   G4bool                        AnalogueMC;
    242   G4bool                        BRBias;
    243   G4bool                        FBeta;
    244   G4int                         NSplit;
    245 
    246   G4double                      halflifethreshold;
    247   G4bool                        applyICM;
    248   G4bool                        applyARM;
    249 
    250   G4int                         NSourceBin;
    251   G4double                      SBin[99];
    252   G4double                      SProfile[99];
    253 
    254   G4int                         NDecayBin;
    255   G4double                      DBin[99];
    256   G4double                      DProfile[99];
    257 
    258   std::vector<G4String>         LoadedNuclei;
    259   std::vector<G4String>         ValidVolumes;
    260 
    261   G4RadioactiveDecayRate        theDecayRate;
    262   G4RadioactiveDecayRates       theDecayRateVector;
    263   G4RadioactiveDecayRateVector  theDecayRateTable;
    264   G4RadioactiveDecayRateTable   theDecayRateTableVector;
    265 
    266   static const G4double         levelTolerance;
    267 
    268   // Remainder of life time at rest
    269   G4double                      fRemainderLifeTime;
    270 
    271   G4int                         verboseLevel;
    272 
    273 
    274   // ParticleChange for decay process
    275   G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
    276 
    277   inline G4double AtRestGetPhysicalInteractionLength
    278     (const G4Track& track, G4ForceCondition* condition)
    279   {fRemainderLifeTime = G4VRestDiscreteProcess::
    280     AtRestGetPhysicalInteractionLength(track, condition );
    281   return fRemainderLifeTime;}
    282 
    283   inline G4VParticleChange* AtRestDoIt
    284     (const G4Track& theTrack, const G4Step& theStep)
    285   {return DecayIt(theTrack, theStep);}
    286 
    287   inline G4VParticleChange* PostStepDoIt
    288     (const G4Track& theTrack, const G4Step& theStep)
    289   {return DecayIt(theTrack, theStep);}
     238        G4RadioactiveDecaymessenger  *theRadioactiveDecaymessenger;
     239
     240        G4PhysicsTable               *aPhysicsTable;
     241
     242        G4RIsotopeTable              *theIsotopeTable;
     243
     244        G4NucleusLimits               theNucleusLimits;
     245
     246        const G4double                HighestBinValue;
     247        const G4double                LowestBinValue;
     248
     249        const G4int                   TotBin;
     250
     251        G4bool                        AnalogueMC;
     252        G4bool                        BRBias;
     253        G4bool                        FBeta;
     254        G4int                         NSplit;
     255
     256        G4double                      halflifethreshold;
     257        G4bool                        applyICM;
     258        G4bool                        applyARM;
     259
     260        G4int                         NSourceBin;
     261        G4double                      SBin[99];
     262        G4double                      SProfile[99];
     263
     264        G4int                         NDecayBin;
     265        G4double                      DBin[99];
     266        G4double                      DProfile[99];
     267
     268        std::vector<G4String>         LoadedNuclei;
     269        std::vector<G4String>         ValidVolumes;
     270
     271        G4RadioactiveDecayRate        theDecayRate;
     272        G4RadioactiveDecayRates       theDecayRateVector;
     273        G4RadioactiveDecayRateVector  theDecayRateTable;
     274        G4RadioactiveDecayRateTable   theDecayRateTableVector;
     275
     276        // for the radioactivity tables
     277        std::vector<G4RadioactivityTable*>      theRadioactivityTables;
     278        G4int                         decayWindows[99];
     279
     280        //
     281        static const G4double             levelTolerance;
     282
     283        // Remainder of life time at rest
     284        G4double                      fRemainderLifeTime;
     285
     286        G4int                         verboseLevel;
     287
     288
     289        // ParticleChange for decay process
     290        G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
     291
     292        inline G4double AtRestGetPhysicalInteractionLength
     293                (const G4Track& track, G4ForceCondition* condition)
     294        {fRemainderLifeTime = G4VRestDiscreteProcess::
     295        AtRestGetPhysicalInteractionLength(track, condition );
     296        return fRemainderLifeTime;}
     297
     298        inline G4VParticleChange* AtRestDoIt
     299                (const G4Track& theTrack, const G4Step& theStep)
     300        {return DecayIt(theTrack, theStep);}
     301
     302        inline G4VParticleChange* PostStepDoIt
     303                (const G4Track& theTrack, const G4Step& theStep)
     304        {return DecayIt(theTrack, theStep);}
    290305
    291306};

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

    r1340 r1347  
    121121// Constructor
    122122//
    123 G4RadioactiveDecay::G4RadioactiveDecay
    124   (const G4String& processName)
    125   :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
    126    LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0)
     123G4RadioactiveDecay::G4RadioactiveDecay (const G4String& processName)
     124 :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
     125  LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0)
    127126{
    128127#ifdef G4VERBOSE
    … …  
    136135  aPhysicsTable                = 0;
    137136  pParticleChange              = &fParticleChangeForRadDecay;
    138  
     137
    139138  //
    140139  // Now register the Isotopetable with G4IonTable.
    … …  
    156155  NSourceBin  = 1;
    157156  SBin[0]     = 0.* s;
    158   SBin[1]     = 1e10 * s;
     157  SBin[1]     = 1.* s;
    159158  SProfile[0] = 1.;
    160   SProfile[1] = 1.;
     159  SProfile[1] = 0.;
    161160  NDecayBin   = 1;
    162   DBin[0]     = 9.9e9 * s ;
    163   DBin[1]     = 1e10 * s;
     161  DBin[0]     = 0. * s ;
     162  DBin[1]     = 1. * s;
    164163  DProfile[0] = 1.;
    165164  DProfile[1] = 0.;
     165  decayWindows[0] = 0;
     166  G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
     167  theRadioactivityTables.push_back(rTable);
    166168  NSplit      = 1;
    167169  AnalogueMC  = true ;
    … …  
    187189      delete aPhysicsTable;
    188190    }
    189   //  delete theIsotopeTable;
    190   delete theRadioactiveDecaymessenger;
     191  //    delete theIsotopeTable;
     192    delete theRadioactiveDecaymessenger;
    191193}
    192194
    … …  
    197199//
    198200G4bool G4RadioactiveDecay::IsApplicable( const G4ParticleDefinition &
    199   aParticle)
     201                                        aParticle)
    200202{
    201203  //
    … …  
    224226//
    225227G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition &
    226   aParticle)
     228                                    aParticle)
    227229{
    228230  //
    … …  
    232234  //
    233235  return std::binary_search(LoadedNuclei.begin(),
    234                        LoadedNuclei.end(),
    235                        aParticle.GetParticleName());
     236                            LoadedNuclei.end(),
     237                            aParticle.GetParticleName());
    236238}
    237239////////////////////////////////////////////////////////////////////////////////
    … …  
    242244void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
    243245{
    244  
     246
    245247  G4LogicalVolumeStore *theLogicalVolumes;
    246248  G4LogicalVolume *volume;
    … …  
    331333    G4cout << " RDM removed from all volumes" << G4endl;
    332334#endif
    333  
     335
    334336}
    335337
    … …  
    340342//
    341343G4bool G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition &
    342   aParticle)
     344                                            aParticle)
    343345{
    344346  //
    … …  
    351353    {
    352354      if (theDecayRateTableVector[i].GetIonName() == aParticleName)
    353         return true ;
     355        return true ;
    354356    }
    355357  return false;
    … …  
    363365//
    364366void G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition &
    365   aParticle)
     367                                           aParticle)
    366368{
    367369
    … …  
    376378    }
    377379#ifdef G4VERBOSE
    378         if (GetVerboseLevel()>0)
    379           {
    380             G4cout <<"The DecayRate Table for "
    381                    << aParticleName << " is selected." <<  G4endl;
    382           }
     380  if (GetVerboseLevel()>0)
     381    {
     382      G4cout <<"The DecayRate Table for "
     383             << aParticleName << " is selected." <<  G4endl;
     384    }
    383385#endif
    384386}
    … …  
    388390// GetTaoTime
    389391//
    390 // to perform the convilution of the sourcetimeprofile function with the
     392// to perform the convolution of the source time profile function with the
    391393// decay constants in the decay chain.
    392394//
    … …  
    408410  }
    409411  taotime +=  SProfile[nbin] * (1-std::exp(-(t-SBin[nbin])/tao));
     412  if (taotime < 0.)  {
     413    G4cout <<" Tao time: " <<taotime << " reset to zero!"<<G4endl;
     414    taotime = 0.;
     415  }
     416
    410417#ifdef G4VERBOSE
    411418  if (GetVerboseLevel()>1)
    412419    {G4cout <<" Tao time: " <<taotime <<G4endl;}
    413420#endif
    414   return  taotime ;
    415 }
    416  
     421  return taotime ;
     422}
     423
    417424////////////////////////////////////////////////////////////////////////////////
    418425//
    … …  
    433440#ifdef G4VERBOSE
    434441  if (GetVerboseLevel()>1)
    435     {G4cout <<" Decay time: " <<decaytime <<"[ns]" <<G4endl;}
     442    {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
    436443#endif
    437444  return  decaytime;       
    … …  
    447454G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
    448455{
    449   for (G4int i = 0; i < NDecayBin; i++)
    450     {
    451       if ( aDecayTime > DBin[i]) return i+1;     
    452     }
    453   return  1;
     456  G4int i = 0;
     457  while ( aDecayTime > DBin[i] ) i++;
     458  return  i;
    454459}
    455460////////////////////////////////////////////////////////////////////////////////
    … …  
    472477    G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
    473478    G4double theLife = theParticleDef->GetPDGLifeTime();
    474      
     479
    475480#ifdef G4VERBOSE
    476481    if (GetVerboseLevel()>2)
    … …  
    485490    else if (theLife < 0.0) {meanlife = DBL_MAX;}
    486491    else {meanlife = theLife;}
    487   // set the meanlife to zero for excited istopes which is not in the RDM database
     492    // set the meanlife to zero for excited istopes which is not in the RDM database
    488493    if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
    489494  }
    490495#ifdef G4VERBOSE
    491496  if (GetVerboseLevel()>1)
    492     {G4cout <<"mean life time: " <<meanlife <<"[ns]" <<G4endl;}
    493 #endif
    494  
     497    {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
     498#endif
     499
    495500  return  meanlife;
    496501}
    … …  
    507512  // constants
    508513  G4bool isOutRange ;
    509  
     514
    510515  // get particle
    511516  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
    512  
     517
    513518  // returns the mean free path in GEANT4 internal units
    514519  G4double pathlength;
    … …  
    516521  G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
    517522  G4double aMass = aParticle->GetMass();
    518  
     523
    519524#ifdef G4VERBOSE
    520525  if (GetVerboseLevel()>2) {
    … …  
    525530  }
    526531#endif
    527  
     532
    528533  // check if the particle is stable?
    529534  if (aParticleDef->GetPDGStable()) {
    530535    pathlength = DBL_MAX;
    531    
     536
    532537  } else if (aCtau < 0.0) {
    533538    pathlength =  DBL_MAX;
    534    
     539
    535540    //check if the particle has very short life time ?
    536541  } else if (aCtau < DBL_MIN) {
    537542    pathlength =  DBL_MIN;
    538    
     543
    539544    //check if zero mass
    540545  } else if (aMass <  DBL_MIN)  {
    … …  
    545550    }
    546551#endif
    547    } else {
    548      //calculate the mean free path
    549      // by using normalized kinetic energy (= Ekin/mass)
    550      G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
    551      if ( rKineticEnergy > HighestBinValue) {
    552        // beta >> 1
    553        pathlength = ( rKineticEnergy + 1.0)* aCtau;
    554      } else if ( rKineticEnergy > LowestBinValue) {
    555        // check if aPhysicsTable exists
    556        if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
    557        // beta is in the range valid for PhysicsTable
    558        pathlength = aCtau *
    559          ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
    560      } else if ( rKineticEnergy < DBL_MIN ) {
    561        // too slow particle
    562 #ifdef G4VERBOSE
    563        if (GetVerboseLevel()>2) {
    564          G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
    565          G4cout << aParticleDef->GetParticleName() << G4endl;
    566          G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
    567        }
    568 #endif
    569        pathlength = DBL_MIN;
    570      } else {
    571        // beta << 1
    572        pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
    573      }
    574    }
    575 #ifdef G4VERBOSE
    576    if (GetVerboseLevel()>1) {
    577      G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
    578    }
    579 #endif
    580    return  pathlength;
     552  } else {
     553    //calculate the mean free path
     554    // by using normalized kinetic energy (= Ekin/mass)
     555    G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
     556    if ( rKineticEnergy > HighestBinValue) {
     557      // beta >> 1
     558      pathlength = ( rKineticEnergy + 1.0)* aCtau;
     559    } else if ( rKineticEnergy > LowestBinValue) {
     560      // check if aPhysicsTable exists
     561      if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
     562      // beta is in the range valid for PhysicsTable
     563      pathlength = aCtau *
     564        ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
     565    } else if ( rKineticEnergy < DBL_MIN ) {
     566      // too slow particle
     567#ifdef G4VERBOSE
     568      if (GetVerboseLevel()>2) {
     569        G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
     570        G4cout << aParticleDef->GetParticleName() << G4endl;
     571        G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
     572      }
     573#endif
     574      pathlength = DBL_MIN;
     575    } else {
     576      // beta << 1
     577      pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
     578    }
     579  }
     580#ifdef G4VERBOSE
     581  if (GetVerboseLevel()>1) {
     582    G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
     583  }
     584#endif
     585  return  pathlength;
    581586}
    582587////////////////////////////////////////////////////////////////////////////////
    … …  
    604609  G4int i;
    605610  for ( i = 0 ; i < TotBin ; i++ ) {
    606       gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
    607       beta  = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
    608       aVector->PutValue(i, beta/gammainv);
     611    gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
     612    beta  = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
     613    aVector->PutValue(i, beta/gammainv);
    609614  }
    610615  aPhysicsTable->insert(aVector);
    … …  
    636641    G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files." << G4endl;
    637642    throw G4HadronicException(__FILE__, __LINE__,
    638               "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
     643                              "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
    639644  }
    640645  G4String dirName = getenv("G4RADIOACTIVEDATA");
    … …  
    647652  G4String file = os.str();
    648653
    649  
     654
    650655  std::ifstream DecaySchemeFile(file);
    651656  G4bool found(false);
    … …  
    664669      modeSumBR[i]       = 0.0;
    665670    }
    666    
     671
    667672    G4bool complete(false);
    668673    char inputChars[80]={' '};
    … …  
    718723            a/=1000.;
    719724            c/=1000.;
    720                  
     725
    721726            //  cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
    722                  
     727
    723728            switch (theDecayMode) {
    724729            case IT:
    … …  
    750755                    e0 = c*MeV/0.511;
    751756                    n = aBetaFermiFunction->GetFFN(e0);
    752                                
     757
    753758                    // now to work out the histogram and initialise the random generator
    754759                    G4int npti = 100;                           
    … …  
    765770                      // Special treatment for K-40 (problem #1068) (flei,06/05/2010)
    766771                      if (theParentNucleus.GetParticleName() == "K40[0.0]") f *=
    767                         (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));
     772                                                                              (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));
    768773                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
    769774                    }             
    … …  
    798803                  if (e0 > 0.) {
    799804                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
    800                                
     805
    801806                    n = aBetaFermiFunction->GetFFN(e0);
    802                                
     807
    803808                    // now to work out the histogram and initialise the random generator
    804809                    G4int npti = 100;                           
    … …  
    824829                    theDecayTable->Insert(aBetaPlusChannel);
    825830                    modeSumBR[2] += b;
    826                                
     831
    827832                    delete[] pdf;
    828833                    delete aBetaFermiFunction;       
    … …  
    909914                G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
    910915                            FatalException, "Error in decay mode selection");
    911                
     916
    912917              }
    913918            }
    … …  
    972977//
    973978void G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE,
    974                                        G4int theG, std::vector<G4double> theRates,
    975                                        std::vector<G4double> theTaos)
     979                                      G4int theG, std::vector<G4double> theRates,
     980                                      std::vector<G4double> theTaos)
    976981{
    977982  //fill the decay rate vector
    … …  
    993998  // 2) Add the coefficiencies to the decay rate table vector
    994999  //
    995  
     1000
    9961001  //
    9971002  // Create and initialise variables used in the method.
    … …  
    9991004
    10001005  theDecayRateVector.clear();
    1001  
     1006
    10021007  G4int nGeneration = 0;
    10031008  std::vector<G4double> rates;
    10041009  std::vector<G4double> taos;
    1005  
     1010
    10061011  // start rate is -1.
     1012  // Eq.4.26 of the Technical Note
    10071013  rates.push_back(-1.);
    10081014  //
    … …  
    10121018  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
    10131019  G4double tao = theParentNucleus.GetPDGLifeTime();
     1020  if (tao < 0.) tao = 1e-30;
    10141021  taos.push_back(tao);
    10151022  G4int nEntry = 0;
    1016  
     1023
    10171024  //fill the decay rate with the intial isotope data
    10181025  SetDecayRate(Z,A,E,nGeneration,rates,taos);
    … …  
    10341041  G4AlphaDecayChannel *theAlphaChannel = 0;
    10351042  G4RadioactiveDecayMode theDecayMode;
    1036   //  G4NuclearLevelManager levelManager;
    1037   //const G4NuclearLevel* level;
    10381043  G4double theBR = 0.0;
    10391044  G4int AP = 0;
    … …  
    10561061  //
    10571062  theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
    1058  
     1063
    10591064  while (!stable) {
    10601065    nGeneration++;
    … …  
    10771082        aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
    10781083      }
     1084       
     1085      G4DecayTable *theDecayTable = new G4DecayTable();
    10791086      aTempDecayTable = aParentNucleus->GetDecayTable();
     1087      for (i=0; i< 7; i++) brs[i] = 0.0;
     1088
    10801089      //
    10811090      //
    10821091      // Go through the decay table and to combine the same decay channels
    10831092      //
    1084       for (i=0; i< 7; i++) brs[i] = 0.0;
    1085      
    1086       G4DecayTable *theDecayTable = new G4DecayTable();
    1087      
    10881093      for (i=0; i<aTempDecayTable->entries(); i++) {
    10891094        theChannel             = aTempDecayTable->GetDecayChannel(i);
    … …  
    10991104
    11001105          if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
    1101            
    1102             // Level hafe life is in ns and the gate as 1 micros by default
    1103             if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){   
    1104               // some further though may needed here
     1106            // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command
     1107            if (level->HalfLife()*ns >= halflifethreshold ){   
     1108              // save the metastable nucleus
    11051109              theDecayTable->Insert(theChannel);
    11061110            }
    … …  
    11201124      brs[3] = brs[4] =brs[5] =  0.0;
    11211125      for (i= 0; i<7; i++){
    1122         if (brs[i] > 0) {
     1126        if (brs[i] > 0.) {
    11231127          switch ( i ) {
    11241128          case 0:
    … …  
    11271131            // Decay mode is isomeric transition.
    11281132            //
    1129            
     1133
    11301134            theITChannel =  new G4ITDecayChannel
    11311135              (0, (const G4Ions*) aParentNucleus, brs[0]);
    1132            
     1136
    11331137            theDecayTable->Insert(theITChannel);
    11341138            break;
    1135            
     1139
    11361140          case 1:
    11371141            //
    … …  
    11421146                                                               brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
    11431147            theDecayTable->Insert(theBetaMinusChannel);
    1144            
     1148
    11451149            break;
    1146            
     1150
    11471151          case 2:
    11481152            //
    … …  
    11541158            theDecayTable->Insert(theBetaPlusChannel);
    11551159            break;                   
    1156            
     1160
    11571161          case 6:
    11581162            //
    … …  
    11641168            theDecayTable->Insert(theAlphaChannel);
    11651169            break;
    1166            
     1170
    11671171          default:
    11681172            break;
    … …  
    11701174        }
    11711175      }
    1172        
    11731176      //
    11741177      // loop over all braches in theDecayTable
    … …  
    11791182        theBR = theChannel->GetBR();
    11801183        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
    1181        
    1182         //
    1183         // now test if the daughterNucleus is a valid one
    1184         //
    1185         if (IsApplicable(*theDaughterNucleus) && theBR
    1186             && aParentNucleus != theDaughterNucleus ) {
     1184        //  first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
     1185        if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1 ) {
     1186            A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
     1187            Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
     1188            theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
     1189        }
     1190        if (IsApplicable(*theDaughterNucleus) && theBR && aParentNucleus != theDaughterNucleus) {
    11871191          // need to make sure daugher has decaytable
    11881192          if (!IsLoaded(*theDaughterNucleus)){
    … …  
    11941198            Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
    11951199            E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
    1196          
     1200
    11971201            TaoPlus = theDaughterNucleus->GetPDGLifeTime();
    11981202            //          cout << TaoPlus <<G4endl;
    1199             if (TaoPlus > 0.) {
    1200               // first set the taos, one simply need to add to the parent ones
    1201               taos.clear();
    1202               taos = TP;
    1203               taos.push_back(TaoPlus);
    1204               // now calculate the coefficiencies
    1205               //
    1206               // they are in two parts, first the les than n ones
    1207               rates.clear();
    1208               size_t k;
    1209               for (k = 0; k < RP.size(); k++){
     1203            if (TaoPlus <= 0.)  TaoPlus = 1e-30;
     1204
     1205            // first set the taos, one simply need to add to the parent ones
     1206            taos.clear();
     1207            taos = TP;
     1208            taos.push_back(TaoPlus);
     1209            // now calculate the coefficiencies
     1210            //
     1211            // they are in two parts, first the less than n ones
     1212            // Eq 4.24 of the TN
     1213            rates.clear();
     1214            size_t k;
     1215            for (k = 0; k < RP.size(); k++){
     1216              if ((TP[k]-TaoPlus) == 0.) {
     1217                theRate =  1e30;
     1218              } else {
    12101219                theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k];
    1211                 rates.push_back(theRate);
    12121220              }
    1213               //
    1214               // the sencond part: the n:n coefficiency
    1215               theRate = 0.;
    1216               for (k = 0; k < RP.size(); k++){
    1217                 theRate -=TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
     1221              rates.push_back(theRate);
     1222            }
     1223            //
     1224            // the sencond part: the n:n coefficiency
     1225            // Eq 4.25 of the TN.  Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013
     1226            //
     1227            theRate = 0.;
     1228            G4double aRate;
     1229            for (k = 0; k < RP.size(); k++){
     1230              if ((TP[k]-TaoPlus) == 0.) {
     1231                aRate =  1e30;
     1232              } else {
     1233                aRate = TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
    12181234              }
    1219               rates.push_back(theRate);               
    1220               SetDecayRate (Z,A,E,nGeneration,rates,taos);
    1221               theDecayRateVector.push_back(theDecayRate);
    1222               nEntry++;
    1223             }
     1235              theRate -= aRate;
     1236            }
     1237            rates.push_back(theRate);         
     1238            SetDecayRate (Z,A,E,nGeneration,rates,taos);
     1239            theDecayRateVector.push_back(theDecayRate);
     1240            nEntry++;
    12241241          }
    12251242        } 
    … …  
    12331250    if (nS == nT) stable = true;
    12341251  }
    1235  
     1252
    12361253  //end of while loop
    12371254  // the calculation completed here
    1238  
    1239  
     1255
     1256
    12401257  // fill the first part of the decay rate table
    12411258  // which is the name of the original particle (isotope)
    … …  
    12461263  // now fill the decay table with the newly completed decay rate vector
    12471264  //
    1248  
     1265
    12491266  theDecayRateTable.SetItsRates(theDecayRateVector);
    1250  
     1267
    12511268  //
    12521269  // finally add the decayratetable to the tablevector
    … …  
    12541271  theDecayRateTableVector.push_back(theDecayRateTable);
    12551272}
    1256  
     1273
    12571274////////////////////////////////////////////////////////////////////////////////
    12581275//
    … …  
    12671284  std::ifstream infile ( filename, std::ios::in );
    12681285  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open source data file" );
    1269  
    1270   float bin, flux;
     1286
     1287  G4double bin, flux;
    12711288  NSourceBin = -1;
    12721289  while (infile >> bin >> flux ) {
    … …  
    12771294  }
    12781295  SetAnalogueMonteCarlo(0);
     1296  infile.close();
     1297
    12791298#ifdef G4VERBOSE
    12801299  if (GetVerboseLevel()>1)
    … …  
    12941313  std::ifstream infile ( filename, std::ios::in);
    12951314  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open bias data file" );
    1296  
    1297   float bin, flux;
     1315
     1316  G4double bin, flux;
     1317  G4int dWindows = 0;
     1318  G4int i ;
     1319
     1320  theRadioactivityTables.clear();
     1321  //  for (i = 0; i<100; i++) decayWindows[i] = -1;
     1322
    12981323  NDecayBin = -1;
    12991324  while (infile >> bin >> flux ) {
    … …  
    13021327    DBin[NDecayBin] = bin * s;
    13031328    DProfile[NDecayBin] = flux;
     1329    if (flux > 0.) {
     1330      decayWindows[NDecayBin] = dWindows;
     1331      dWindows++;
     1332      G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
     1333      theRadioactivityTables.push_back(rTable);
     1334    }
    13041335  }
    1305   G4int i ;
    13061336  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
    13071337  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
     1338  // converted to accumulated probabilities
     1339  //
    13081340  SetAnalogueMonteCarlo(0);
     1341  infile.close();
     1342
    13091343#ifdef G4VERBOSE
    13101344  if (GetVerboseLevel()>1)
    … …  
    13151349////////////////////////////////////////////////////////////////////////////////
    13161350//
    1317 //
    13181351// DecayIt
    13191352//
    1320 G4VParticleChange* G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step& )
    1321 {
    1322   //
     1353G4VParticleChange*
     1354G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
     1355{
    13231356  // Initialize the G4ParticleChange object. Get the particle details and the
    13241357  // decay table.
    1325   //
     1358
    13261359  fParticleChangeForRadDecay.Initialize(theTrack);
    13271360  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
    … …  
    13291362
    13301363  // First check whether RDM applies to the current logical volume
    1331   //
    1332   if(!std::binary_search(ValidVolumes.begin(),
    1333                     ValidVolumes.end(),
    1334                     theTrack.GetVolume()->GetLogicalVolume()->GetName()))
    1335     {
    1336 #ifdef G4VERBOSE
    1337       if (GetVerboseLevel()>0)
    1338         {
    1339           G4cout <<"G4RadioactiveDecay::DecayIt : "
    1340                  << theTrack.GetVolume()->GetLogicalVolume()->GetName()
    1341                  << " is not selected for the RDM"<< G4endl;
    1342           G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
    1343           G4cout << " The Valid volumes are " << G4endl;
    1344           for (size_t i = 0; i< ValidVolumes.size(); i++)
    1345             G4cout << ValidVolumes[i] << G4endl;
    1346         }
    1347 #endif
    1348       fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
    1349       //
    1350       //
    1351       // Kill the parent particle.
    1352       //
    1353       fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
    1354       fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
    1355       ClearNumberOfInteractionLengthLeft();
    1356       return &fParticleChangeForRadDecay;
     1364
     1365  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
     1366                          theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
     1367#ifdef G4VERBOSE
     1368    if (GetVerboseLevel()>0) {
     1369      G4cout <<"G4RadioactiveDecay::DecayIt : "
     1370             << theTrack.GetVolume()->GetLogicalVolume()->GetName()
     1371             << " is not selected for the RDM"<< G4endl;
     1372      G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
     1373      G4cout << " The Valid volumes are " << G4endl;
     1374      for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
    13571375    }
    1358    
     1376#endif
     1377    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
     1378
     1379    // Kill the parent particle.
     1380
     1381    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
     1382    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
     1383    ClearNumberOfInteractionLengthLeft();
     1384    return &fParticleChangeForRadDecay;
     1385  }
     1386
    13591387  // now check is the particle is valid for RDM
    1360   //
    1361   if (!(IsApplicable(*theParticleDef)))
    1362     {
    1363       //
    1364       // The particle is not a Ion or outside the nucleuslimits for decay
    1365       //
    1366 #ifdef G4VERBOSE
    1367       if (GetVerboseLevel()>0)
    1368         {
    1369           G4cerr <<"G4RadioactiveDecay::DecayIt : "
    1370                  <<theParticleDef->GetParticleName()
    1371                  << " is not a valid nucleus for the RDM"<< G4endl;
    1372         }
    1373 #endif
    1374       fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
    1375       //
    1376       //
    1377       // Kill the parent particle.
    1378       //
    1379       fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
    1380       fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
    1381       ClearNumberOfInteractionLengthLeft();
    1382       return &fParticleChangeForRadDecay;
     1388
     1389  if (!(IsApplicable(*theParticleDef))) {
     1390    //
     1391    // The particle is not a Ion or outside the nucleuslimits for decay
     1392    //
     1393#ifdef G4VERBOSE
     1394    if (GetVerboseLevel()>0) {
     1395      G4cerr <<"G4RadioactiveDecay::DecayIt : "
     1396             <<theParticleDef->GetParticleName()
     1397             << " is not a valid nucleus for the RDM"<< G4endl;
    13831398    }
    1384  
     1399#endif
     1400    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
     1401
     1402    //
     1403    // Kill the parent particle.
     1404    //
     1405    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
     1406    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
     1407    ClearNumberOfInteractionLengthLeft();
     1408    return &fParticleChangeForRadDecay;
     1409  }
     1410
    13851411  if (!IsLoaded(*theParticleDef))
    13861412    {
    … …  
    13881414    }
    13891415  G4DecayTable *theDecayTable = theParticleDef->GetDecayTable();
    1390  
     1416
    13911417  if  (theDecayTable == 0 || theDecayTable->entries() == 0 )
    13921418    {
    … …  
    14231449      G4ThreeVector currentPosition;
    14241450      currentPosition = theTrack.GetPosition();
    1425      
     1451
    14261452      // check whether use Analogue or VR implementation
    14271453      //
    … …  
    14541480        G4double ParentEnergy = theParticle->GetTotalEnergy();
    14551481        G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
    1456        
     1482
    14571483        if (theTrack.GetTrackStatus() == fStopButAlive )
    14581484          {
    … …  
    15011527              (products->PopProducts(), finalGlobalTime, currentPosition);
    15021528            secondary->SetGoodForTrackingFlag();
    1503                         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
     1529            secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
    15041530            fParticleChangeForRadDecay.AddSecondary(secondary);
    15051531          }
    … …  
    15211547        if (!IsRateTableReady(*theParticleDef)) {
    15221548          // if the decayrates are not ready, calculate them and
    1523           // add to the rate table vector
     1549          // add to the rate table vector
    15241550          AddDecayRateTable(*theParticleDef);
    15251551        }
    … …  
    15381564        G4double taotime;
    15391565        G4double decayRate;
    1540        
     1566
    15411567        size_t i;
    15421568        size_t j;
    … …  
    15571583        for (G4int n = 0; n < NSplit; n++)
    15581584          {
    1559             /*
    1560             //
    15611585            // Get the decay time following the decay probability function
    15621586            // suppllied by user
    1563             //
     1587           
    15641588            G4double theDecayTime = GetDecayTime();
    1565            
    15661589            G4int nbin = GetDecayTimeBin(theDecayTime);
    15671590           
    1568             // claculate the first part of the weight function
     1591            // calculate the first part of the weight function
    15691592           
    1570             G4double weight1 =1./DProfile[nbin-1]
    1571               *(DBin[nbin]-DBin[nbin-1])
    1572               /NSplit;
    1573             if (nbin > 1) {
    1574                weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
    1575                  *(DBin[nbin]-DBin[nbin-1])
    1576                  /NSplit;}
     1593            G4double weight1 = 1.;
     1594            if (nbin == 1) {
     1595              weight1 = 1./DProfile[nbin-1]
     1596                *(DBin[nbin]-DBin[nbin-1])/NSplit;
     1597            } else if (nbin > 1) {
     1598              weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
     1599                *(DBin[nbin]-DBin[nbin-1])/NSplit;
     1600            }
     1601
    15771602            // it should be calculated in seconds
    15781603            weight1 /= s ;
    1579             */
    1580             //
     1604           
    15811605            // loop over all the possible secondaries of the nucleus
    15821606            // the first one is itself.
    1583             //
    1584             for ( i = 0; i<theDecayRateVector.size(); i++){
     1607
     1608            for (i = 0; i<theDecayRateVector.size(); i++){
    15851609              PZ = theDecayRateVector[i].GetZ();
    15861610              PA = theDecayRateVector[i].GetA();
    … …  
    15891613              PR = theDecayRateVector[i].GetDecayRateC();
    15901614
    1591               //
    1592               // Get the decay time following the decay probability function
    1593               // suppllied by user
    1594               //
    1595               G4double theDecayTime = GetDecayTime();
    1596              
    1597               G4int nbin = GetDecayTimeBin(theDecayTime);
    1598              
    1599               // claculate the first part of the weight function
    1600              
    1601               G4double weight1 =1./DProfile[nbin-1]
    1602                 *(DBin[nbin]-DBin[nbin-1])
    1603                 /NSplit;
    1604               if (nbin > 1) {
    1605                 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
    1606                   *(DBin[nbin]-DBin[nbin-1])
    1607                   /NSplit;}
    1608               // it should be calculated in seconds
    1609               weight1 /= s ;
    1610              
    16111615              // a temprary products buffer and its contents is transfered to
    16121616              // the products at the end of the loop
    1613               //
     1617
    16141618              G4DecayProducts *tempprods = 0;
    1615              
     1619
    16161620              // calculate the decay rate of the isotope
    1617               // one need to fold the the source bias function with the decaytime
    1618               //
     1621              // decayRate is the radioactivity of isotope (PZ,PA,PE) at the
     1622              // time 'theDecayTime'
     1623              // it will be used to calculate the statistical weight of the
     1624              // decay products of this isotope
     1625
     1626              //              G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE  <<G4endl;
    16191627              decayRate = 0.;
    16201628              for ( j = 0; j < PT.size(); j++){
    16211629                taotime = GetTaoTime(theDecayTime,PT[j]);
    16221630                decayRate -= PR[j] * taotime;
     1631                // Eq.4.23 of of the TN
     1632                // note the negative here is required as the rate in the eqation is defined to be negative,
     1633                // i.e. decay away, but we need pasitive value here.
     1634
     1635                //              G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t" << decayRate << G4endl ;             
    16231636              }
    1624              
    1625               // decayRatehe radioactivity of isotope (PZ,PA,PE) at the
    1626               // time 'theDecayTime'
    1627               // it will be used to calculate the statistical weight of the
    1628               // decay products of this isotope
    1629              
    1630              
    1631               //
     1637
    16321638              // now calculate the statistical weight
    1633               //
    1634              
    1635               G4double weight = weight1*decayRate;
     1639
     1640              // one need to fold the the source bias function with the decaytime
     1641              // also need to include the track weight! (F.Lei, 28/10/10)
     1642              G4double weight = weight1*decayRate*theTrack.GetWeight();
     1643               
     1644              // add the isotope to the radioactivity tables
     1645              //                            G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
     1646              //G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
     1647              theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight);
     1648                               
    16361649              // decay the isotope
    16371650              theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
    16381651              parentNucleus = theIonTable->GetIon(PZ,PA,PE);
    1639              
     1652
    16401653              // decide whther to apply branching ratio bias or not
    16411654              //
    … …  
    16441657                ndecaych = G4int(theDecayTable->entries()*G4UniformRand());
    16451658                G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(ndecaych);
    1646                 if (theDecayChannel == 0)
    1647                   {
    1648                     // Decay channel not found.
    1649 #ifdef G4VERBOSE
    1650                     if (GetVerboseLevel()>0)
    1651                       {
    1652                         G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
    1653                         G4cerr <<G4endl;
    1654                         theDecayTable ->DumpInfo();
    1655                       }
    1656 #endif
     1659                if (theDecayChannel == 0) {
     1660                  // Decay channel not found.
     1661#ifdef G4VERBOSE
     1662                  if (GetVerboseLevel()>0) {
     1663                    G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
     1664                    G4cerr <<G4endl;
     1665                    theDecayTable ->DumpInfo();
    16571666                  }
    1658                 else
    1659                   {
    1660                     // A decay channel has been identified, so execute the DecayIt.
    1661                     G4double tempmass = parentNucleus->GetPDGMass();     
    1662                     tempprods = theDecayChannel->DecayIt(tempmass);
    1663                     weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
    1664                   }
    1665               }
    1666               else {
     1667#endif
     1668                } else {
     1669                  // A decay channel has been identified, so execute the DecayIt.
     1670                  G4double tempmass = parentNucleus->GetPDGMass();     
     1671                  tempprods = theDecayChannel->DecayIt(tempmass);
     1672                  weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
     1673                }
     1674              } else {
    16671675                tempprods = DoDecay(*parentNucleus);
    16681676              }
    1669               //
     1677
    16701678              // save the secondaries for buffers
    1671               //
     1679
    16721680              numberOfSecondaries = tempprods->entries();
    16731681              currentTime = finalGlobalTime + theDecayTime;
    1674               for (index=0; index < numberOfSecondaries; index++)
    1675                 {
    1676                   asecondaryparticle = tempprods->PopProducts();
    1677                   if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){
    1678                     pw.push_back(weight);
    1679                     ptime.push_back(currentTime);
    1680                     secondaryparticles.push_back(asecondaryparticle);
    1681                   }
     1682              for (index=0; index < numberOfSecondaries; index++) {
     1683                asecondaryparticle = tempprods->PopProducts();
     1684                if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){
     1685                  pw.push_back(weight);
     1686                  ptime.push_back(currentTime);
     1687                  secondaryparticles.push_back(asecondaryparticle);
    16821688                }
    1683               //
     1689              }
     1690
    16841691              delete tempprods;
    1685              
     1692
    16861693              //end of i loop
    16871694            }
    1688            
     1695
    16891696            // end of n loop
    1690           }
     1697          }
     1698
    16911699        // now deal with the secondaries in the two stl containers
    16921700        // and submmit them back to the tracking manager
    … …  
    16991707                                             secondaryparticles[index], ptime[index], currentPosition);
    17001708            secondary->SetGoodForTrackingFlag();           
    1701                         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
     1709            secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
    17021710            secondary->SetWeight(pw[index]);       
    1703             fParticleChangeForRadDecay.AddSecondary(secondary);
     1711            fParticleChangeForRadDecay.AddSecondary(secondary);
    17041712          }
    17051713        //
    … …  
    17081716        //theTrack.SetTrackStatus(fStopButAlive);
    17091717        //energyDeposit += theParticle->GetKineticEnergy();
    1710        
     1718
    17111719      }
    1712    
     1720
    17131721      //
    17141722      // Kill the parent particle.
    … …  
    17221730      //
    17231731      ClearNumberOfInteractionLengthLeft();
    1724      
     1732
    17251733      return &fParticleChangeForRadDecay ;
    17261734    }
    17271735}
    17281736
    1729 ////////////////////////////////////////////////////////////////////////////////
    1730 //
     1737///////////////////////////////////////////////////////////////////
    17311738//
    17321739// DoDecay
    17331740//
    1734 G4DecayProducts* G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
    1735 {
    1736   G4DecayProducts *products = 0;
    1737   //
    1738   //
     1741G4DecayProducts*
     1742G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
     1743{
     1744  G4DecayProducts* products = 0;
     1745
    17391746  // follow the decaytable and generate the secondaries...
    1740   //
    1741   //
    1742 #ifdef G4VERBOSE
    1743   if (GetVerboseLevel()>0)
    1744     {
    1745       G4cout<<"Begin of DoDecay..."<<G4endl;
     1747
     1748#ifdef G4VERBOSE
     1749  if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
     1750#endif
     1751
     1752  G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();
     1753
     1754  // Choose a decay channel.
     1755
     1756#ifdef G4VERBOSE
     1757  if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
     1758#endif
     1759
     1760  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
     1761  if (theDecayChannel == 0) {
     1762    // Decay channel not found.
     1763
     1764    G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
     1765    G4cerr <<G4endl;
     1766    theDecayTable ->DumpInfo();
     1767  } else {
     1768
     1769    // A decay channel has been identified, so execute the DecayIt.
     1770
     1771#ifdef G4VERBOSE
     1772    if (GetVerboseLevel()>1) {
     1773      G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
     1774      G4cerr <<theDecayChannel <<G4endl;
    17461775    }
    17471776#endif
    1748   G4DecayTable *theDecayTable = theParticleDef.GetDecayTable();
    1749   //
    1750   // Choose a decay channel.
    1751   //
    1752 #ifdef G4VERBOSE
    1753   if (GetVerboseLevel()>0)
    1754     {
    1755       G4cout <<"Selecte a channel..."<<G4endl;
    1756     }
    1757 #endif
    1758   G4VDecayChannel *theDecayChannel = theDecayTable->SelectADecayChannel();
    1759   if (theDecayChannel == 0)
    1760     {
    1761       // Decay channel not found.
    1762       //
    1763       G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
    1764       G4cerr <<G4endl;
    1765       theDecayTable ->DumpInfo();
    1766     }
    1767       else
    1768     {
    1769       //
    1770       // A decay channel has been identified, so execute the DecayIt.
    1771       //
    1772 #ifdef G4VERBOSE
    1773       if (GetVerboseLevel()>1)
    1774         {
    1775           G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
    1776           G4cerr <<theDecayChannel <<G4endl;
    1777         }
    1778 #endif
    1779      
    1780       G4double tempmass = theParticleDef.GetPDGMass();
    1781       //
    1782      
    1783       products = theDecayChannel->DecayIt(tempmass);
    1784      
    1785     }
     1777
     1778    G4double tempmass = theParticleDef.GetPDGMass();
     1779    products = theDecayChannel->DecayIt(tempmass);
     1780  }
     1781
    17861782  return products;
    1787 
    1788 }
    1789 
    1790 
    1791 
    1792 
    1793 
    1794 
    1795 
    1796 
    1797 
     1783}

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

    r1340 r1347  
    123123  icmCmd->SetParameterName("applyICM",true);
    124124  icmCmd->SetDefaultValue(true);
    125   icmCmd->AvailableForStates(G4State_PreInit);
     125  //icmCmd->AvailableForStates(G4State_PreInit);
    126126  //
    127127  // Command contols whether ARM will be applied or not
    … …  
    131131  armCmd->SetParameterName("applyARM",true);
    132132  armCmd->SetDefaultValue(true);
    133   armCmd->AvailableForStates(G4State_PreInit);
     133  //armCmd->AvailableForStates(G4State_PreInit);
    134134  //
    135135  // Command to set the h-l thresold for isomer production
    … …  
    140140  // hlthCmd->SetRange("hlThreshold>0.");
    141141  hlthCmd->SetUnitCategory("Time");
    142   hlthCmd->AvailableForStates(G4State_PreInit);
     142  //  hlthCmd->AvailableForStates(G4State_PreInit);
    143143  //
    144144  // Command to define the incident particle source time profile.
    … …  
    209209                                    SetNucleusLimits(nucleuslimitsCmd->GetNewNucleusLimitsValue(newValues));}
    210210  else if  (command==analoguemcCmd) {
    211     G4int vl;
    212     const char* t = newValues;
    213     std::istringstream is(t);
    214     is >> vl;
    215     theRadioactiveDecayContainer->SetAnalogueMonteCarlo(vl!=0);}
     211    theRadioactiveDecayContainer->SetAnalogueMonteCarlo(analoguemcCmd->GetNewBoolValue(newValues));}
    216212  else if  (command==fbetaCmd) {
    217     G4int vl;
    218     const char* t = newValues;
    219     std::istringstream is(t);
    220     is >> vl;
    221     theRadioactiveDecayContainer->SetFBeta(vl!=0);}
     213    theRadioactiveDecayContainer->SetFBeta(fbetaCmd->GetNewBoolValue(newValues));}
    222214  else if  (command==avolumeCmd) {theRadioactiveDecayContainer->
    223215                                   SelectAVolume(newValues);}
    … …  
    229221                                   DeselectAllVolumes();}
    230222  else if  (command==brbiasCmd) {
    231     G4int vl;
    232     const char* t = newValues;
    233     std::istringstream is(t);
    234     is >> vl;
    235     theRadioactiveDecayContainer->SetBRBias(vl!=0);}
     223    theRadioactiveDecayContainer->SetBRBias(brbiasCmd->GetNewBoolValue(newValues));}
    236224  else if (command==sourcetimeprofileCmd) {theRadioactiveDecayContainer->
    237225                                             SetSourceTimeProfile(newValues);}
Note: See TracChangeset for help on using the changeset viewer.