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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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          }
Note: See TracChangeset for help on using the changeset viewer.