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

geant4 tag 9.4

File:
1 edited

Legend:

Unmodified
Added
Removed

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