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

Legend:

Unmodified
Added
Removed

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

    r1196 r1315  
    1515     ----------------------------------------------------------
    1616
     1726 May 2010 Dennis Wright (radioactive_decay-V09-03-00)
     18-------------------------------------------------------
     19tag HEAD to submit for testing
     20
     2111 may 2010 F.Lei
     22- G4RadioactiveDecay.cc::LoadDecayTable:line 931 create a decaytable for isomers not included in RDM database and assume
     23  they all under go IT decay.
     24- G4RadioactiveDecay.cc: line 1434 after DoDecay() check if there is any secondary produced. Kill the track if not to
     25  prevent it entering a infinite loop
     26- G4RadioactiveDecay.cc: line 1459. make sure the propertime is positive. negative case occures when the isomer is not in
     27  RDM database and its f-l is set to -1 by default. 
     28- G4NuclearDecayChannel.cc: Added initialisation to all three constructors:
     29  halflifethreshold = 1e-6*second;
     30  applyICM = true;
     31  applyARM = true;
     32
     3310 May 2010 F.Lei
     34- G4RadioactiveDecay.hh:line 264 corrected typo as pointed out by Luciano Pandola
     35- G4RadioactiveDecay.cc: insert a special treatment for K-40 beta decay at line 774 as suggested by Mauro Taiuti.
     36  This fixes bug #1068
     37  applied SetICM() SetARM() and SetHLThreshold() to a decay channel is created.
     38- G4NuclearDecayChannel.cc: line 396. Limit the shell index to < 7, as required by the ARM.
     39  line 329: change to use BreakUp rather than BreakItUp so to limit to one transition per step when ICM is applied (bug #1001)
     40  (Note: changes have also been made in 30/08/2008 to use the correct atomic shell index and to conserve energy atfer ARM,
     41   these problems were pointed out by Andreas Zoglauer in the bug report).
     42- G4RadioactiveDecay.hh & G4NuclearDecayChannel.hh: added private member data "applyICM" "applyARM" "halflifethresold"
     43  and their public Set methods SetICM(), SetARM() and SetHLThreshold().
     44- G4RadioactiveDecaymessenger: added the UICommands for SetICM, SETARM and SetHLThreshold.
     45
    174624 July 2009 V.Ivanchenko (radioactive_decay-V09-02-00)
    1847- GNUmakefile - added dependence on electromagnetic/utils
    1948
    204909 July 2008 Dennis Wright (radioactive_decay-V09-01-02)
     50--------------------------------------------------------
    2151- replace exit() with G4Exception in G4RadioactiveDecay and G4NuclearDecayChannel
    2252

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

    r819 r1315  
    8888  // Returns the decay products
    8989  //
     90  void SetHLThreshold (G4double hl) {halflifethreshold = hl;}
     91  // Set the half-life threshold for isomer production
     92  //
     93  void SetICM (G4bool icm) {applyICM = icm;}
     94  // Enable/disable ICM
     95  //
     96  void SetARM (G4bool arm) {applyARM = arm;}
     97  // Enable/disable ARM
     98  //
    9099  inline G4RadioactiveDecayMode GetDecayMode () {return decayMode;}
    91100  // Returns the decay mode
    … …  
    136145  G4double               FermiFN;
    137146  G4bool                 BetaSimple;
     147  G4double               halflifethreshold;
     148  G4bool                 applyICM;
     149  G4bool                 applyARM;
    138150  CLHEP::RandGeneral*    RandomEnergy;   
    139151};

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

    r819 r1315  
    114114  //   Sets the decay biasing scheme using the data in "filename"
    115115  //
     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  //
    116125  void SetSourceTimeProfile (G4String filename) ;
    117126  //  Sets source exposure function using histograms in "filename"
    … …  
    235244  G4int                         NSplit;
    236245
     246  G4double                      halflifethreshold;
     247  G4bool                        applyICM;
     248  G4bool                        applyARM;
     249
    237250  G4int                         NSourceBin;
    238251  G4double                      SBin[99];
    … …  
    243256  G4double                      DProfile[99];
    244257
    245   std::vector<G4String>              LoadedNuclei;
    246   std::vector<G4String>              ValidVolumes;
     258  std::vector<G4String>         LoadedNuclei;
     259  std::vector<G4String>         ValidVolumes;
    247260
    248261  G4RadioactiveDecayRate        theDecayRate;
    … …  
    262275  G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
    263276
    264   inline G4double AtRestGetPhysicalInteractionLengt
     277  inline G4double AtRestGetPhysicalInteractionLength
    265278    (const G4Track& track, G4ForceCondition* condition)
    266279  {fRemainderLifeTime = G4VRestDiscreteProcess::

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

    r819 r1315  
    107107  G4UIcmdWithoutParameter        *allvolumesCmd;
    108108  G4UIcmdWithoutParameter        *deallvolumesCmd;
     109  G4UIcmdWithABool               *icmCmd;
     110  G4UIcmdWithABool               *armCmd;
     111  G4UIcmdWithADoubleAndUnit      *hlthCmd;
    109112};
    110113

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

    r962 r1315  
    101101  FillDaughterNucleus (0, A, Z, theDaughterExcitation);
    102102  Qtransition = theQtransition;
     103  halflifethreshold = 1e-6*second;
     104  applyICM = true;
     105  applyARM = true;
    103106}
    104107///////////////////////////////////////////////////////////////////////////////
    … …  
    131134  FillDaughterNucleus (1, A, Z, theDaughterExcitation);
    132135  Qtransition = theQtransition;
     136  halflifethreshold = 1e-6*second;
     137  applyICM = true;
     138  applyARM = true;
    133139}
    134140///////////////////////////////////////////////////////////////////////////////
    … …  
    171177  Qtransition = theQtransition;
    172178  FermiFN = theFFN;
     179  halflifethreshold = 1e-6*second;
     180  applyICM = true;
     181  applyARM = true;
    173182}
    174183
    … …  
    288297  //
    289298  // If the decay is to an excited state of the daughter nuclide, we need
    290   // to apply the photo-evaporation process.
     299  // to apply the photo-evaporation process. This includes the IT decay mode itself.
    291300  //
    292301  // needed to hold the shell idex after ICM
    293302  G4int shellIndex = -1;
    294303  //
    295   if (daughterExcitation > 0.0)
     304  if (daughterExcitation > 0.0) 
    296305    {
    297306      //
    … …  
    311320      deexcitation->SetVerboseLevel(GetVerboseLevel());
    312321      // switch on/off internal electron conversion
    313       deexcitation->SetICM(true);
     322      deexcitation->SetICM(applyICM);
    314323      // set the maximum life-time for a level that will be treated. Level with life-time longer than this
    315324      // will be outputed as meta-stable isotope
    316325      //
    317       deexcitation->SetMaxHalfLife(1e-6*second);
     326      deexcitation->SetMaxHalfLife(halflifethreshold);
    318327      // but in IT mode, we need to force the transition
    319328      if (decayMode == 0) {
    … …  
    323332      }
    324333      //
    325       // Get the gammas by deexciting the nucleus.
    326       //
    327       G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus);
    328       // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide
     334      // Now apply the photo-evaporation
     335      // Use BreakUp() so limit to one transition at a time, if ICM is requested
     336      // this change is realted to bug#1001  (F.Lei 07/05/2010)
     337      G4FragmentVector* gammas = 0;     
     338      if (applyICM) {
     339        gammas = deexcitation->BreakUp(nucleus);       
     340      } else {
     341        gammas = deexcitation->BreakItUp(nucleus);
     342      }
     343      // the returned G4FragmentVector contains the residual nuclide
    329344      // as its last entry.
    330345      G4int nGammas=gammas->size()-1;
    … …  
    343358        }
    344359      //
    345       //      now the nucleus
     360      // now the nucleus
    346361      G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
    347362      // f.lei (03/01/03) this is needed to fix the crach in test18
    … …  
    356371         daughterMomentum1);
    357372      products->PushProducts (dynamicDaughter);
     373     
    358374      // retrive the ICM shell index
    359375      shellIndex = deexcitation->GetVacantShellNumber();
    … …  
    392408        //
    393409        {
    394           eShell = G4int(G4UniformRand()*5)+4;
     410        // limit the shell index to 6 as specified by the ARM (F.Lei 06/05/2010)
     411        // eShell = G4int(G4UniformRand()*5)+4;
     412          eShell = G4int(G4UniformRand()*3)+4;
    395413        }
    396414        break;
    … …  
    401419      }
    402420  }
    403   // now deal with the IT case where ICM may have been applied
     421  // Also have to deal with the IT case where ICM may have been applied
    404422  //
    405423  if (decayMode == 0) {
    406424    eShell = shellIndex;
    407425  }
    408   // now apply ARM if there is a vaccancy
    409   //
    410   if (eShell != -1) {
     426  //
     427  // now apply ARM if it is requested and there is a vaccancy
     428  //
     429  if (applyARM && eShell != -1) {
    411430    G4int aZ = daughterZ;
    412431    if (aZ > 5 && aZ < 100) {  // only applies to 5< Z <100

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

    r962 r1315  
    168168  FBeta       = false ;
    169169  BRBias      = true ;
     170  applyICM    = true ;
     171  applyARM    = true ;
     172  halflifethreshold = 1e-6*second;
    170173  //
    171174  // RDM applies to xall logical volumes as default
    … …  
    200203  // All particles, other than G4Ions, are rejected by default.
    201204  //
     205  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
    202206  if (aParticle.GetParticleName() == "GenericIon")      {return true;}
    203207  else if (!(aParticle.GetParticleType() == "nucleus") || aParticle.GetPDGLifeTime() < 0. ) {return false;}
    204 
    205208  //
    206209  //
    … …  
    482485    else if (theLife < 0.0) {meanlife = DBL_MAX;}
    483486    else {meanlife = theLife;}
     487  // set the meanlife to zero for excited istopes which is not in the RDM database
     488    if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
    484489  }
    485490#ifdef G4VERBOSE
    … …  
    612617//     theParentNucleus.
    613618//
    614 G4DecayTable *G4RadioactiveDecay::LoadDecayTable (G4ParticleDefinition
    615   &theParentNucleus)
     619G4DecayTable *G4RadioactiveDecay::LoadDecayTable (G4ParticleDefinition &theParentNucleus)
    616620{
    617621  //
    … …  
    645649 
    646650  std::ifstream DecaySchemeFile(file);
    647 
    648   if (!DecaySchemeFile)
     651  G4bool found(false);
     652  if (DecaySchemeFile) {
    649653    //
     654    //
     655    // Initialise variables used for reading in radioactive decay data.
     656    //
     657    G4int    nMode = 7;
     658    G4bool   modeFirstRecord[7];
     659    G4double modeTotalBR[7];
     660    G4double modeSumBR[7];
     661    G4int i;
     662    for (i=0; i<nMode; i++) {
     663      modeFirstRecord[i] = true;
     664      modeSumBR[i]       = 0.0;
     665    }
     666   
     667    G4bool complete(false);
     668    char inputChars[80]={' '};
     669    G4String inputLine;
     670    G4String recordType("");
     671    G4RadioactiveDecayMode theDecayMode;
     672    G4double a(0.0);
     673    G4double b(0.0);
     674    G4double c(0.0);
     675    G4double n(1.0);
     676    G4double e0;
     677    //
     678    //
     679    // Go through each record in the data file until you identify the decay
     680    // data relating to the nuclide of concern.
     681    //
     682    //  while (!complete && -DecaySchemeFile.getline(inputChars, 80).eof() != EOF)
     683    while (!complete && !DecaySchemeFile.getline(inputChars, 80).eof()) {
     684      inputLine = inputChars;
     685      //    G4String::stripType stripend(1);
     686      //    inputLine = inputLine.strip(stripend);
     687      inputLine = inputLine.strip(1);
     688      if (inputChars[0] != '#' && inputLine.length() != 0) {
     689        std::istringstream tmpStream(inputLine);
     690        if (inputChars[0] == 'P') {
     691          //
     692          //
     693          // Nucleus is a parent type.  Check the excitation level to see if it matches
     694          // that of theParentNucleus
     695          //
     696          tmpStream >>recordType >>a >>b;
     697          if (found) {complete = true;}
     698          else {found = (std::abs(a*keV - E)<levelTolerance);}
     699        } else if (found) {
     700          //
     701          //
     702          // The right part of the radioactive decay data file has been found.  Search
     703          // through it to determine the mode of decay of the subsequent records.
     704          //
     705          if (inputChars[0] == 'W') {
     706#ifdef G4VERBOSE
     707            if (GetVerboseLevel()>0) {
     708              // a comment line identified and print out the message
     709              //
     710              G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
     711              G4cout << "   In data file " << file << G4endl;
     712              G4cout << "   " << inputLine << G4endl;
     713            }
     714#endif
     715          }     
     716          else {
     717            tmpStream >>theDecayMode >>a >>b >>c;
     718            a/=1000.;
     719            c/=1000.;
     720                 
     721            //  cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
     722                 
     723            switch (theDecayMode) {
     724            case IT:
     725              {
     726                //
     727                //
     728                // Decay mode is isomeric transition.
     729                //
     730                G4ITDecayChannel *anITChannel = new G4ITDecayChannel
     731                  (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, b);
     732                anITChannel->SetICM(applyICM);
     733                anITChannel->SetARM(applyARM);
     734                anITChannel->SetHLThreshold(halflifethreshold);
     735                theDecayTable->Insert(anITChannel);
     736                break;
     737              }
     738            case BetaMinus:
     739              {
     740                //
     741                //
     742                // Decay mode is beta-.
     743                //
     744                if (modeFirstRecord[1])
     745                  {modeFirstRecord[1] = false; modeTotalBR[1] = b;}
     746                else {
     747                  if (c > 0.) {
     748                    // to work out the Fermi function normalization factor first
     749                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, (Z+1));
     750                    e0 = c*MeV/0.511;
     751                    n = aBetaFermiFunction->GetFFN(e0);
     752                               
     753                    // now to work out the histogram and initialise the random generator
     754                    G4int npti = 100;                           
     755                    G4double* pdf = new G4double[npti];
     756                    G4int ptn;
     757                    G4double g,e,ee,f;
     758                    ee = e0+1.;
     759                    for (ptn=0; ptn<npti; ptn++) {
     760                      // e =e0*(ptn+1.)/102.;
     761                      // bug fix (#662) (flei, 22/09/2004)
     762                      e =e0*(ptn+0.5)/100.;
     763                      g = e+1.;
     764                      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
     765                      // Special treatment for K-40 (problem #1068) (flei,06/05/2010)
     766                      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)));
     768                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
     769                    }             
     770                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
     771                    G4BetaMinusDecayChannel *aBetaMinusChannel = new
     772                      G4BetaMinusDecayChannel (GetVerboseLevel(), &theParentNucleus,
     773                                               b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
     774                    theDecayTable->Insert(aBetaMinusChannel);
     775                    modeSumBR[1] += b;
     776                    delete[] pdf;
     777                    delete aBetaFermiFunction;
     778                  }
     779                }
     780                break;
     781              case BetaPlus:
     782                //
     783                //
     784                // Decay mode is beta+.
     785                //
     786                if (modeFirstRecord[2])
     787                  {modeFirstRecord[2] = false; modeTotalBR[2] = b;}
     788                else {
     789                  //  e0 = c*MeV/0.511;
     790                  // bug fix (#662) (flei, 22/09/2004)
     791                  // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
     792                  // with Q < 2Me
     793                  //
     794                  e0 = c*MeV/0.511 -2.;
     795                  if (e0 > 0.) {
     796                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
     797                               
     798                    n = aBetaFermiFunction->GetFFN(e0);
     799                               
     800                    // now to work out the histogram and initialise the random generator
     801                    G4int npti = 100;                           
     802                    G4double* pdf = new G4double[npti];
     803                    G4int ptn;
     804                    G4double g,e,ee,f;
     805                    ee = e0+1.;
     806                    for (ptn=0; ptn<npti; ptn++) {
     807                      // e =e0*(ptn+1.)/102.;
     808                      // bug fix (#662) (flei, 22/09/2004)
     809                      e =e0*(ptn+0.5)/100.;
     810                      g = e+1.;
     811                      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
     812                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
     813                    }
     814                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
     815                    G4BetaPlusDecayChannel *aBetaPlusChannel = new
     816                      G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
     817                                              b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
     818                    theDecayTable->Insert(aBetaPlusChannel);
     819                    modeSumBR[2] += b;
     820                               
     821                    delete[] pdf;
     822                    delete aBetaFermiFunction;       
     823                  }
     824                }
     825                break;
     826              case KshellEC:
     827                //
     828                //
     829                // Decay mode is K-electron capture.
     830                //
     831                if (modeFirstRecord[3])
     832                  {modeFirstRecord[3] = false; modeTotalBR[3] = b;}
     833                else {
     834                  G4KshellECDecayChannel *aKECChannel = new
     835                    G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
     836                                            b, c*MeV, a*MeV);
     837                  aKECChannel->SetICM(applyICM);
     838                  aKECChannel->SetARM(applyARM);
     839                  aKECChannel->SetHLThreshold(halflifethreshold);
     840                  theDecayTable->Insert(aKECChannel);
     841                  //delete aKECChannel;
     842                  modeSumBR[3] += b;
     843                }
     844                break;
     845              case LshellEC:
     846                //
     847                //
     848                // Decay mode is L-electron capture.
     849                //
     850                if (modeFirstRecord[4])
     851                  {modeFirstRecord[4] = false; modeTotalBR[4] = b;}
     852                else {
     853                  G4LshellECDecayChannel *aLECChannel = new
     854                    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
     855                                            b, c*MeV, a*MeV);
     856                  aLECChannel->SetICM(applyICM);
     857                  aLECChannel->SetARM(applyARM);
     858                  aLECChannel->SetHLThreshold(halflifethreshold);
     859                  theDecayTable->Insert(aLECChannel);
     860                  //delete aLECChannel;
     861                  modeSumBR[4] += b;
     862                }
     863                break;
     864              case MshellEC:
     865                //
     866                //
     867                // Decay mode is M-electron capture. In this implementation it is added to L-shell case
     868                //
     869                if (modeFirstRecord[5])
     870                  {modeFirstRecord[5] = false; modeTotalBR[5] = b;}
     871                else {
     872                  G4MshellECDecayChannel *aMECChannel = new
     873                    G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
     874                                            b, c*MeV, a*MeV);
     875                  aMECChannel->SetICM(applyICM);
     876                  aMECChannel->SetARM(applyARM);
     877                  aMECChannel->SetHLThreshold(halflifethreshold);
     878                  theDecayTable->Insert(aMECChannel);
     879                  modeSumBR[5] += b;
     880                }
     881                break;
     882              case Alpha:
     883                //
     884                //
     885                // Decay mode is alpha.
     886                //
     887                if (modeFirstRecord[6])
     888                  {modeFirstRecord[6] = false; modeTotalBR[6] = b;}
     889                else {
     890                  G4AlphaDecayChannel *anAlphaChannel = new
     891                    G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
     892                                         b, c*MeV, a*MeV);
     893                  theDecayTable->Insert(anAlphaChannel);
     894                  //          delete anAlphaChannel;
     895                  modeSumBR[6] += b;
     896                }
     897                break;
     898              case ERROR:
     899              default:
     900                G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
     901                            FatalException, "Error in decay mode selection");
     902               
     903              }
     904            }
     905          }
     906        }
     907      }
     908    }   
     909    //
     910    //
     911    // Go through the decay table and make sure that the branching ratios are
     912    // correctly normalised.
     913    //
     914    G4VDecayChannel       *theChannel             = 0;
     915    G4NuclearDecayChannel *theNuclearDecayChannel = 0;
     916    G4String mode                     = "";
     917    G4int j                           = 0;
     918    G4double theBR                    = 0.0;
     919    for (i=0; i<theDecayTable->entries(); i++) {
     920      theChannel             = theDecayTable->GetDecayChannel(i);
     921      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
     922      theDecayMode           = theNuclearDecayChannel->GetDecayMode();
     923      j          = 0;
     924      if (theDecayMode != IT) {
     925        theBR = theChannel->GetBR();
     926        theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
     927      }
     928    }
     929  }             
     930  DecaySchemeFile.close();             
     931  if (!found && E > 0.) {
     932    // cases where IT cascade for exited isotopes without entry in RDM database
     933    // Decay mode is isomeric transition.
     934    //
     935    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
     936      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
     937    anITChannel->SetICM(applyICM);
     938    anITChannel->SetARM(applyARM);
     939    anITChannel->SetHLThreshold(halflifethreshold);
     940    theDecayTable->Insert(anITChannel);
     941  }
     942  if (!theDecayTable) {
    650943    //
    651944    // There is no radioactive decay data for this nucleus.  Return a null
    652945    // decay table.
    653946    //
    654     {
    655       G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
    656       theDecayTable = 0;
    657       return theDecayTable;
    658     }
    659   //
    660   //
    661   // Initialise variables used for reading in radioactive decay data.
    662   //
    663   G4int    nMode = 7;
    664   G4bool   modeFirstRecord[7];
    665   G4double modeTotalBR[7];
    666   G4double modeSumBR[7];
    667   G4int i;
    668   for (i=0; i<nMode; i++)
    669     {
    670       modeFirstRecord[i] = true;
    671       modeSumBR[i]       = 0.0;
    672     }
    673 
    674   G4bool complete(false);
    675   G4bool found(false);
    676   char inputChars[80]={' '};
    677   G4String inputLine;
    678   G4String recordType("");
    679   G4RadioactiveDecayMode theDecayMode;
    680   G4double a(0.0);
    681   G4double b(0.0);
    682   G4double c(0.0);
    683   G4double n(1.0);
    684   G4double e0;
    685   //
    686   //
    687   // Go through each record in the data file until you identify the decay
    688   // data relating to the nuclide of concern.
    689   //
    690 //  while (!complete && -DecaySchemeFile.getline(inputChars, 80).eof() != EOF)
    691   while (!complete && !DecaySchemeFile.getline(inputChars, 80).eof())    {
    692       inputLine = inputChars;
    693       //    G4String::stripType stripend(1);
    694       //    inputLine = inputLine.strip(stripend);
    695       inputLine = inputLine.strip(1);
    696       if (inputChars[0] != '#' && inputLine.length() != 0)
    697         {
    698           std::istringstream tmpStream(inputLine);
    699           if (inputChars[0] == 'P')
    700             //
    701             //
    702             // Nucleus is a parent type.  Check the excitation level to see if it matches
    703             // that of theParentNucleus
    704             //
    705             {
    706               tmpStream >>recordType >>a >>b;
    707               if (found) {complete = true;}
    708               else {found = (std::abs(a*keV - E)<levelTolerance);}
    709             }
    710           else if (found)
    711             {
    712               //
    713               //
    714               // The right part of the radioactive decay data file has been found.  Search
    715               // through it to determine the mode of decay of the subsequent records.
    716               //
    717               if (inputChars[0] == 'W') {
    718 #ifdef G4VERBOSE
    719                 if (GetVerboseLevel()>0) {
    720                   // a comment line identified and print out the message
    721                   //
    722                   G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
    723                   G4cout << "   In data file " << file << G4endl;
    724                   G4cout << "   " << inputLine << G4endl;
    725                 }
    726 #endif
    727               }
    728               else
    729                 {
    730                   tmpStream >>theDecayMode >>a >>b >>c;
    731                   a/=1000.;
    732                   c/=1000.;
    733          
    734                   //    cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
    735          
    736                   switch (theDecayMode)
    737                     {
    738                     case IT:
    739                       //
    740                       //
    741                       // Decay mode is isomeric transition.
    742                       //
    743                       {
    744                         G4ITDecayChannel *anITChannel = new G4ITDecayChannel
    745                           (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, b);
    746                         theDecayTable->Insert(anITChannel);
    747                         break;
    748                       }
    749                     case BetaMinus:
    750                       //
    751                       //
    752                       // Decay mode is beta-.
    753                       //
    754                       if (modeFirstRecord[1])
    755                         {modeFirstRecord[1] = false; modeTotalBR[1] = b;}
    756                       else
    757                         {
    758                           if (c > 0.) {
    759                             // to work out the Fermi function normalization factor first
    760                             G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, (Z+1));
    761                             e0 = c*MeV/0.511;
    762                             n = aBetaFermiFunction->GetFFN(e0);
    763                            
    764                             // now to work out the histogram and initialise the random generator
    765                             G4int npti = 100;                           
    766                             G4double* pdf = new G4double[npti];
    767                             G4int ptn;
    768                             G4double g,e,ee,f;
    769                             ee = e0+1.;
    770                             for (ptn=0; ptn<npti; ptn++) {
    771                               // e =e0*(ptn+1.)/102.;
    772                               // bug fix (#662) (flei, 22/09/2004)
    773                               e =e0*(ptn+0.5)/100.;
    774                               g = e+1.;
    775                               f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
    776                               pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
    777                             }             
    778                             RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
    779                             G4BetaMinusDecayChannel *aBetaMinusChannel = new
    780                               G4BetaMinusDecayChannel (GetVerboseLevel(), &theParentNucleus,
    781                                                        b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
    782                             theDecayTable->Insert(aBetaMinusChannel);
    783                             modeSumBR[1] += b;
    784                             delete[] pdf;
    785                             delete aBetaFermiFunction;
    786                           }
    787                         }
    788                       break;
    789                     case BetaPlus:
    790                       //
    791                       //
    792                       // Decay mode is beta+.
    793                       //
    794                       if (modeFirstRecord[2])
    795                         {modeFirstRecord[2] = false; modeTotalBR[2] = b;}
    796                       else
    797                         {
    798                           //  e0 = c*MeV/0.511;
    799                           // bug fix (#662) (flei, 22/09/2004)
    800                           // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
    801                           // with Q < 2Me
    802                           //
    803                           e0 = c*MeV/0.511 -2.;
    804                           if (e0 > 0.) {
    805                             G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
    806                            
    807                             n = aBetaFermiFunction->GetFFN(e0);
    808 
    809                             // now to work out the histogram and initialise the random generator
    810                             G4int npti = 100;                           
    811                             G4double* pdf = new G4double[npti];
    812                             G4int ptn;
    813                             G4double g,e,ee,f;
    814                             ee = e0+1.;
    815                             for (ptn=0; ptn<npti; ptn++) {
    816                               // e =e0*(ptn+1.)/102.;
    817                               // bug fix (#662) (flei, 22/09/2004)
    818                               e =e0*(ptn+0.5)/100.;
    819                               g = e+1.;
    820                               f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
    821                               pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
    822                             }
    823                             RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
    824                             G4BetaPlusDecayChannel *aBetaPlusChannel = new
    825                               G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
    826                                                       b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
    827                             theDecayTable->Insert(aBetaPlusChannel);
    828                             modeSumBR[2] += b;
    829                            
    830                             delete[] pdf;
    831                             delete aBetaFermiFunction;       
    832                           }
    833                         }
    834                       break;
    835                     case KshellEC:
    836                       //
    837                       //
    838                       // Decay mode is K-electron capture.
    839                       //
    840                       if (modeFirstRecord[3])
    841                         {modeFirstRecord[3] = false; modeTotalBR[3] = b;}
    842                       else
    843                         {
    844                           G4KshellECDecayChannel *aKECChannel = new
    845                             G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
    846                                                     b, c*MeV, a*MeV);
    847                           theDecayTable->Insert(aKECChannel);
    848                           //delete aKECChannel;
    849                           modeSumBR[3] += b;
    850                         }
    851                       break;
    852                     case LshellEC:
    853                       //
    854                       //
    855                       // Decay mode is L-electron capture.
    856                       //
    857                       if (modeFirstRecord[4])
    858                         {modeFirstRecord[4] = false; modeTotalBR[4] = b;}
    859                       else
    860                         {
    861                           G4LshellECDecayChannel *aLECChannel = new
    862                             G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
    863                                                     b, c*MeV, a*MeV);
    864                           theDecayTable->Insert(aLECChannel);
    865                           //delete aLECChannel;
    866                           modeSumBR[4] += b;
    867                         }
    868                       break;
    869                     case MshellEC:
    870                       //
    871                       //
    872                       // Decay mode is M-electron capture. In this implementation it is added to L-shell case
    873                       //
    874                       if (modeFirstRecord[5])
    875                         {modeFirstRecord[5] = false; modeTotalBR[5] = b;}
    876                       else
    877                         {
    878                           G4MshellECDecayChannel *aMECChannel = new
    879                             G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
    880                                                     b, c*MeV, a*MeV);
    881                           theDecayTable->Insert(aMECChannel);
    882                           modeSumBR[5] += b;
    883                         }
    884                       break;
    885                     case Alpha:
    886                       //
    887                       //
    888                       // Decay mode is alpha.
    889                       //
    890                       if (modeFirstRecord[6])
    891                         {modeFirstRecord[6] = false; modeTotalBR[6] = b;}
    892                       else
    893                         {
    894                           G4AlphaDecayChannel *anAlphaChannel = new
    895                             G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
    896                                                  b, c*MeV, a*MeV);
    897                           theDecayTable->Insert(anAlphaChannel);
    898                           //          delete anAlphaChannel;
    899                           modeSumBR[6] += b;
    900                         }
    901                       break;
    902                     case ERROR:
    903                     default:
    904                       G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
    905                                   FatalException, "Error in decay mode selection");
    906                      
    907                     }
    908                 }
    909             }
    910         }
    911 
    912     } 
    913   DecaySchemeFile.close();
    914   if (!found && E > 0.) {
    915     // cases where IT cascade follows immediately after a decay
    916    
    917     // Decay mode is isomeric transition.
    918     //
    919     G4ITDecayChannel *anITChannel = new G4ITDecayChannel
    920       (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
    921     theDecayTable->Insert(anITChannel);
    922   }
    923   //
    924   //
    925   // Go through the decay table and make sure that the branching ratios are
    926   // correctly normalised.
    927   //
    928   G4VDecayChannel       *theChannel             = 0;
    929   G4NuclearDecayChannel *theNuclearDecayChannel = 0;
    930   G4String mode                     = "";
    931   G4int j                           = 0;
    932   G4double theBR                    = 0.0;
    933   for (i=0; i<theDecayTable->entries(); i++)
    934     {
    935       theChannel             = theDecayTable->GetDecayChannel(i);
    936       theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
    937       theDecayMode           = theNuclearDecayChannel->GetDecayMode();
    938       j          = 0;
    939       if (theDecayMode != IT)
    940         {
    941           theBR = theChannel->GetBR();
    942           theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
    943         }
    944     } 
    945 
     947    G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
     948    theDecayTable = 0;
     949    return theDecayTable;
     950  }     
    946951  if (theDecayTable && GetVerboseLevel()>1)
    947952    {
    … …  
    10861091          if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
    10871092           
    1088             // Level hafe life is in ns and I want to set the gate as 1 micros
    1089             if ( theDecayMode == 0 && level->HalfLife() >= 1000.){   
     1093            // Level hafe life is in ns and the gate as 1 micros by default
     1094            if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){   
    10901095              // some further though may needed here
    10911096              theDecayTable->Insert(theChannel);
    … …  
    14241429        G4DecayProducts *products = DoDecay(*theParticleDef);
    14251430        //
     1431        // check if the product is the same as input and kill the track if necessary to prevent infinite loop
     1432        // (11/05/10, F.Lei)
     1433        //
     1434        if ( products->entries() == 1) {
     1435          fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
     1436          fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
     1437          fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
     1438          ClearNumberOfInteractionLengthLeft();
     1439          return &fParticleChangeForRadDecay;
     1440        }       
    14261441        //
    14271442        // Get parent particle information and boost the decay products to the
    … …  
    14401455            // time lapsed between the particle come to rest and the actual decay. This time
    14411456            // is simply sampled with the mean-life of the particle.
    1442             //
    1443             finalGlobalTime += -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime() ;
     1457            // but we need to protect the case PDGTime < 0. (F.Lei 11/05/10)
     1458            G4double temptime =  -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime();
     1459            if (temptime <0.) temptime =0.;
     1460            finalGlobalTime += temptime ;
    14441461            energyDeposit += theParticle->GetKineticEnergy();
    14451462          }

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

    r819 r1315  
    116116  brbiasCmd->SetParameterName("BRBias",true);
    117117  brbiasCmd->SetDefaultValue(true);
    118 
     118  //
     119  // Command contols whether ICM will be applied or not
     120  //
     121  icmCmd = new G4UIcmdWithABool ("/grdm/applyICM",this);
     122  icmCmd->SetGuidance("True: ICM is applied; false: no");
     123  icmCmd->SetParameterName("applyICM",true);
     124  icmCmd->SetDefaultValue(true);
     125  //
     126  // Command contols whether ARM will be applied or not
     127  //
     128  armCmd = new G4UIcmdWithABool ("/grdm/applyARM",this);
     129  armCmd->SetGuidance("True: ARM is applied; false: no");
     130  armCmd->SetParameterName("applyARM",true);
     131  armCmd->SetDefaultValue(true);
     132  //
     133  // Command to set the h-l thresold for isomer production
     134  //
     135  hlthCmd = new G4UIcmdWithADoubleAndUnit("/grdm/hlThreshold",this);
     136  hlthCmd->SetGuidance("Set the h-l threshold for isomer production");
     137  hlthCmd->SetParameterName("hlThreshold",false);
     138  hlthCmd->SetRange("hlThreshold>0.");
     139  hlthCmd->SetUnitCategory("Time");
     140  hlthCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
    119141  //
    120142  // Command to define the incident particle source time profile.
    … …  
    174196  delete allvolumesCmd;
    175197  delete deallvolumesCmd;
     198  delete icmCmd;
     199  delete armCmd;
     200  delete hlthCmd;
    176201}
    177202////////////////////////////////////////////////////////////////////////////////
    … …  
    215240  else if (command==verboseCmd) {theRadioactiveDecayContainer->
    216241                                       SetVerboseLevel(verboseCmd->GetNewIntValue(newValues));}
     242  else if (command==icmCmd ) {theRadioactiveDecayContainer->
     243      SetICM(icmCmd->GetNewBoolValue(newValues));}
     244  else if (command==armCmd ) {theRadioactiveDecayContainer->
     245      SetARM(armCmd->GetNewBoolValue(newValues));}
     246  else if (command==hlthCmd ) {theRadioactiveDecayContainer->
     247      SetHLThreshold(hlthCmd->GetNewDoubleValue(newValues));}
    217248}
    218249
Note: See TracChangeset for help on using the changeset viewer.