Ignore:
Timestamp:
Nov 25, 2009, 5:13:58 PM (15 years ago)
Author:
garnier
Message:

update CVS release candidate geant4.9.3.01

Location:
trunk/source/processes/cuts/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/cuts/src/G4MCCIndexConversionTable.cc

    r1007 r1196  
    2525//
    2626// $Id: G4MCCIndexConversionTable.cc,v 1.3 2006/06/29 19:30:08 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//

  • trunk/source/processes/cuts/src/G4MaterialCutsCouple.cc

    r1007 r1196  
    2626//
    2727// $Id: G4MaterialCutsCouple.cc,v 1.3 2006/06/29 19:30:10 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//

  • trunk/source/processes/cuts/src/G4PhysicsTableHelper.cc

    r1007 r1196  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PhysicsTableHelper.cc,v 1.5 2006/06/29 19:30:12 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4PhysicsTableHelper.cc,v 1.6 2009/08/01 07:57:13 kurasige Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//
    … …  
    7373      // enlarge physcis table
    7474      physTable->resize(numberOfMCC, (G4PhysicsVector*)(0));
     75#ifdef G4VERBOSE 
     76      if (verboseLevel>2) {
     77        G4cerr << "G4PhysicsTableHelper::PreparePhysicsTable  ";
     78        G4cerr << "Physics Table "<< physTable ;
     79        G4cerr << " is resized to " << numberOfMCC << G4endl;
     80      }
     81#endif
    7582    } else if ( physTable->size() > numberOfMCC){
    7683      // ERROR: this situation should not occur 
    … …  
    162169      size_t i = converter->GetIndex(idx);
    163170      G4PhysicsVector* vec = (*physTable)[i];
    164       if (vec !=0 ) delete vec;
     171       if (vec !=0 ) delete vec;
    165172      (*physTable)[i] =  (*tempTable)[idx];
    166173      physTable->ClearFlag(i);

  • trunk/source/processes/cuts/src/G4ProductionCuts.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4ProductionCuts.cc,v 1.5 2006/06/29 19:30:14 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4ProductionCuts.cc,v 1.6 2009/08/01 07:57:13 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    4040const G4ParticleDefinition* G4ProductionCuts::electDef = 0;
    4141const G4ParticleDefinition* G4ProductionCuts::positDef = 0;
     42const G4ParticleDefinition* G4ProductionCuts::protonDef = 0;
    4243
    4344G4ProductionCuts::G4ProductionCuts() :
    … …  
    8788G4int  G4ProductionCuts::GetIndex(const G4String& name)
    8889{
     90  static G4String gamma("gamma");
     91  static G4String electron("e-");
     92  static G4String positron("e+");
     93  static G4String proton("proton");
     94 
    8995  G4int index;
    90   if       ( name == "gamma" )        { index =  0; }
    91   else  if ( name == "e-" )           { index =  1; }
    92   else  if ( name == "e+" )           { index =  2; }
    93   else                                { index = -1; }
     96  if       ( name == gamma )        { index =  0; }
     97  else  if ( name == electron )     { index =  1; }
     98  else  if ( name == positron )     { index =  2; }
     99  else  if ( name == proton )       { index =  3; }
     100  else                              { index = -1; }
    94101
    95102  return index;
    … …  
    100107{
    101108  if(!ptcl) return -1;
    102   if(gammaDef==0 && ptcl->GetParticleName()=="gamma") { gammaDef = ptcl; }
    103   if(electDef==0 && ptcl->GetParticleName()=="e-") { electDef = ptcl; }
    104   if(positDef==0 && ptcl->GetParticleName()=="e+") { positDef = ptcl; }
     109  // In the first call, pointers are set
     110  if(gammaDef==0  && ptcl->GetParticleName()=="gamma")  { gammaDef = ptcl; }
     111  if(electDef==0  && ptcl->GetParticleName()=="e-")     { electDef = ptcl; }
     112  if(positDef==0  && ptcl->GetParticleName()=="e+")     { positDef = ptcl; }
     113  if(protonDef==0 && ptcl->GetParticleName()=="proton") { protonDef = ptcl; }
     114
    105115  G4int index;
    106   if(ptcl==gammaDef)      { index = 0; }
    107   else if(ptcl==electDef) { index = 1; }
    108   else if(ptcl==positDef) { index = 2; }
    109   else                    { index = -1; }
     116  if(ptcl==gammaDef)       { index = 0;  }
     117  else if(ptcl==electDef)  { index = 1;  }
     118  else if(ptcl==positDef)  { index = 2;  }
     119  else if(ptcl==protonDef) { index = 3;  }
     120  else                     { index = -1; }
    110121
    111122  return index;

  • trunk/source/processes/cuts/src/G4ProductionCutsTable.cc

    r1055 r1196  
    2525//
    2626//
    27 // $Id: G4ProductionCutsTable.cc,v 1.19 2009/04/02 02:43:42 kurasige Exp $
    28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     27// $Id: G4ProductionCutsTable.cc,v 1.25 2009/11/11 03:20:22 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    4747#include "G4RToEConvForGamma.hh"
    4848#include "G4RToEConvForPositron.hh"
     49#include "G4RToEConvForProton.hh"
    4950#include "G4MaterialTable.hh"
    5051#include "G4Material.hh"
    5152#include "G4UnitsTable.hh"
    5253
     54#include "G4Timer.hh"
     55
    5356#include "G4ios.hh"
    5457#include <iomanip>               
    … …  
    5861G4ProductionCutsTable* G4ProductionCutsTable::fG4ProductionCutsTable = 0;
    5962
     63/////////////////////////////////////////////////////////////
    6064G4ProductionCutsTable* G4ProductionCutsTable::GetProductionCutsTable()
    6165{
    … …  
    6771}
    6872
     73/////////////////////////////////////////////////////////////
    6974G4ProductionCutsTable::G4ProductionCutsTable()
    7075  : firstUse(true),verboseLevel(1),fMessenger(0)
    … …  
    8590}
    8691
     92/////////////////////////////////////////////////////////////
    8793G4ProductionCutsTable::G4ProductionCutsTable(const G4ProductionCutsTable& )
    8894{;}
    8995
     96/////////////////////////////////////////////////////////////
    9097G4ProductionCutsTable::~G4ProductionCutsTable()
    9198{
    … …  
    127134      converters[2] = new G4RToEConvForPositron();
    128135      converters[2]->SetVerboseLevel(GetVerboseLevel());
     136    }
     137    if(G4ParticleTable::GetParticleTable()->FindParticle("proton")){
     138      converters[3] = new G4RToEConvForProton();
     139      converters[3]->SetVerboseLevel(GetVerboseLevel());
    129140    }
    130141    firstUse = false;
    … …  
    217228  // Update RangeEnergy cuts tables
    218229  size_t idx = 0;
     230  G4Timer timer;
     231  if (verboseLevel>2) {
     232    timer.Start();
     233  }
    219234  for(CoupleTableIterator cItr=coupleTable.begin();
    220235      cItr!=coupleTable.end();cItr++){
    … …  
    235250    idx++; 
    236251  }
     252  if (verboseLevel>2) {
     253    timer.Stop();
     254    std::cout << "G4ProductionCutsTable::UpdateCoupleTable "
     255              << "  elapsed time for calculation of  energy cuts " << G4endl;
     256    std::cout << timer <<G4endl;
     257  }
    237258
    238259  // resize Range/Energy cuts double vectors if new couple is made
    … …  
    258279
    259280
     281/////////////////////////////////////////////////////////////
    260282G4double G4ProductionCutsTable::ConvertRangeToEnergy(
    261                                                       const G4ParticleDefinition* particle,
    262                                                       const G4Material*           material,
    263                                                       G4double                    range     
     283                                  const G4ParticleDefinition* particle,
     284                                  const G4Material*           material,
     285                                  G4double                    range     
    264286                                                     )
    265287{
    … …  
    272294
    273295  // check particle
    274   G4int index = -1;
    275   if (particle->GetParticleName() == "gamma") index = 0;
    276   if (particle->GetParticleName() == "e-")    index = 1;
    277   if (particle->GetParticleName() == "e+")    index = 2;
     296  G4int index = G4ProductionCuts::GetIndex(particle);
     297
    278298  if (index<0) {
    279299#ifdef G4VERBOSE 
    … …  
    290310}
    291311
    292 
    293 
     312/////////////////////////////////////////////////////////////
     313void G4ProductionCutsTable::ResetConverters()
     314{
     315  for(size_t i=0;i< NumberOfG4CutIndex;i++){
     316    if (converters[i]!=0) converters[i]->Reset();
     317  }
     318}
     319
     320
     321/////////////////////////////////////////////////////////////
    294322void G4ProductionCutsTable::SetEnergyRange(G4double lowedge, G4double highedge)
    295323{
    … …  
    297325}
    298326
     327/////////////////////////////////////////////////////////////
    299328G4double  G4ProductionCutsTable::GetLowEdgeEnergy() const
    300329{
    … …  
    302331}
    303332
     333/////////////////////////////////////////////////////////////
    304334G4double G4ProductionCutsTable::GetHighEdgeEnergy() const
    305335{
    … …  
    308338 
    309339
     340/////////////////////////////////////////////////////////////
    310341void G4ProductionCutsTable::ScanAndSetCouple(G4LogicalVolume* aLV,
    311342                                             G4MaterialCutsCouple* aCouple,
    … …  
    330361}
    331362
     363/////////////////////////////////////////////////////////////
    332364void G4ProductionCutsTable::DumpCouples() const
    333365{
    … …  
    350382    G4cout << " Material : " << aCouple->GetMaterial()->GetName() << G4endl;
    351383    G4cout << " Range cuts        : "
    352            << " gamma " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
    353            << "    e- " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
    354            << "    e+ " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
     384           << " gamma  " << G4BestUnit(aCut->GetProductionCut("gamma"),"Length")
     385           << "    e-  " << G4BestUnit(aCut->GetProductionCut("e-"),"Length")
     386           << "    e+  " << G4BestUnit(aCut->GetProductionCut("e+"),"Length")
     387           << " proton " << G4BestUnit(aCut->GetProductionCut("proton"),"Length")
    355388           << G4endl;
    356389    G4cout << " Energy thresholds : " ;
    … …  
    358391      G4cout << " is not ready to print";
    359392    } else {
    360       G4cout << " gamma " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
    361              << "    e- " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
    362              << "    e+ " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy");
     393      G4cout << " gamma  " << G4BestUnit((*(energyCutTable[0]))[aCouple->GetIndex()],"Energy")
     394             << "    e-  " << G4BestUnit((*(energyCutTable[1]))[aCouple->GetIndex()],"Energy")
     395             << "    e+  " << G4BestUnit((*(energyCutTable[2]))[aCouple->GetIndex()],"Energy")
     396             << " proton " << G4BestUnit((*(energyCutTable[3]))[aCouple->GetIndex()],"Energy");
    363397    }
    364398    G4cout << G4endl;
    … …  
    381415           
    382416
     417/////////////////////////////////////////////////////////////
    383418// Store cuts and material information in files under the specified directory.
    384419G4bool  G4ProductionCutsTable::StoreCutsTable(const G4String& dir,
    … …  
    404439}
    405440 
     441/////////////////////////////////////////////////////////////
    406442G4bool  G4ProductionCutsTable::RetrieveCutsTable(const G4String& dir,
    407443                                                 G4bool          ascii)
    … …  
    424460}
    425461
     462/////////////////////////////////////////////////////////////
    426463// check stored material and cut values are consistent
    427464// with the current detector setup.
    … …  
    444481}
    445482 
     483/////////////////////////////////////////////////////////////
    446484// Store material information in files under the specified directory.
    447485//
    … …  
    450488{
    451489  const G4String fileName = directory + "/" + "material.dat";
    452   const G4String key = "MATERIAL-V2.0";
     490  const G4String key = "MATERIAL-V3.0";
    453491  std::ofstream fOut; 
    454492
    … …  
    528566}
    529567
     568/////////////////////////////////////////////////////////////
    530569// check stored material is consistent with the current detector setup.
    531570G4bool  G4ProductionCutsTable::CheckMaterialInfo(const G4String& directory,
    … …  
    533572{
    534573  const G4String fileName = directory + "/" + "material.dat";
    535   const G4String key = "MATERIAL-V2.0";
     574  const G4String key = "MATERIAL-V3.0";
    536575  std::ifstream fIn; 
    537576
    … …  
    648687
    649688 
     689/////////////////////////////////////////////////////////////
    650690// Store materialCutsCouple information in files under the specified directory.
    651691//
    … …  
    655695{ 
    656696  const G4String fileName = directory + "/" + "couple.dat";
    657   const G4String key = "COUPLE-V2.0";
     697  const G4String key = "COUPLE-V3.0";
    658698  std::ofstream fOut; 
    659699  char temp[FixedStringLengthForStore];
    … …  
    774814
    775815
     816/////////////////////////////////////////////////////////////
    776817// check stored materialCutsCouple is consistent
    777818// with the current detector setup.
    … …  
    782823{
    783824  const G4String fileName = directory + "/" + "couple.dat";
    784   const G4String key = "COUPLE-V2.0";
     825  const G4String key = "COUPLE-V3.0";
    785826  std::ifstream fIn; 
    786827
    … …  
    944985
    945986
     987/////////////////////////////////////////////////////////////
    946988// Store cut values information in files under the specified directory.
    947989//
    … …  
    950992{
    951993  const G4String fileName = directory + "/" + "cut.dat";
    952   const G4String key = "CUT-V2.0";
     994  const G4String key = "CUT-V3.0";
    953995  std::ofstream fOut; 
    954996  char temp[FixedStringLengthForStore];
    … …  
    10181060}
    10191061 
     1062/////////////////////////////////////////////////////////////
    10201063// Retrieve cut values information in files under the specified directory.
    10211064//
    … …  
    10241067{
    10251068  const G4String fileName = directory + "/" + "cut.dat";
    1026   const G4String key = "CUT-V2.0";
     1069  const G4String key = "CUT-V3.0";
    10271070  std::ifstream fIn; 
    10281071
    … …  
    10951138}
    10961139
     1140/////////////////////////////////////////////////////////////
    10971141// Set Verbose Level
    10981142//   set same verbosity to all registered RangeToEnergyConverters 
    … …  
    11071151}
    11081152
     1153/////////////////////////////////////////////////////////////
     1154G4double G4ProductionCutsTable::GetMaxEnergyCut()
     1155{
     1156  return G4VRangeToEnergyConverter::GetMaxEnergyCut();
     1157}
     1158
     1159
     1160///////////////////////////////////////////////////////////// 
     1161void G4ProductionCutsTable::SetMaxEnergyCut(G4double value)
     1162{
     1163  G4VRangeToEnergyConverter::SetMaxEnergyCut(value);
     1164}

  • trunk/source/processes/cuts/src/G4ProductionCutsTableMessenger.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4ProductionCutsTableMessenger.cc,v 1.1 2008/03/02 10:52:55 kurasige Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4ProductionCutsTableMessenger.cc,v 1.3 2009/11/12 00:20:03 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    8686  setHighEdgeCmd->AvailableForStates(G4State_PreInit);
    8787 
     88  // /cuts/setMaxCutEnergy command
     89  setMaxEnergyCutCmd = new G4UIcmdWithADoubleAndUnit("/cuts/setMaxCutEnergy",this);
     90  setMaxEnergyCutCmd->SetGuidance("Set maximum of cut energy value ");
     91  setMaxEnergyCutCmd->SetParameterName("cut",false);
     92  setMaxEnergyCutCmd->SetDefaultValue(10.0);
     93  setMaxEnergyCutCmd->SetRange("cut >0.0");
     94  setMaxEnergyCutCmd->SetDefaultUnit("GeV");
     95  setMaxEnergyCutCmd->AvailableForStates(G4State_PreInit);
     96 
    8897 // /cuts/dump command
    8998  dumpCmd = new G4UIcmdWithoutParameter("/cuts/dump",this);
    … …  
    95104{
    96105  delete dumpCmd;
     106  delete setMaxEnergyCutCmd;
    97107  delete setHighEdgeCmd;
    98108  delete setLowEdgeCmd;
    … …  
    120130    theCutsTable->SetEnergyRange(lowEdge, highEdge);
    121131   
     132 } else if( command==setMaxEnergyCutCmd ){
     133    G4double cut = setHighEdgeCmd->GetNewDoubleValue(newValue);
     134    theCutsTable->SetMaxEnergyCut(cut);
     135   
    122136  }
    123137}
    … …  
    137151    G4double highEdge = theCutsTable->GetHighEdgeEnergy();
    138152    cv = setHighEdgeCmd->ConvertToString( highEdge, "TeV" );
     153
     154 } else if( command==setMaxEnergyCutCmd ){
     155    G4double cut = theCutsTable->GetMaxEnergyCut();
     156    cv = setMaxEnergyCutCmd->ConvertToString( cut, "GeV" );
    139157 }
     158
    140159
    141160  return cv;

  • trunk/source/processes/cuts/src/G4RToEConvForElectron.cc

    r1055 r1196  
    2525//
    2626//
    27 // $Id: G4RToEConvForElectron.cc,v 1.6 2009/04/02 02:43:42 kurasige Exp $
    28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     27// $Id: G4RToEConvForElectron.cc,v 1.7 2009/09/11 15:21:39 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    7171  const  G4double Tlow=10.*keV, Thigh=1.*GeV;
    7272  static G4double bremfactor= 0.1 ;
    73  
    74   G4double Mass = theParticle->GetPDGMass();       
     73
     74  static G4double Mass= theParticle->GetPDGMass();
     75
    7576  //  calculate dE/dx for electrons
    7677  if( std::fabs(AtomicNumber-Z)>0.1 ) {
    … …  
    119120  return dEdx;
    120121}
    121 
    122 
    123 void G4RToEConvForElectron::BuildRangeVector(const G4Material* aMaterial,
    124                                              G4double       maxEnergy,
    125                                              G4double       aMass,
    126                                              G4PhysicsLogVector* rangeVector)
    127 {
    128   //  create range vector for a material
    129   const G4double tlim = 10.*keV;
    130   const G4int maxnbint = 100;
    131 
    132   const G4ElementVector* elementVector = aMaterial->GetElementVector();
    133   const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
    134   G4int NumEl = aMaterial->GetNumberOfElements();
    135 
    136   // calculate parameters of the low energy part first
    137   size_t i;
    138   G4double loss=0.;
    139   for (i=0; i<size_t(NumEl); i++) {
    140     G4bool isOut;
    141     G4int IndEl = (*elementVector)[i]->GetIndex();
    142     loss += atomicNumDensityVector[i]*
    143            (*theLossTable)[IndEl]->GetValue(tlim,isOut);
    144   }
    145   G4double taulim = tlim/aMass;
    146   G4double clim = std::sqrt(taulim)*loss;
    147   G4double taumax = maxEnergy/aMass;
    148 
    149   // now the range vector can be filled
    150   for ( i=0; i<size_t(TotBin); i++) {
    151     G4double LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
    152     G4double tau = LowEdgeEnergy/aMass;
    153 
    154     if ( tau <= taulim ) {
    155       G4double Value = 2.*aMass*tau*std::sqrt(tau)/(3.*clim);
    156       rangeVector->PutValue(i,Value);
    157     } else {
    158       G4double rangelim = 2.*aMass*taulim*std::sqrt(taulim)/(3.*clim);
    159       G4double ltaulow = std::log(taulim);
    160       G4double ltauhigh = std::log(tau);
    161       G4double ltaumax = std::log(taumax);
    162       G4int    nbin = G4int(maxnbint*(ltauhigh-ltaulow)/(ltaumax-ltaulow));
    163       if( nbin < 1 ) nbin = 1;
    164       G4double Value = RangeLogSimpson( NumEl, elementVector,
    165                                         atomicNumDensityVector, aMass,
    166                                         ltaulow, ltauhigh, nbin)
    167                      + rangelim;
    168       rangeVector->PutValue(i,Value);
    169     }
    170   }
    171 }

  • trunk/source/processes/cuts/src/G4RToEConvForGamma.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4RToEConvForGamma.cc,v 1.4 2006/06/29 19:30:24 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4RToEConvForGamma.cc,v 1.6 2009/09/12 12:09:42 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    5353#endif
    5454  }
    55   TotBin = 100;
    5655}
    5756
    … …  
    6665void G4RToEConvForGamma::BuildAbsorptionLengthVector(
    6766                            const G4Material* aMaterial,
    68                             G4double       ,     
    69                             G4double       ,
    7067                            G4RangeVector* absorptionLengthVector )
    7168{
    … …  
    8279  G4int NumEl = aMaterial->GetNumberOfElements();
    8380  G4double absorptionLengthMax = 0.0;
     81  G4double previous = 0.;
    8482  for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
    85     G4double lowEdgeEnergy = absorptionLengthVector->GetLowEdgeEnergy(ibin);
    86    
    8783    G4double SIGMA = 0. ;
    88    
    8984    for (size_t iel=0; iel<size_t(NumEl); iel++) {
    90       G4bool isOut;
    9185      G4int IndEl = (*elementVector)[iel]->GetIndex();
    9286      SIGMA +=  atomicNumDensityVector[iel]*
    93                    (*aCrossSectionTable)[IndEl]->GetValue(lowEdgeEnergy,isOut);
     87                   (*((*aCrossSectionTable)[IndEl]))[ibin];
    9488    }
    9589    //  absorption length=5./SIGMA
    9690    absorptionLengthVector->PutValue(ibin, 5./SIGMA);
    9791    if (absorptionLengthMax < 5./SIGMA ) absorptionLengthMax = 5./SIGMA;
     92 
     93    //if (previous > 5./SIGMA) {
     94    //  G4cout << "G4RToEConvForGamma::BuildAbsorptionLengthVector"
     95    //     << ": WARNING absorptionVector "
     96    //     << ibin << ":" <<  5./SIGMA  << " <-- "
     97    //     << ibin -1 <<  ":" << previous <<G4endl;
     98    //}
     99
     100    previous = 5./SIGMA;
    98101  }
    99102}

  • trunk/source/processes/cuts/src/G4RToEConvForPositron.cc

    r1055 r1196  
    2525//
    2626//
    27 // $Id: G4RToEConvForPositron.cc,v 1.6 2009/04/02 02:43:42 kurasige Exp $
    28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     27// $Id: G4RToEConvForPositron.cc,v 1.7 2009/09/11 15:21:39 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    122122
    123123
    124 
    125 void G4RToEConvForPositron::BuildRangeVector(const G4Material* aMaterial,
    126                                              G4double       maxEnergy,
    127                                              G4double       aMass,
    128                                              G4PhysicsLogVector* rangeVector)
    129 {
    130   //  create range vector for a material
    131   const G4double tlim = 10.*keV;
    132   const G4int maxnbint = 100;
    133 
    134   const G4ElementVector* elementVector = aMaterial->GetElementVector();
    135   const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
    136   G4int NumEl = aMaterial->GetNumberOfElements();
    137 
    138   // calculate parameters of the low energy part first
    139   size_t i;
    140   G4double loss=0.;
    141   for (i=0; i<size_t(NumEl); i++) {
    142     G4bool isOut;
    143     G4int IndEl = (*elementVector)[i]->GetIndex();
    144     loss += atomicNumDensityVector[i]*
    145            (*theLossTable)[IndEl]->GetValue(tlim,isOut);
    146   }
    147   G4double taulim = tlim/aMass;
    148   G4double clim = std::sqrt(taulim)*loss;
    149   G4double taumax = maxEnergy/aMass;
    150 
    151   // now the range vector can be filled
    152   for ( i=0; i<size_t(TotBin); i++) {
    153     G4double LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
    154     G4double tau = LowEdgeEnergy/aMass;
    155 
    156     if ( tau <= taulim ) {
    157       G4double Value = 2.*aMass*tau*std::sqrt(tau)/(3.*clim);
    158       rangeVector->PutValue(i,Value);
    159     } else {
    160       G4double rangelim = 2.*aMass*taulim*std::sqrt(taulim)/(3.*clim);
    161       G4double ltaulow = std::log(taulim);
    162       G4double ltauhigh = std::log(tau);
    163       G4double ltaumax = std::log(taumax);
    164       G4int    nbin = G4int(maxnbint*(ltauhigh-ltaulow)/(ltaumax-ltaulow));
    165       if( nbin < 1 ) nbin = 1;
    166       G4double Value = RangeLogSimpson( NumEl, elementVector,
    167                                         atomicNumDensityVector, aMass,
    168                                         ltaulow,       ltauhigh, nbin)
    169                      + rangelim;
    170       rangeVector->PutValue(i,Value);
    171     }
    172   }
    173 }

  • trunk/source/processes/cuts/src/G4RToEConvForProton.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4RToEConvForProton.cc,v 1.3 2006/06/29 19:30:30 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4RToEConvForProton.cc,v 1.6 2009/09/14 07:27:46 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    5858}
    5959
     60
     61G4double G4RToEConvForProton::Convert(G4double rangeCut, const G4Material* )
     62{
     63  // Simple formula
     64  //   range = Ekin/(100*keV)*(1*mm);
     65  return (rangeCut/(1.0*mm)) * (100.0*keV);
     66}
     67
     68
     69// **********************************************************************
     70// ************************* ComputeLoss ********************************
     71// **********************************************************************
     72G4double G4RToEConvForProton::ComputeLoss(G4double AtomicNumber,
     73                                                G4double KineticEnergy) const
     74{
     75  //  calculate dE/dx
     76
     77  static G4double Z; 
     78  static G4double ionpot, tau0, taum, taul, ca, cba, cc;
     79
     80  G4double  z2Particle = theParticle->GetPDGCharge()/eplus;
     81  z2Particle *=  z2Particle;
     82  if (z2Particle < 0.1) return 0.0;
     83
     84  if( std::fabs(AtomicNumber-Z)>0.1 ){
     85    // recalculate constants
     86    Z = AtomicNumber;
     87    G4double Z13 = std::exp(std::log(Z)/3.);
     88    tau0 = 0.1*Z13*MeV/proton_mass_c2;
     89    taum = 0.035*Z13*MeV/proton_mass_c2;
     90    taul = 2.*MeV/proton_mass_c2;
     91    ionpot = 1.6e-5*MeV*std::exp(0.9*std::log(Z));
     92   cc = (taul+1.)*(taul+1.)*std::log(2.*electron_mass_c2*taul*(taul+2.)/ionpot)/(taul*(taul+2.))-1.;
     93    cc = 2.*twopi_mc2_rcl2*Z*cc*std::sqrt(taul);
     94    ca = cc/((1.-0.5*std::sqrt(tau0/taum))*tau0);
     95    cba = -0.5/std::sqrt(taum);
     96  }
     97
     98  G4double tau = KineticEnergy/theParticle->GetPDGMass();
     99  G4double dEdx;
     100  if ( tau <= tau0 ) {
     101    dEdx = ca*(std::sqrt(tau)+cba*tau);
     102  } else {
     103    if( tau <= taul ) {
     104      dEdx = cc/std::sqrt(tau);
     105    } else {
     106      dEdx = (tau+1.)*(tau+1.)*
     107             std::log(2.*electron_mass_c2*tau*(tau+2.)/ionpot)/(tau*(tau+2.))-1.;
     108      dEdx = 2.*twopi_mc2_rcl2*Z*dEdx;
     109    }
     110  }
     111  return dEdx*z2Particle ;
     112}
     113
     114
     115// **********************************************************************
     116// ************************* Reset       ********************************
     117// **********************************************************************
     118void G4RToEConvForProton::Reset()
     119{
     120  // do nothing because loss tables and range vectors are not used
     121  return;
     122}

  • trunk/source/processes/cuts/src/G4VRangeToEnergyConverter.cc

    r1055 r1196  
    2525//
    2626//
    27 // $Id: G4VRangeToEnergyConverter.cc,v 1.10 2009/04/02 02:43:42 kurasige Exp $
    28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
     27// $Id: G4VRangeToEnergyConverter.cc,v 1.15 2009/09/14 07:27:46 kurasige Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    3737#include "G4ParticleTable.hh"
    3838#include "G4Material.hh"
     39#include "G4MaterialTable.hh"
    3940#include "G4PhysicsLogVector.hh"
    4041
    … …  
    4546G4double  G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
    4647
     48// max energy cut
     49G4double  G4VRangeToEnergyConverter::MaxEnergyCut = 10.0*GeV;
     50
    4751G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
    48   theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(200),
     52  theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(300),
    4953  verboseLevel(1)
    5054{
    51 }
    52 
    53 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right)
     55  fMaxEnergyCut = 0.;
     56}
     57
     58G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right) : TotBin(right.TotBin)
    5459{
    5560  *this = right;
    … …  
    6671
    6772  NumberOfElements = right.NumberOfElements;
    68   TotBin = right.TotBin;
     73  //TotBin = right.TotBin;
    6974  theParticle = right.theParticle;
    7075  verboseLevel = right.verboseLevel;
    … …  
    7681  for (size_t j=0; j<size_t(NumberOfElements); j++){
    7782    G4LossVector* aVector= new
    78             G4LossVector(LowestEnergy, HighestEnergy, TotBin);
     83            G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
    7984    for (size_t i=0; i<size_t(TotBin); i++) {
    8085      G4double Value = (*((*right.theLossTable)[j]))[i];
    … …  
    8388    theLossTable->insert(aVector);
    8489  }
     90
     91  // clean up range vector store
     92  for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
     93    delete fRangeVectorStore.at(idx);
     94  }
     95  fRangeVectorStore.clear();
     96
     97  // copy range vector store
     98  for (size_t j=0; j<((right.fRangeVectorStore).size()); j++){
     99    G4RangeVector* vector = (right.fRangeVectorStore).at(j);
     100    G4RangeVector* rangeVector = 0;
     101    if (vector !=0 ) {
     102      rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
     103      for (size_t i=0; i<size_t(TotBin); i++) {
     104        G4double Value = (*vector)[i];
     105        rangeVector->PutValue(i,Value);
     106      }
     107    }
     108    fRangeVectorStore.push_back(rangeVector);
     109  }
    85110  return *this;
    86111}
    … …  
    89114G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
    90115{
    91   if (theLossTable) { 
    92     theLossTable->clearAndDestroy();
    93     delete theLossTable;
    94   }
    95   theLossTable=0;
     116  Reset();
    96117}
    97118
    … …  
    113134                                            const G4Material* material)
    114135{
    115   G4double Mass   = theParticle->GetPDGMass();
     136#ifdef G4VERBOSE
     137    if (GetVerboseLevel()>3) {
     138      G4cout << "G4VRangeToEnergyConverter::Convert() ";
     139      G4cout << "Convert for " << material->GetName()
     140             << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
     141    }
     142#endif
     143
    116144  G4double theKineticEnergyCuts = 0.;
     145
     146  if (fMaxEnergyCut != MaxEnergyCut) {
     147    fMaxEnergyCut = MaxEnergyCut;     
     148    // clear loss table and renge vectors
     149    Reset();
     150  }
    117151 
    118152  // Build the energy loss table
    … …  
    123157  G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
    124158
     159  // check density
     160  G4double density = material->GetDensity() ;
     161  if(density <= 0.) {
     162 #ifdef G4VERBOSE
     163    if (GetVerboseLevel()>0) {
     164      G4cout << "G4VRangeToEnergyConverter::Convert() ";
     165      G4cout << material->GetName() << "has zero density "
     166             << "( " << density << ")" << G4endl;
     167    }
     168#endif
     169    return 0.;
     170  }
     171 
     172   // initialize RangeVectorStore
     173  const G4MaterialTable* table = G4Material::GetMaterialTable();
     174  G4int ext_size = table->size() - fRangeVectorStore.size();
     175  for (int i=0; i<ext_size; i++) fRangeVectorStore.push_back(0);
     176 
     177  // Build Range Vector
    125178  G4int idx = material->GetIndex();
    126   G4double density = material->GetDensity() ;
    127   if(density > 0.) {
    128     G4RangeVector* rangeVector = new G4RangeVector(LowestEnergy, HighestEnergy, TotBin);
    129     BuildRangeVector(material, HighestEnergy, Mass, rangeVector);
    130     theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
    131 
    132     if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
    133            && (theKineticEnergyCuts < lowen) )
    134 
     179  G4RangeVector* rangeVector = fRangeVectorStore.at(idx);
     180  if (rangeVector == 0) {
     181    rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
     182    BuildRangeVector(material, rangeVector);
     183    fRangeVectorStore.at(idx) = rangeVector;
     184  }
     185
     186  // Convert Range Cut ro Kinetic Energy Cut
     187  theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
     188 
     189  if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
     190      && (theKineticEnergyCuts < lowen) ) {
    135191    //  corr. should be switched on smoothly   
    136     { theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
    137                                     tune/(rangeCut*density)); }
    138     if(theKineticEnergyCuts < LowestEnergy) {
    139       theKineticEnergyCuts = LowestEnergy ;
    140     }
    141     delete rangeVector;
    142   }
    143 
     192    theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
     193                             tune/(rangeCut*density));
     194  }
     195 
     196  if(theKineticEnergyCuts < LowestEnergy) {
     197    theKineticEnergyCuts = LowestEnergy ;
     198  } else if(theKineticEnergyCuts > MaxEnergyCut) {
     199    theKineticEnergyCuts = MaxEnergyCut;
     200  }
     201 
    144202  return theKineticEnergyCuts;
    145203}
    … …  
    175233
    176234// **********************************************************************
    177 // ************************ RangeLinSimpson *****************************
    178 // **********************************************************************
    179 G4double G4VRangeToEnergyConverter::RangeLinSimpson(
    180                                      G4int numberOfElement,
    181                                      const G4ElementVector* elementVector,
    182                                      const G4double* atomicNumDensityVector,
    183                                      G4double aMass,   
    184                                      G4double taulow, G4double tauhigh, G4int nbin)
    185 {
    186   // Simpson numerical integration, linear binning
    187   G4double dtau = (tauhigh-taulow)/nbin;
    188   G4double Value=0.;
    189   for (size_t i=0; i<=size_t(nbin); i++){
    190     G4double taui=taulow+dtau*i;
    191     G4double ti=aMass*taui;
    192     G4double lossi=0.;
    193     size_t nEl = (size_t)(numberOfElement);
    194     for (size_t j=0; j<nEl; j++) {
    195       G4bool isOut;
    196       G4int IndEl = (*elementVector)[j]->GetIndex();
    197       lossi += atomicNumDensityVector[j]*
    198               (*theLossTable)[IndEl]->GetValue(ti,isOut);
    199    }
    200     if ( i==0 ) {
    201       Value += 0.5/lossi;
    202     } else {
    203       if ( i<size_t(nbin) ) Value += 1./lossi;
    204       else            Value += 0.5/lossi;
    205     }
    206   }
    207   Value *= aMass*dtau;
    208 
    209   return Value;
    210 }
    211 
    212 
    213 // **********************************************************************
    214 // ************************ RangeLogSimpson *****************************
    215 // **********************************************************************
    216 G4double G4VRangeToEnergyConverter::RangeLogSimpson(
    217                                      G4int numberOfElement,
    218                                      const G4ElementVector* elementVector,
    219                                      const G4double* atomicNumDensityVector,
    220                                      G4double aMass,   
    221                                      G4double ltaulow, G4double ltauhigh,
    222                                      G4int nbin)
    223 {
    224   // Simpson numerical integration, logarithmic binning
    225   if(nbin<0) nbin = TotBin;
    226   G4double ltt = ltauhigh-ltaulow;
    227   G4double dltau = ltt/nbin;
    228   G4double Value = 0.;
    229   for (size_t i=0; i<=size_t(nbin); i++){
    230     G4double ui = ltaulow+dltau*i;
    231     G4double taui = std::exp(ui);
    232     G4double ti = aMass*taui;
    233     G4double lossi = 0.;
    234     size_t nEl = (size_t)(numberOfElement);
    235 
    236     for (size_t j=0; j<nEl; j++) {
    237       G4bool isOut;
    238       G4int IndEl = (*elementVector)[j]->GetIndex();
    239       lossi += atomicNumDensityVector[j]*
    240               (*theLossTable)[IndEl]->GetValue(ti,isOut);
    241     }
    242     if ( i==0 ) {
    243       Value +=  0.5*taui/lossi;
    244     } else {
    245       if ( i<size_t(nbin) ) Value += taui/lossi;
    246       else Value +=  0.5*taui/lossi;
    247     }
    248   }
    249   Value *= aMass*dltau;
    250 
    251   return Value;
    252 }
     235// ******************* Get/SetMaxEnergyCut  *****************************
     236// **********************************************************************
     237G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
     238{
     239  return MaxEnergyCut;
     240}
     241
     242void G4VRangeToEnergyConverter::SetMaxEnergyCut(G4double value)
     243{
     244  MaxEnergyCut = value;
     245}
     246
     247// **********************************************************************
     248// ************************ Reset  **************************************
     249// **********************************************************************
     250void G4VRangeToEnergyConverter::Reset()
     251{
     252  // delete loss table
     253  if (theLossTable) { 
     254    theLossTable->clearAndDestroy();
     255    delete theLossTable;
     256  }
     257  theLossTable=0;
     258  NumberOfElements=0;
     259 
     260  //clear RangeVectorStore
     261  for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
     262    delete fRangeVectorStore.at(idx);
     263  }
     264  fRangeVectorStore.clear();
     265}
     266
    253267
    254268// **********************************************************************
    … …  
    259273void G4VRangeToEnergyConverter::BuildLossTable()
    260274{
    261    //  Build dE/dx tables for elements
    262   if (size_t(NumberOfElements) != G4Element::GetNumberOfElements()) {
    263     if (theLossTable!=0) {
    264       theLossTable->clearAndDestroy();
    265       delete theLossTable;
    266     }
    267     theLossTable =0;
    268     NumberOfElements = 0;
    269   }
    270 
    271   if (NumberOfElements ==0) {
    272     NumberOfElements = G4Element::GetNumberOfElements();
    273     theLossTable = new G4LossTable();
    274     theLossTable->reserve(G4Element::GetNumberOfElements());
     275  if (size_t(NumberOfElements) == G4Element::GetNumberOfElements()) return;
     276 
     277  // clear Loss table and Range vectors
     278  Reset();
     279
     280  //  Build dE/dx tables for elements
     281  NumberOfElements = G4Element::GetNumberOfElements();
     282  theLossTable = new G4LossTable();
     283  theLossTable->reserve(G4Element::GetNumberOfElements());
    275284#ifdef G4VERBOSE
    276     if (GetVerboseLevel()>3) {
    277       G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
    278       G4cout << "Create theLossTable[" << theLossTable << "]";
    279       G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
    280     }
     285  if (GetVerboseLevel()>3) {
     286    G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
     287    G4cout << "Create theLossTable[" << theLossTable << "]";
     288    G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
     289  }
    281290#endif
    282  
    283 
    284     // fill the loss table
    285     for (size_t j=0; j<size_t(NumberOfElements); j++){
    286       G4double Value;
    287       G4LossVector* aVector= new
    288         G4LossVector(LowestEnergy, HighestEnergy, TotBin);
    289       for (size_t i=0; i<size_t(TotBin); i++) {
    290         Value = ComputeLoss(  (*G4Element::GetElementTable())[j]->GetZ(),
    291                               aVector->GetLowEdgeEnergy(i)
    292                               );
    293         aVector->PutValue(i,Value);
    294       }
    295       theLossTable->insert(aVector);
    296     }
    297   }
    298 }
    299 
    300 // **********************************************************************
    301 // ************************** ComputeLoss *******************************
    302 // **********************************************************************
    303 G4double G4VRangeToEnergyConverter::ComputeLoss(G4double AtomicNumber,
    304                                                 G4double KineticEnergy) const
    305 {
    306   //  calculate dE/dx
    307 
    308   static G4double Z; 
    309   static G4double ionpot, tau0, taum, taul, ca, cba, cc;
    310 
    311   G4double  z2Particle = theParticle->GetPDGCharge()/eplus;
    312   z2Particle *=  z2Particle;
    313   if (z2Particle < 0.1) return 0.0;
    314 
    315   if( std::fabs(AtomicNumber-Z)>0.1 ){
    316     // recalculate constants
    317     Z = AtomicNumber;
    318     G4double Z13 = std::exp(std::log(Z)/3.);
    319     tau0 = 0.1*Z13*MeV/proton_mass_c2;
    320     taum = 0.035*Z13*MeV/proton_mass_c2;
    321     taul = 2.*MeV/proton_mass_c2;
    322     ionpot = 1.6e-5*MeV*std::exp(0.9*std::log(Z));
    323     cc = (taul+1.)*(taul+1.)*std::log(2.*electron_mass_c2*taul*(taul+2.)/ionpot)/(taul*(taul+2.))-1.;
    324     cc = 2.*twopi_mc2_rcl2*Z*cc*std::sqrt(taul);
    325     ca = cc/((1.-0.5*std::sqrt(tau0/taum))*tau0);
    326     cba = -0.5/std::sqrt(taum);
    327   }
    328 
    329   G4double tau = KineticEnergy/theParticle->GetPDGMass();
    330   G4double dEdx;
    331   if ( tau <= tau0 ) {
    332     dEdx = ca*(std::sqrt(tau)+cba*tau);
    333   } else {
    334     if( tau <= taul ) {
    335       dEdx = cc/std::sqrt(tau);
    336     } else {
    337       dEdx = (tau+1.)*(tau+1.)*
    338              std::log(2.*electron_mass_c2*tau*(tau+2.)/ionpot)/(tau*(tau+2.))-1.;
    339       dEdx = 2.*twopi_mc2_rcl2*Z*dEdx;
    340     }
    341   }
    342   return dEdx*z2Particle ;
     291 
     292 
     293  // fill the loss table
     294  for (size_t j=0; j<size_t(NumberOfElements); j++){
     295    G4double Value;
     296    G4LossVector* aVector= 0;
     297    aVector= new G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
     298    for (size_t i=0; i<size_t(TotBin); i++) {
     299      Value = ComputeLoss(  (*G4Element::GetElementTable())[j]->GetZ(),
     300                            aVector->GetLowEdgeEnergy(i)
     301                            );
     302      aVector->PutValue(i,Value);
     303    }
     304    theLossTable->insert(aVector);
     305  }
    343306}
    344307
    … …  
    346309// ************************ BuildRangeVector ****************************
    347310// **********************************************************************
    348 void G4VRangeToEnergyConverter::BuildRangeVector(
    349                                   const G4Material* aMaterial,
    350                                   G4double       maxEnergy,
    351                                   G4double       aMass,
    352                                   G4RangeVector* rangeVector)
     311void G4VRangeToEnergyConverter::BuildRangeVector(const G4Material* aMaterial,
     312                                             G4PhysicsLogVector* rangeVector)
    353313{
    354314  //  create range vector for a material
    355   const G4double tlim=2.*MeV, t1=0.1*MeV, t2=0.025*MeV;
    356   const G4int  maxnbint=100;
    357  
    358315  const G4ElementVector* elementVector = aMaterial->GetElementVector();
    359316  const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
    360 
    361317  G4int NumEl = aMaterial->GetNumberOfElements();
    362318
    363319  // calculate parameters of the low energy part first
    364   G4double loss1=0.;
    365   G4double loss2=0.;
    366320  size_t i;
    367   for (i=0; i<size_t(NumEl); i++) {
    368     G4bool isOut;
    369     G4int IndEl = (*elementVector)[i]->GetIndex();
    370     loss1 += atomicNumDensityVector[i]*
    371             (*theLossTable)[IndEl]->GetValue(t1,isOut);
    372     loss2 += atomicNumDensityVector[i]*
    373             (*theLossTable)[IndEl]->GetValue(t2,isOut);
    374   }
    375   G4double tau1 = t1/proton_mass_c2;
    376   G4double sqtau1 = std::sqrt(tau1);
    377   G4double ca = (4.*loss2-loss1)/sqtau1;
    378   G4double cb = (2.*loss1-4.*loss2)/tau1;
    379   G4double cba = cb/ca;
    380   G4double taulim = tlim/proton_mass_c2;
    381   G4double taumax = maxEnergy/aMass;
    382   G4double ltaumax = std::log(taumax);
    383 
    384   // now we can fill the range vector....
    385   G4double  rmax = 0.0;
    386   for (i=0; i<size_t(TotBin); i++) {
    387     G4double  LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
    388     G4double  tau = LowEdgeEnergy/aMass;
    389     G4double  Value;
    390  
    391     if ( tau <= tau1 ){
    392       Value =2.*aMass*std::log(1.+cba*std::sqrt(tau))/cb;
     321  std::vector<G4double> lossV;
     322  for ( size_t ib=0; ib<size_t(TotBin); ib++) {
     323    G4double loss=0.;
     324    for (i=0; i<size_t(NumEl); i++) {
     325      G4int IndEl = (*elementVector)[i]->GetIndex();
     326      loss += atomicNumDensityVector[i]*
     327                (*((*theLossTable)[IndEl]))[ib];
     328    }
     329    lossV.push_back(loss);
     330  }
     331   
     332  // Integrate with Simpson formula with logarithmic binning
     333  G4double ltt = std::log(MaxEnergyCut/LowestEnergy);
     334  G4double dltau = ltt/TotBin;
     335
     336  G4double s0 = 0.;
     337  G4double Value;
     338  for ( i=0; i<size_t(TotBin); i++) {
     339    G4double t = rangeVector->GetLowEdgeEnergy(i);
     340    G4double s = t/lossV[i];
     341    if (i==0) s0 += 0.5*s;
     342    else s0 += s;
     343   
     344    if (i==0) {
     345       Value = (s0 + 0.5*s)*dltau ;
    393346    } else {
    394       Value = 2.*aMass*std::log(1.+cba*sqtau1)/cb;
    395       if ( tau <= taulim ) {
    396         G4int nbin = (G4int)(maxnbint*(tau-tau1)/(taulim-tau1));
    397         if ( nbin<1 ) nbin = 1;
    398         Value += RangeLinSimpson( NumEl, elementVector,
    399                                   atomicNumDensityVector, aMass,
    400                                   tau1, tau, nbin);
    401       } else {
    402         Value += RangeLinSimpson( NumEl, elementVector,
    403                                   atomicNumDensityVector, aMass,
    404                                   tau1, taulim, maxnbint);
    405         G4double ltaulow  = std::log(taulim);
    406         G4double ltauhigh = std::log(tau);
    407         G4int nbin = (G4int)(maxnbint*(ltauhigh-ltaulow)/(ltaumax-ltaulow));
    408         if ( nbin<1 ) nbin = 1;
    409         Value += RangeLogSimpson(NumEl, elementVector,
    410                                  atomicNumDensityVector, aMass,
    411                                  ltaulow, ltauhigh, nbin);
    412       }
    413     }
    414     rangeVector->PutValue(i,Value);
    415     if (rmax < Value) rmax = Value;
    416   }
    417 }
     347      Value = (s0 - 0.5*s)*dltau ;
     348    }
     349    rangeVector->PutValue(i,Value);
     350  }
     351}
    418352
    419353// **********************************************************************
    … …  
    430364  //  find max. range and the corresponding energy (rmax,Tmax)
    431365  G4double rmax= -1.e10*mm;
    432   G4double Tmax= HighestEnergy;
    433   G4double fac = std::exp( std::log(HighestEnergy/LowestEnergy)/TotBin );
    434   G4double T=LowestEnergy/fac;
    435   G4bool isOut;
    436 
     366
     367  G4double T1 = LowestEnergy;
     368  G4double r1 =(*rangeVector)[0] ;
     369
     370  G4double T2 = MaxEnergyCut;
     371
     372  // check theCutInLength < r1
     373  if ( theCutInLength <= r1 ) {  return T1; }
     374
     375  // scan range vector to find nearest bin
     376  // ( suppose that r(Ti) > r(Tj) if Ti >Tj )
    437377  for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
    438     T *= fac;
    439     G4double r=rangeVector->GetValue(T,isOut);
    440     if ( r>rmax )    {
    441        Tmax=T;
    442        rmax=r;
     378    G4double T=rangeVector->GetLowEdgeEnergy(ibin);
     379    G4double r=(*rangeVector)[ibin];
     380    if ( r>rmax )   rmax=r;
     381    if (r <theCutInLength ) {
     382      T1 = T;
     383      r1 = r;
     384    } else if (r >theCutInLength ) {
     385      T2 = T;
     386      break;
    443387    }
    444388  }
    … …  
    453397      G4cout << " is too big  " ;
    454398      G4cout << " for material  idx=" << materialIndex <<G4endl;
    455       G4cout << "The cut in energy is set" << DBL_MAX/GeV << "GeV " <<G4endl;
    456399    }
    457400#endif
    458     return  DBL_MAX;
     401    return  MaxEnergyCut;
    459402  }
    460403 
    461404  // convert range to energy
    462   G4double T1 = LowestEnergy;
    463   G4double r1 = rangeVector->GetValue(T1,isOut);
    464   if ( theCutInLength <= r1 )
    465   {
    466     return T1;
    467   }
    468 
    469   G4double T2 = Tmax ;
    470405  G4double T3 = std::sqrt(T1*T2);
    471   G4double r3 = rangeVector->GetValue(T3,isOut);
     406  G4double r3 = rangeVector->Value(T3);
    472407  while ( std::fabs(1.-r3/theCutInLength)>epsilon ) {
    473408    if ( theCutInLength <= r3 ) {
    … …  
    477412    }
    478413    T3 = std::sqrt(T1*T2);
    479     r3 = rangeVector->GetValue(T3,isOut);
     414    r3 = rangeVector->Value(T3);
    480415  }
    481416
Note: See TracChangeset for help on using the changeset viewer.