Ignore:
Timestamp:
Apr 20, 2009, 4:53:50 PM (17 years ago)
Author:
garnier
Message:

fichier oublies

Location:
trunk/source/processes/electromagnetic/utils/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.cc,v 1.44 2008/08/03 18:47:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmCalculator.cc,v 1.46 2009/02/24 09:56:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    755755    if(currentProcess) currentProcessName = currentProcess->GetProcessName();
    756756
    757     if(p->GetParticleType() == "nucleus" &&
    758        currentParticleName  != "deuteron" && currentParticleName != "triton") {
     757    if(p->GetParticleType() == "nucleus"
     758       && currentParticleName != "deuteron" 
     759       && currentParticleName != "triton"
     760       && currentParticleName != "alpha+"
     761       && currentParticleName != "helium"
     762       && currentParticleName != "hydrogen"
     763      ) {
    759764      baseParticle = theGenericIon;
    760765      massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();

  • trunk/source/processes/electromagnetic/utils/src/G4EmConfigurator.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.cc,v 1.3 2008/11/21 12:30:29 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmConfigurator.cc,v 1.4 2009/02/26 11:33:33 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    117117  G4String fname = "";
    118118  if(fm) fname = fm->GetName();
    119   AddModelForRegion(particleName, processName, mod->GetName(), regionName,
     119  G4String mname = "";
     120  if(mod) mname = mod->GetName();
     121  AddModelForRegion(particleName, processName, mname, regionName,
    120122                    emin, emax, fname);
    121123}
    … …  
    203205
    204206        for(G4int i=0; i<nm; i++) {
    205           if(modelName == modelList[i]->GetName() &&
     207          G4String mname = "";
     208          if(modelList[i]) mname = modelList[i]->GetName();
     209          G4String fname = "";
     210          if(flucModelList[i]) fname = flucModelList[i]->GetName();
     211          if(modelName == mname && flucModelName == fname &&
    206212             (particleList[i] == "" || particleList[i] == particleName) ) {
    207213            mod  = modelList[i];
    … …  
    214220
    215221        if(!mod) {
    216           G4cout << "### G4EmConfigurator WARNING: fails to find a model <"
    217                  << modelName << "> for process <"
    218                  << processName << "> and " << particleName
    219                  << G4endl;
    220           if(flucModelName != "") 
    221             G4cout << "                            fluctuation model <"
    222                    << flucModelName << G4endl;
     222
     223          // set fluctuation model for ionisation processes
     224          if(fluc && ptype == eloss) {
     225            G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc);
     226            p->SetFluctModel(fluc);
     227         
     228          } else {
     229            G4cout << "### G4EmConfigurator WARNING: fails to find a model <"
     230                   << modelName << "> for process <"
     231                   << processName << "> and " << particleName
     232                   << G4endl;
     233            if(flucModelName != "") {
     234              G4cout << "                            fluctuation model <"
     235                     << flucModelName << G4endl;
     236            }
     237          }
    223238        } else {
    224239

  • trunk/source/processes/electromagnetic/utils/src/G4EmProcessOptions.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmProcessOptions.cc,v 1.24 2008/04/17 10:33:27 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmProcessOptions.cc,v 1.26 2009/02/18 14:43:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5959#include "G4VMultipleScattering.hh"
    6060#include "G4Region.hh"
     61#include "G4RegionStore.hh"
    6162
    6263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    392393//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    393394
    394 void G4EmProcessOptions::ActivateDeexcitation(G4bool val, const G4Region* r)
    395 {
    396   const std::vector<G4VEnergyLossProcess*>& v =
    397         theManager->GetEnergyLossProcessVector();
    398   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    399   for(itr = v.begin(); itr != v.end(); itr++) {
    400     G4VEnergyLossProcess* p = *itr;
    401     if(p) p->ActivateDeexcitation(val,r);
    402   }
    403   const std::vector<G4VEmProcess*>& w =
    404         theManager->GetEmProcessVector();
    405   std::vector<G4VEmProcess*>::const_iterator itp;
    406   for(itp = w.begin(); itp != w.end(); itp++) {
    407     G4VEmProcess* q = *itp;
    408     if(q) q->ActivateDeexcitation(val,r);
     395void G4EmProcessOptions::ActivateDeexcitation(const G4String& pname,
     396                                              G4bool val,
     397                                              const G4String& reg)
     398{
     399  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     400  const G4Region* r = 0;
     401  if(reg == "" || reg == "World") {
     402    r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
     403  } else {
     404    r = regionStore->GetRegion(reg, false);
     405  }
     406  if(!r) {
     407    G4cout << "G4EmProcessOptions::ActivateDeexcitation ERROR: G4Region <"
     408           << reg << "> not found, the command ignored" << G4endl;
     409    return;
     410  }
     411
     412  const std::vector<G4VEnergyLossProcess*>& v =
     413        theManager->GetEnergyLossProcessVector();
     414  std::vector<G4VEnergyLossProcess*>::const_iterator itr;
     415  for(itr = v.begin(); itr != v.end(); itr++) {
     416    G4VEnergyLossProcess* p = *itr;
     417    if(p) {
     418      if(pname == p->GetProcessName()) p->ActivateDeexcitation(val,r);
     419    }
     420  }
     421  const std::vector<G4VEmProcess*>& w =
     422        theManager->GetEmProcessVector();
     423  std::vector<G4VEmProcess*>::const_iterator itp;
     424  for(itp = w.begin(); itp != w.end(); itp++) {
     425    G4VEmProcess* q = *itp;
     426    if(q) {
     427      if(pname == q->GetProcessName()) q->ActivateDeexcitation(val,r);
     428    }
    409429  }
    410430}

  • trunk/source/processes/electromagnetic/utils/src/G4EnergyLossMessenger.cc

    r991 r1005  
    2525//
    2626//
    27 // $Id: G4EnergyLossMessenger.cc,v 1.35 2008/10/20 13:27:45 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4EnergyLossMessenger.cc,v 1.37 2009/02/18 14:43:27 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    178178  aplCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
    179179
     180  deexCmd = new G4UIcommand("/process/em/deexcitation",this);
     181  deexCmd->SetGuidance("Set deexcitation flag per process and G4Region.");
     182  deexCmd->SetGuidance("  procName  : process name");
     183  deexCmd->SetGuidance("  flag      : flag");
     184  deexCmd->SetGuidance("  regName   : G4Region name");
     185
     186  G4UIparameter* pName = new G4UIparameter("pName",'s',false);
     187  deexCmd->SetParameter(pName);
     188
     189  G4UIparameter* flag = new G4UIparameter("flag",'s',false);
     190  deexCmd->SetParameter(flag);
     191
     192  G4UIparameter* regName = new G4UIparameter("regName",'s',false);
     193  deexCmd->SetParameter(regName);
     194
    180195  dedxCmd = new G4UIcmdWithAnInteger("/process/eLoss/binsDEDX",this);
    181196  dedxCmd->SetGuidance("Set number of bins for DEDX tables");
    … …  
    259274  delete MinSubSecCmd;
    260275  delete StepFuncCmd;
     276  delete deexCmd;
    261277  delete eLossDirectory;
    262278  delete mscDirectory;
    … …  
    320336  }
    321337
     338  if (command == deexCmd) {
     339    G4String s1 (""), s2(""), s3("");
     340    G4bool b = false;
     341    std::istringstream is(newValue);
     342    is >> s1 >> s2 >> s3;
     343    if(s2 == "true") b = true;
     344    opt->ActivateDeexcitation(s1,b,s3);
     345  }
     346
    322347  if (command == mscCmd) {
    323348    if(newValue == "Minimal")

  • trunk/source/processes/electromagnetic/utils/src/G4LossTableBuilder.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableBuilder.cc,v 1.27 2008/07/22 15:55:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LossTableBuilder.cc,v 1.28 2009/02/18 16:24:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    132132      G4double dedx1  = pv->GetValue(elow, b);
    133133
     134      //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;
     135
    134136      // protection for specific cases dedx=0
    135137      if(dedx1 == 0.0) {
    136         for (size_t k=1; k<nbins; k++) {
     138        for (size_t k=1; k<nbins-1; k++) {
    137139          bin0++;
    138140          elow  = pv->GetLowEdgeEnergy(k);
    … …  
    143145      }
    144146
     147      //G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh << G4endl;
    145148      // initialisation of a new vector
    146149      G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
     150      // dedx is exect zero
     151      if(nbins <= 2) {
     152        v->PutValue(0,1000.);
     153        v->PutValue(1,2000.);
     154        G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
     155        return;
     156      }
    147157      v->SetSpline(splineFlag);
    148158

  • trunk/source/processes/electromagnetic/utils/src/G4VEmFluctuationModel.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmFluctuationModel.cc,v 1.3 2008/07/15 16:56:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmFluctuationModel.cc,v 1.4 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6565}
    6666
     67void G4VEmFluctuationModel::InitialiseMe(const G4ParticleDefinition*)
     68{}
    6769
     70void G4VEmFluctuationModel::SetParticleAndCharge(const G4ParticleDefinition*,
     71                                                 G4double)
     72{}
     73
     74

  • trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.cc,v 1.20 2008/11/13 23:13:18 schaelic Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmModel.cc,v 1.25 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4141// 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
    4242// 06.02.2006 add method ComputeMeanFreePath() (mma)
     43// 16.02.2009 Move implementations of virtual methods to source (VI)
    4344//
    4445//
    … …  
    6061  fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
    6162  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
    62   pParticleChange(0),nuclearStopping(false),nsec(5)
     63  pParticleChange(0),nuclearStopping(false),
     64  currentCouple(0),currentElement(0),
     65  nsec(5),flagDeexcitation(false)
    6366{
    6467  xsec.resize(nsec);
    … …  
    8285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    8386
    84 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
    85                                            const G4ParticleDefinition* p,
    86                                            G4double ekin,
    87                                            G4double emin,
    88                                            G4double emax)
    89 {
    90   SetupForMaterial(p, material, ekin);
    91   G4double cross = 0.0;
    92   const G4ElementVector* theElementVector = material->GetElementVector();
    93   const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
    94   G4int nelm = material->GetNumberOfElements();
    95   if(nelm > nsec) {
    96     xsec.resize(nelm);
    97     nsec = nelm;
    98   }
    99   for (G4int i=0; i<nelm; i++) {
    100     cross += theAtomNumDensityVector[i]*
    101       ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
    102     xsec[i] = cross;
    103   }
    104   return cross;
    105 }
    106 
    107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    108 
    109 G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
    110                                          G4double ekin,
    111                                          const G4Material* material,     
    112                                          G4double emin,
    113                                          G4double emax)
    114 {
    115   G4double mfp = DBL_MAX;
    116   G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
    117   if (cross > DBL_MIN) mfp = 1./cross;
    118   return mfp;
    119 }
    120 
    121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    122 
    12387void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
    12488                                            const G4DataVector& cuts)
    12589{
     90  // initialise before run
     91  flagDeexcitation = false;
     92
    12693  G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
    12794  if(nbins < 3) nbins = 3;
    … …  
    162129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    163130
    164 
     131G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
     132                                          const G4ParticleDefinition*,
     133                                          G4double,G4double)
     134{
     135  return 0.0;
     136}
     137
     138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     139
     140G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
     141                                           const G4ParticleDefinition* p,
     142                                           G4double ekin,
     143                                           G4double emin,
     144                                           G4double emax)
     145{
     146  SetupForMaterial(p, material, ekin);
     147  G4double cross = 0.0;
     148  const G4ElementVector* theElementVector = material->GetElementVector();
     149  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
     150  G4int nelm = material->GetNumberOfElements();
     151  if(nelm > nsec) {
     152    xsec.resize(nelm);
     153    nsec = nelm;
     154  }
     155  for (G4int i=0; i<nelm; i++) {
     156    cross += theAtomNumDensityVector[i]*
     157      ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
     158    xsec[i] = cross;
     159  }
     160  return cross;
     161}
     162
     163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     164
     165G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     166                                                G4double, G4double, G4double,
     167                                                G4double, G4double)
     168{
     169  return 0.0;
     170}
     171
     172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     173
     174G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
     175                                  const G4MaterialCutsCouple*)
     176{
     177  return 0.0;
     178}
     179
     180//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     181
     182G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
     183                                          const G4Material*, G4double)
     184{
     185  G4double q = p->GetPDGCharge()/CLHEP::eplus;
     186  return q*q;
     187}
     188
     189//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     190
     191G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
     192                                       const G4Material*, G4double)
     193{
     194  return p->GetPDGCharge();
     195}
     196
     197//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     198
     199void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
     200                                      const G4DynamicParticle*,
     201                                      G4double&,G4double&,G4double)
     202{}
     203
     204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     205
     206void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,
     207                                             const G4Track&,
     208                                             G4double& )
     209{}
     210
     211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     212
     213G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     214                                        G4double kineticEnergy)
     215{
     216  return kineticEnergy;
     217}
     218
     219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     220
     221void G4VEmModel::SampleScattering(const G4DynamicParticle*, G4double)
     222{}
     223
     224//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     225
     226G4double G4VEmModel::ComputeTruePathLengthLimit(const G4Track&,
     227                                                G4PhysicsTable*,
     228                                                G4double)
     229{
     230  return DBL_MAX;
     231}
     232
     233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     234
     235G4double G4VEmModel::ComputeGeomPathLength(G4double truePathLength)
     236{
     237  return truePathLength;
     238}
     239
     240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     241
     242G4double G4VEmModel::ComputeTrueStepLength(G4double geomPathLength)
     243{
     244  return geomPathLength;
     245}
     246
     247//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     248
     249void G4VEmModel::DefineForRegion(const G4Region*)
     250{}
     251
     252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     253
     254void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
     255                                  const G4Material*, G4double)
     256{}
     257
     258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.cc,v 1.60 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmProcess.cc,v 1.62 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9191  applyCuts(false),
    9292  startFromNull(true),
    93   nRegions(0),
    94   selectedModel(0),
     93  useDeexcitation(false),
     94  nDERegions(0),
     95  idxDERegions(0),
     96  currentModel(0),
    9597  particle(0),
    9698  currentCouple(0)
    … …  
    135137  delete modelManager;
    136138  (G4LossTableManager::Instance())->DeRegister(this);
     139}
     140
     141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     142
     143void G4VEmProcess::Clear()
     144{
     145  delete [] theEnergyOfCrossSectionMax;
     146  delete [] theCrossSectionMax;
     147  delete [] idxDERegions;
     148  theEnergyOfCrossSectionMax = 0;
     149  theCrossSectionMax = 0;
     150  idxDERegions = 0;
     151  currentCouple = 0;
     152  preStepLambda = 0.0;
     153  mfpKinEnergy  = DBL_MAX;
     154  deRegions.clear();
     155  nDERegions = 0;
    137156}
    138157
    … …  
    162181      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
    163182  }
    164 }
    165 
    166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    167 
    168 void G4VEmProcess::Clear()
    169 {
    170   if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
    171   if(theCrossSectionMax) delete [] theCrossSectionMax;
    172   theEnergyOfCrossSectionMax = 0;
    173   theCrossSectionMax = 0;
    174   currentCouple = 0;
    175   preStepLambda = 0.0;
    176   mfpKinEnergy  = DBL_MAX;
     183  // Sub Cutoff and Deexcitation
     184  if (nDERegions>0) {
     185
     186    const G4ProductionCutsTable* theCoupleTable=
     187          G4ProductionCutsTable::GetProductionCutsTable();
     188    size_t numOfCouples = theCoupleTable->GetTableSize();
     189
     190    idxDERegions = new G4bool[numOfCouples];
     191
     192    for (size_t j=0; j<numOfCouples; j++) {
     193
     194      const G4MaterialCutsCouple* couple =
     195        theCoupleTable->GetMaterialCutsCouple(j);
     196      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
     197      G4bool reg = false;
     198      for(G4int i=0; i<nDERegions; i++) {
     199        if(deRegions[i]) {
     200          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     201        }
     202      }
     203      idxDERegions[j] = reg;
     204    }
     205  }
     206  if (1 < verboseLevel && nDERegions>0) {
     207    G4cout << " Deexcitation is activated for regions: " << G4endl;
     208    for (G4int i=0; i<nDERegions; i++) {
     209      const G4Region* r = deRegions[i];
     210      G4cout << "           " << r->GetName() << G4endl;
     211    }
     212  }
    177213}
    178214
    … …  
    237273      G4cout << *theLambdaTable << G4endl;
    238274    }
     275  }
     276}
     277
     278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     279
     280void G4VEmProcess::PrintInfoDefinition()
     281{
     282  if(verboseLevel > 0) {
     283    G4cout << G4endl << GetProcessName() << ":   for  "
     284           << particle->GetParticleName();
     285    if(integral) G4cout << ", integral: 1 ";
     286    if(applyCuts) G4cout << ", applyCuts: 1 ";
     287    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
     288    if(buildLambdaTable) {
     289      G4cout << "      Lambda tables from "
     290             << G4BestUnit(minKinEnergy,"Energy")
     291             << " to "
     292             << G4BestUnit(maxKinEnergy,"Energy")
     293             << " in " << nLambdaBins << " bins, spline: "
     294             << (G4LossTableManager::Instance())->SplineFlag()
     295             << G4endl;
     296    }
     297    PrintInfo();
     298    modelManager->DumpModelList(verboseLevel);
     299  }
     300
     301  if(verboseLevel > 2 && buildLambdaTable) {
     302    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     303    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
    239304  }
    240305}
    … …  
    304369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    305370
    306 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
    307                                        G4double,
    308                                        G4ForceCondition* condition)
    309 {
    310   *condition = NotForced;
    311   return G4VEmProcess::MeanFreePath(track);
    312 }
    313 
    314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    315 
    316371G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
    317372                                              const G4Step&)
    … …  
    342397  }
    343398
    344   G4VEmModel* currentModel = SelectModel(finalT);
    345 
     399  SelectModel(finalT);
     400  if(useDeexcitation) {
     401    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     402  }
    346403  /* 
    347404  if(0 < verboseLevel) {
    … …  
    404461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    405462
    406 void G4VEmProcess::PrintInfoDefinition()
    407 {
    408   if(verboseLevel > 0) {
    409     G4cout << G4endl << GetProcessName() << ":   for  "
    410            << particle->GetParticleName();
    411     if(integral) G4cout << ", integral: 1 ";
    412     if(applyCuts) G4cout << ", applyCuts: 1 ";
    413     G4cout << "    SubType= " << GetProcessSubType() << G4endl;
    414     if(buildLambdaTable) {
    415       G4cout << "      Lambda tables from "
    416              << G4BestUnit(minKinEnergy,"Energy")
    417              << " to "
    418              << G4BestUnit(maxKinEnergy,"Energy")
    419              << " in " << nLambdaBins << " bins, spline: "
    420              << (G4LossTableManager::Instance())->SplineFlag()
    421              << G4endl;
    422     }
    423     PrintInfo();
    424     modelManager->DumpModelList(verboseLevel);
    425   }
    426 
    427   if(verboseLevel > 2 && buildLambdaTable) {
    428     G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
    429     if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
    430   }
    431 }
    432 
    433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    434 
    435 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
    436                                              const G4MaterialCutsCouple* couple)
    437 {
    438   // Cross section per atom is calculated
    439   DefineMaterial(couple);
    440   G4double cross = 0.0;
    441   G4bool b;
    442   if(theLambdaTable) {
    443     cross = (((*theLambdaTable)[currentMaterialIndex])->
    444                            GetValue(kineticEnergy, b));
    445   } else {
    446     G4VEmModel* model = SelectModel(kineticEnergy);
    447     cross =
    448       model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
    449   }
    450 
    451   return cross;
    452 }
    453 
    454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    455 
    456463G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
    457464                                       const G4String& directory,
    … …  
    521528
    522529  return yes;
     530}
     531
     532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     533
     534void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     535{
     536  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     537  const G4Region* reg = r;
     538  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     539
     540  // the region is in the list
     541  if (nDERegions) {
     542    for (G4int i=0; i<nDERegions; i++) {
     543      if (reg == deRegions[i]) {
     544        if(!val) deRegions[i] = 0;
     545        return;
     546      }
     547    }
     548  }
     549
     550  // new region
     551  if(val) {
     552    useDeexcitation = true;
     553    deRegions.push_back(reg);
     554    nDERegions++;
     555  } else {
     556    useDeexcitation = false;
     557  }
     558}
     559
     560//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     561
     562G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
     563                                             const G4MaterialCutsCouple* couple)
     564{
     565  // Cross section per atom is calculated
     566  DefineMaterial(couple);
     567  G4double cross = 0.0;
     568  G4bool b;
     569  if(theLambdaTable) {
     570    cross = (((*theLambdaTable)[currentMaterialIndex])->
     571                           GetValue(kineticEnergy, b));
     572  } else {
     573    SelectModel(kineticEnergy);
     574    cross = currentModel->CrossSectionPerVolume(currentMaterial,
     575                                                particle,kineticEnergy);
     576  }
     577
     578  return cross;
     579}
     580
     581//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     582
     583G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
     584                                       G4double,
     585                                       G4ForceCondition* condition)
     586{
     587  *condition = NotForced;
     588  return G4VEmProcess::MeanFreePath(track);
    523589}
    524590

  • trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.143 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.146 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    150150  secondaryParticle(0),
    151151  nSCoffRegions(0),
     152  nDERegions(0),
    152153  idxSCoffRegions(0),
     154  idxDERegions(0),
    153155  nProcesses(0),
    154156  theDEDXTable(0),
    … …  
    176178  isIonisation(true),
    177179  useSubCutoff(false),
     180  useDeexcitation(false),
    178181  particle(0),
    179182  currentCouple(0),
    … …  
    237240           << G4endl;
    238241  delete vstrag;
    239   Clear();
     242  Clean();
    240243
    241244  if ( !baseParticle ) {
    … …  
    291294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    292295
    293 void G4VEnergyLossProcess::Clear()
     296void G4VEnergyLossProcess::Clean()
    294297{
    295298  if(1 < verboseLevel) {
    … …  
    302305  delete [] theCrossSectionMax;
    303306  delete [] idxSCoffRegions;
     307  delete [] idxDERegions;
    304308
    305309  theDEDXAtMaxEnergy = 0;
    … …  
    312316  scProcesses.clear();
    313317  nProcesses = 0;
     318}
     319
     320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     321
     322G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
     323                                                const G4Material*,
     324                                                G4double cut)
     325{
     326  return cut;
    314327}
    315328
    … …  
    364377  }
    365378
    366   Clear();
     379  Clean();
    367380
    368381  // Base particle and set of models can be defined here
    … …  
    407420                                     minSubRange, verboseLevel);
    408421
    409   // Sub Cutoff Regime
    410   if (nSCoffRegions>0) {
     422  // Sub Cutoff and Deexcitation
     423  if (nSCoffRegions>0 || nDERegions>0) {
    411424    theSubCuts = modelManager->SubCutoff();
    412425
    … …  
    414427          G4ProductionCutsTable::GetProductionCutsTable();
    415428    size_t numOfCouples = theCoupleTable->GetTableSize();
    416     idxSCoffRegions = new G4int[numOfCouples];
     429
     430    if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
     431    if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
    417432 
    418433    for (size_t j=0; j<numOfCouples; j++) {
    … …  
    421436        theCoupleTable->GetMaterialCutsCouple(j);
    422437      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    423       G4int reg = 0;
    424       for(G4int i=0; i<nSCoffRegions; i++) {
    425         if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
    426       }
    427       idxSCoffRegions[j] = reg;
     438     
     439      if(nSCoffRegions>0) {
     440        G4bool reg = false;
     441        for(G4int i=0; i<nSCoffRegions; i++) {
     442          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
     443        }
     444        idxSCoffRegions[j] = reg;
     445      }
     446      if(nDERegions>0) {
     447        G4bool reg = false;
     448        for(G4int i=0; i<nDERegions; i++) {
     449          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     450        }
     451        idxDERegions[j] = reg;
     452      }
    428453    }
    429454  }
    … …  
    445470      }
    446471    }
     472    if (nDERegions) {
     473      G4cout << " Deexcitation is ON for regions: " << G4endl;
     474      for (G4int i=0; i<nDERegions; i++) {
     475        const G4Region* r = deRegions[i];
     476        G4cout << "           " << r->GetName() << G4endl;
     477      }
     478    }
    447479  }
    448480}
    … …  
    479511    if(isIonisation) G4cout << "  isIonisation  flag = 1";
    480512    G4cout << G4endl;
    481   }
    482 }
    483 
    484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    485 
    486 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
    487 {
    488   G4RegionStore* regionStore = G4RegionStore::GetInstance();
    489   if(val) {
    490     useSubCutoff = true;
    491     if (!r) {r = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
    492     if (nSCoffRegions) {
    493       for (G4int i=0; i<nSCoffRegions; i++) {
    494         if (r == scoffRegions[i]) return;
    495       }
    496     }
    497     scoffRegions.push_back(r);
    498     nSCoffRegions++;
    499   } else {
    500     useSubCutoff = false;
    501513  }
    502514}
    … …  
    640652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    641653
    642 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
    643                 const G4Track&,
    644                 G4double, G4double, G4double&)
    645 {
    646   return DBL_MAX;
     654void G4VEnergyLossProcess::PrintInfoDefinition()
     655{
     656  if(0 < verboseLevel) {
     657    G4cout << G4endl << GetProcessName() << ":   for  "
     658           << particle->GetParticleName()
     659           << "    SubType= " << GetProcessSubType()
     660           << G4endl
     661           << "      dE/dx and range tables from "
     662           << G4BestUnit(minKinEnergy,"Energy")
     663           << " to " << G4BestUnit(maxKinEnergy,"Energy")
     664           << " in " << nBins << " bins" << G4endl
     665           << "      Lambda tables from threshold to "
     666           << G4BestUnit(maxKinEnergy,"Energy")
     667           << " in " << nBins << " bins, spline: "
     668           << (G4LossTableManager::Instance())->SplineFlag()
     669           << G4endl;
     670    if(theRangeTableForLoss && isIonisation) {
     671      G4cout << "      finalRange(mm)= " << finalRange/mm
     672             << ", dRoverRange= " << dRoverRange
     673             << ", integral: " << integral
     674             << ", fluct: " << lossFluctuationFlag
     675             << ", linLossLimit= " << linLossLimit
     676             << G4endl;
     677    }
     678    PrintInfo();
     679    modelManager->DumpModelList(verboseLevel);
     680    if(theCSDARangeTable && isIonisation) {
     681      G4cout << "      CSDA range table up"
     682             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
     683             << " in " << nBinsCSDA << " bins" << G4endl;
     684    }
     685    if(nSCoffRegions>0 && isIonisation) {
     686      G4cout << "      Subcutoff sampling in " << nSCoffRegions
     687             << " regions" << G4endl;
     688    }
     689    if(2 < verboseLevel) {
     690      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
     691      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
     692      G4cout << "non restricted DEDXTable address= "
     693             << theDEDXunRestrictedTable << G4endl;
     694      if(theDEDXunRestrictedTable && isIonisation) {
     695           G4cout << (*theDEDXunRestrictedTable) << G4endl;
     696      }
     697      if(theDEDXSubTable && isIonisation) {
     698        G4cout << (*theDEDXSubTable) << G4endl;
     699      }
     700      G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     701             << G4endl;
     702      if(theCSDARangeTable && isIonisation) {
     703        G4cout << (*theCSDARangeTable) << G4endl;
     704      }
     705      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     706             << G4endl;
     707      if(theRangeTableForLoss && isIonisation) {
     708             G4cout << (*theRangeTableForLoss) << G4endl;
     709      }
     710      G4cout << "      InverseRangeTable address= " << theInverseRangeTable
     711             << G4endl;
     712      if(theInverseRangeTable && isIonisation) {
     713             G4cout << (*theInverseRangeTable) << G4endl;
     714      }
     715      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     716      if(theLambdaTable && isIonisation) {
     717        G4cout << (*theLambdaTable) << G4endl;
     718      }
     719      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
     720      if(theSubLambdaTable && isIonisation) {
     721        G4cout << (*theSubLambdaTable) << G4endl;
     722      }
     723    }
     724  }
     725}
     726
     727//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     728
     729void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
     730{
     731  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     732  const G4Region* reg = r;
     733  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     734
     735  // the region is in the list
     736  if (nSCoffRegions) {
     737    for (G4int i=0; i<nSCoffRegions; i++) {
     738      if (reg == scoffRegions[i]) {
     739        if(!val) deRegions[i] = 0;
     740        return;
     741      }
     742    }
     743  }
     744
     745  // new region
     746  if(val) {
     747    useSubCutoff = true;
     748    scoffRegions.push_back(reg);
     749    nSCoffRegions++;
     750  } else {
     751    useSubCutoff = false;
     752  }
     753}
     754
     755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     756
     757void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     758{
     759  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     760  const G4Region* reg = r;
     761  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     762
     763  // the region is in the list
     764  if (nDERegions) {
     765    for (G4int i=0; i<nDERegions; i++) {
     766      if (reg == deRegions[i]) {
     767        if(!val) deRegions[i] = 0;
     768        return;
     769      }
     770    }
     771  }
     772
     773  // new region
     774  if(val) {
     775    useDeexcitation = true;
     776    deRegions.push_back(reg);
     777    nDERegions++;
     778  } else {
     779    useDeexcitation = false;
     780  }
    647781}
    648782
    … …  
    674808  //  <<" stepLimit= "<<x<<G4endl;
    675809  return x;
    676 }
    677 
    678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    679 
    680 G4double G4VEnergyLossProcess::GetMeanFreePath(
    681                              const G4Track& track,
    682                              G4double,
    683                              G4ForceCondition* condition)
    684 
    685 {
    686   *condition = NotForced;
    687   return MeanFreePath(track);
    688810}
    689811
    … …  
    806928  if (length >= fRange) {
    807929    eloss = preStepKinEnergy;
    808     currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
    809                                        eloss, esecdep, length);
     930    if (useDeexcitation) {
     931      if(idxDERegions[currentMaterialIndex]) {
     932        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     933        if(eloss < 0.0) eloss = 0.0;
     934      }
     935    }
    810936    fParticleChange.SetProposedKineticEnergy(0.0);
    811     fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
     937    fParticleChange.ProposeLocalEnergyDeposit(eloss);
    812938    return &fParticleChange;
    813939  }
    … …  
    9341060      G4double tmax =
    9351061        std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
     1062      G4double emean = eloss;
    9361063      eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
    937                                        tmax,length,eloss);
     1064                                       tmax,length,emean);
    9381065      /*                           
    9391066      if(-1 < verboseLevel)
    … …  
    9501077  eloss += esecdep;
    9511078  if(eloss < 0.0) eloss = 0.0;
     1079
     1080  // deexcitation
     1081  else if (useDeexcitation) {
     1082    if(idxDERegions[currentMaterialIndex]) {
     1083      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     1084      if(eloss < 0.0) eloss = 0.0;
     1085    }
     1086  }
    9521087
    9531088  // Energy balanse
    … …  
    11061241
    11071242  SelectModel(postStepScaledEnergy);
     1243  if(useDeexcitation) {
     1244    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     1245  }
     1246
    11081247  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    11091248  G4double tcut = (*theCuts)[currentMaterialIndex];
    … …  
    11401279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    11411280
    1142 void G4VEnergyLossProcess::PrintInfoDefinition()
    1143 {
    1144   if(0 < verboseLevel) {
    1145     G4cout << G4endl << GetProcessName() << ":   for  "
    1146            << particle->GetParticleName()
    1147            << "    SubType= " << GetProcessSubType()
    1148            << G4endl
    1149            << "      dE/dx and range tables from "
    1150            << G4BestUnit(minKinEnergy,"Energy")
    1151            << " to " << G4BestUnit(maxKinEnergy,"Energy")
    1152            << " in " << nBins << " bins" << G4endl
    1153            << "      Lambda tables from threshold to "
    1154            << G4BestUnit(maxKinEnergy,"Energy")
    1155            << " in " << nBins << " bins, spline: "
    1156            << (G4LossTableManager::Instance())->SplineFlag()
     1281G4bool G4VEnergyLossProcess::StorePhysicsTable(
     1282       const G4ParticleDefinition* part, const G4String& directory,
     1283       G4bool ascii)
     1284{
     1285  G4bool res = true;
     1286  if ( baseParticle || part != particle ) return res;
     1287
     1288  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
     1289    {res = false;}
     1290
     1291  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
     1292    {res = false;}
     1293
     1294  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
     1295    {res = false;}
     1296
     1297  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
     1298    {res = false;}
     1299
     1300  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
     1301    {res = false;}
     1302
     1303  if(isIonisation &&
     1304     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
     1305    {res = false;}
     1306
     1307  if(isIonisation &&
     1308     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
     1309    {res = false;}
     1310 
     1311  if(isIonisation &&
     1312     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
     1313    {res = false;}
     1314 
     1315  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
     1316    {res = false;}
     1317
     1318  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
     1319    {res = false;}
     1320
     1321  if ( res ) {
     1322    if(0 < verboseLevel) {
     1323      G4cout << "Physics tables are stored for " << particle->GetParticleName()
     1324             << " and process " << GetProcessName()
     1325             << " in the directory <" << directory
     1326             << "> " << G4endl;
     1327    }
     1328  } else {
     1329    G4cout << "Fail to store Physics Tables for "
     1330           << particle->GetParticleName()
     1331           << " and process " << GetProcessName()
     1332           << " in the directory <" << directory
     1333           << "> " << G4endl;
     1334  }
     1335  return res;
     1336}
     1337
     1338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1339
     1340G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
     1341       const G4ParticleDefinition* part, const G4String& directory,
     1342       G4bool ascii)
     1343{
     1344  G4bool res = true;
     1345  const G4String particleName = part->GetParticleName();
     1346
     1347  if(1 < verboseLevel) {
     1348    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
     1349           << particleName << " and process " << GetProcessName()
     1350           << "; tables_are_built= " << tablesAreBuilt
    11571351           << G4endl;
    1158     if(theRangeTableForLoss && isIonisation) {
    1159       G4cout << "      finalRange(mm)= " << finalRange/mm
    1160              << ", dRoverRange= " << dRoverRange
    1161              << ", integral: " << integral
    1162              << ", fluct: " << lossFluctuationFlag
    1163              << ", linLossLimit= " << linLossLimit
    1164              << G4endl;
    1165     }
    1166     PrintInfo();
    1167     modelManager->DumpModelList(verboseLevel);
    1168     if(theCSDARangeTable && isIonisation) {
    1169       G4cout << "      CSDA range table up"
    1170              << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
    1171              << " in " << nBinsCSDA << " bins" << G4endl;
    1172     }
    1173     if(nSCoffRegions>0 && isIonisation) {
    1174       G4cout << "      Subcutoff sampling in " << nSCoffRegions
    1175              << " regions" << G4endl;
    1176     }
    1177     if(2 < verboseLevel) {
    1178       G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
    1179       if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
    1180       G4cout << "non restricted DEDXTable address= "
    1181              << theDEDXunRestrictedTable << G4endl;
    1182       if(theDEDXunRestrictedTable && isIonisation) {
    1183            G4cout << (*theDEDXunRestrictedTable) << G4endl;
    1184       }
    1185       if(theDEDXSubTable && isIonisation) {
    1186         G4cout << (*theDEDXSubTable) << G4endl;
    1187       }
    1188       G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     1352  }
     1353  if(particle == part) {
     1354
     1355    //    G4bool yes = true;
     1356    if ( !baseParticle ) {
     1357
     1358      G4bool fpi = true;
     1359      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
     1360        {fpi = false;}
     1361
     1362      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
     1363        {fpi = false;}
     1364
     1365      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
     1366        {res = false;}
     1367
     1368      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
     1369        {res = false;}
     1370
     1371      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
     1372        {res = false;}
     1373
     1374      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
     1375        {res = false;}
     1376
     1377      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
     1378        {res = false;}
     1379
     1380      G4bool yes = false;
     1381      if(nSCoffRegions > 0) {yes = true;}
     1382
     1383      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
     1384        {res = false;}
     1385
     1386      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
     1387        {res = false;}
     1388
     1389      if(!fpi) yes = false;
     1390      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
     1391        {res = false;}
     1392    }
     1393  }
     1394
     1395  return res;
     1396}
     1397
     1398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1399
     1400G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
     1401                                        G4PhysicsTable* aTable, G4bool ascii,
     1402                                        const G4String& directory,
     1403                                        const G4String& tname)
     1404{
     1405  G4bool res = true;
     1406  if ( aTable ) {
     1407    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
     1408    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
     1409  }
     1410  return res;
     1411}
     1412
     1413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1414
     1415G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
     1416                                           G4PhysicsTable* aTable, G4bool ascii,
     1417                                           const G4String& directory,
     1418                                           const G4String& tname,
     1419                                           G4bool mandatory)
     1420{
     1421  G4bool res = true;
     1422  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
     1423  G4bool yes = aTable->ExistPhysicsTable(filename);
     1424  if(yes) {
     1425    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
     1426    if((G4LossTableManager::Instance())->SplineFlag()) {
     1427      size_t n = aTable->length();
     1428      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
     1429    }
     1430  }
     1431  if(yes) {
     1432    if (0 < verboseLevel) {
     1433      G4cout << tname << " table for " << part->GetParticleName()
     1434             << " is Retrieved from <" << filename << ">"
    11891435             << G4endl;
    1190       if(theCSDARangeTable && isIonisation) {
    1191         G4cout << (*theCSDARangeTable) << G4endl;
    1192       }
    1193       G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     1436    }
     1437  } else {
     1438    if(mandatory) res = false;
     1439    if(mandatory || 1 < verboseLevel) {
     1440      G4cout << tname << " table for " << part->GetParticleName()
     1441             << " from file <"
     1442             << filename << "> is not Retrieved"
    11941443             << G4endl;
    1195       if(theRangeTableForLoss && isIonisation) {
    1196              G4cout << (*theRangeTableForLoss) << G4endl;
    1197       }
    1198       G4cout << "      InverseRangeTable address= " << theInverseRangeTable
    1199              << G4endl;
    1200       if(theInverseRangeTable && isIonisation) {
    1201              G4cout << (*theInverseRangeTable) << G4endl;
    1202       }
    1203       G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
    1204       if(theLambdaTable && isIonisation) {
    1205         G4cout << (*theLambdaTable) << G4endl;
    1206       }
    1207       G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
    1208       if(theSubLambdaTable && isIonisation) {
    1209         G4cout << (*theSubLambdaTable) << G4endl;
    1210       }
     1444    }
     1445  }
     1446  return res;
     1447}
     1448
     1449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1450
     1451G4double G4VEnergyLossProcess::GetDEDXDispersion(
     1452                                  const G4MaterialCutsCouple *couple,
     1453                                  const G4DynamicParticle* dp,
     1454                                        G4double length)
     1455{
     1456  DefineMaterial(couple);
     1457  G4double ekin = dp->GetKineticEnergy();
     1458  SelectModel(ekin*massRatio);
     1459  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
     1460  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
     1461  G4double d = 0.0;
     1462  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
     1463  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
     1464  return d;
     1465}
     1466
     1467//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1468
     1469G4double G4VEnergyLossProcess::CrossSectionPerVolume(
     1470         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
     1471{
     1472  // Cross section per volume is calculated
     1473  DefineMaterial(couple);
     1474  G4double cross = 0.0;
     1475  G4bool b;
     1476  if(theLambdaTable) {
     1477    cross =
     1478      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
     1479  } else {
     1480    SelectModel(kineticEnergy);
     1481    cross =
     1482      currentModel->CrossSectionPerVolume(currentMaterial,
     1483                                          particle, kineticEnergy,
     1484                                          (*theCuts)[currentMaterialIndex]);
     1485  }
     1486
     1487  return cross;
     1488}
     1489
     1490//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1491
     1492G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
     1493{
     1494  DefineMaterial(track.GetMaterialCutsCouple());
     1495  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
     1496  G4double x = DBL_MAX;
     1497  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     1498  return x;
     1499}
     1500
     1501//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1502
     1503G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
     1504                                                   G4double x, G4double y,
     1505                                                   G4double& z)
     1506{
     1507  G4GPILSelection sel;
     1508  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
     1509}
     1510
     1511//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1512
     1513G4double G4VEnergyLossProcess::GetMeanFreePath(
     1514                             const G4Track& track,
     1515                             G4double,
     1516                             G4ForceCondition* condition)
     1517
     1518{
     1519  *condition = NotForced;
     1520  return MeanFreePath(track);
     1521}
     1522
     1523//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1524
     1525G4double G4VEnergyLossProcess::GetContinuousStepLimit(
     1526                const G4Track&,
     1527                G4double, G4double, G4double&)
     1528{
     1529  return DBL_MAX;
     1530}
     1531
     1532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1533
     1534G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
     1535                 const G4MaterialCutsCouple* couple, G4double cut)
     1536{
     1537  //  G4double cut  = (*theCuts)[couple->GetIndex()];
     1538  //  G4int nbins = nLambdaBins;
     1539  G4double tmin =
     1540    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
     1541             minKinEnergy);
     1542  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1543  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1544  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     1545
     1546  return v;
     1547}
     1548
     1549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1550 
     1551void G4VEnergyLossProcess::AddCollaborativeProcess(
     1552            G4VEnergyLossProcess* p)
     1553{
     1554  G4bool add = true;
     1555  if(p->GetProcessName() != "eBrem") add = false;
     1556  if(add && nProcesses > 0) {
     1557    for(G4int i=0; i<nProcesses; i++) {
     1558      if(p == scProcesses[i]) {
     1559        add = false;
     1560        break;
     1561      }
     1562    }
     1563  }
     1564  if(add) {
     1565    scProcesses.push_back(p);
     1566    nProcesses++;
     1567    if (1 < verboseLevel) {
     1568      G4cout << "### The process " << p->GetProcessName()
     1569             << " is added to the list of collaborative processes of "
     1570             << GetProcessName() << G4endl;
    12111571    }
    12121572  }
    … …  
    13771737//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    13781738
    1379 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
    1380                  const G4MaterialCutsCouple* couple, G4double cut)
    1381 {
    1382   //  G4double cut  = (*theCuts)[couple->GetIndex()];
    1383   //  G4int nbins = nLambdaBins;
    1384   G4double tmin =
    1385     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
    1386              minKinEnergy);
    1387   if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
    1388   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
    1389   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
    1390 
    1391   return v;
    1392 }
    1393 
    1394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1395 
    1396 G4double G4VEnergyLossProcess::CrossSectionPerVolume(
    1397          G4double kineticEnergy, const G4MaterialCutsCouple* couple)
    1398 {
    1399   // Cross section per volume is calculated
    1400   DefineMaterial(couple);
    1401   G4double cross = 0.0;
    1402   G4bool b;
    1403   if(theLambdaTable) {
    1404     cross =
    1405       ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
    1406   } else {
    1407     SelectModel(kineticEnergy);
    1408     cross =
    1409       currentModel->CrossSectionPerVolume(currentMaterial,
    1410                                           particle, kineticEnergy,
    1411                                           (*theCuts)[currentMaterialIndex]);
    1412   }
    1413 
    1414   return cross;
    1415 }
    1416 
    1417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1418 
    1419 G4bool G4VEnergyLossProcess::StorePhysicsTable(
    1420        const G4ParticleDefinition* part, const G4String& directory,
    1421        G4bool ascii)
    1422 {
    1423   G4bool res = true;
    1424   if ( baseParticle || part != particle ) return res;
    1425 
    1426   if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
    1427     {res = false;}
    1428 
    1429   if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
    1430     {res = false;}
    1431 
    1432   if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
    1433     {res = false;}
    1434 
    1435   if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
    1436     {res = false;}
    1437 
    1438   if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
    1439     {res = false;}
    1440 
    1441   if(isIonisation &&
    1442      !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
    1443     {res = false;}
    1444 
    1445   if(isIonisation &&
    1446      !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
    1447     {res = false;}
    1448  
    1449   if(isIonisation &&
    1450      !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
    1451     {res = false;}
    1452  
    1453   if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
    1454     {res = false;}
    1455 
    1456   if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
    1457     {res = false;}
    1458 
    1459   if ( res ) {
    1460     if(0 < verboseLevel) {
    1461       G4cout << "Physics tables are stored for " << particle->GetParticleName()
    1462              << " and process " << GetProcessName()
    1463              << " in the directory <" << directory
    1464              << "> " << G4endl;
    1465     }
    1466   } else {
    1467     G4cout << "Fail to store Physics Tables for "
    1468            << particle->GetParticleName()
    1469            << " and process " << GetProcessName()
    1470            << " in the directory <" << directory
    1471            << "> " << G4endl;
    1472   }
    1473   return res;
    1474 }
    1475 
    1476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1477 
    1478 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
    1479                                         G4PhysicsTable* aTable, G4bool ascii,
    1480                                         const G4String& directory,
    1481                                         const G4String& tname)
    1482 {
    1483   G4bool res = true;
    1484   if ( aTable ) {
    1485     const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
    1486     if( !aTable->StorePhysicsTable(name,ascii)) res = false;
    1487   }
    1488   return res;
    1489 }
    1490 
    1491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1492 
    1493 G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
    1494        const G4ParticleDefinition* part, const G4String& directory,
    1495        G4bool ascii)
    1496 {
    1497   G4bool res = true;
    1498   const G4String particleName = part->GetParticleName();
    1499 
    1500   if(1 < verboseLevel) {
    1501     G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
    1502            << particleName << " and process " << GetProcessName()
    1503            << "; tables_are_built= " << tablesAreBuilt
    1504            << G4endl;
    1505   }
    1506   if(particle == part) {
    1507 
    1508     //    G4bool yes = true;
    1509     if ( !baseParticle ) {
    1510 
    1511       G4bool fpi = true;
    1512       if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
    1513         {fpi = false;}
    1514 
    1515       if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
    1516         {fpi = false;}
    1517 
    1518       if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
    1519         {res = false;}
    1520 
    1521       if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
    1522         {res = false;}
    1523 
    1524       if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
    1525         {res = false;}
    1526 
    1527       if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
    1528         {res = false;}
    1529 
    1530       if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
    1531         {res = false;}
    1532 
    1533       G4bool yes = false;
    1534       if(nSCoffRegions > 0) {yes = true;}
    1535 
    1536       if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
    1537         {res = false;}
    1538 
    1539       if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
    1540         {res = false;}
    1541 
    1542       if(!fpi) yes = false;
    1543       if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
    1544         {res = false;}
    1545     }
    1546   }
    1547 
    1548   return res;
    1549 }
    1550 
    1551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1552 
    1553 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
    1554                                            G4PhysicsTable* aTable, G4bool ascii,
    1555                                            const G4String& directory,
    1556                                            const G4String& tname,
    1557                                            G4bool mandatory)
    1558 {
    1559   G4bool res = true;
    1560   G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
    1561   G4bool yes = aTable->ExistPhysicsTable(filename);
    1562   if(yes) {
    1563     yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
    1564     if((G4LossTableManager::Instance())->SplineFlag()) {
    1565       size_t n = aTable->length();
    1566       for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
    1567     }
    1568   }
    1569   if(yes) {
    1570     if (0 < verboseLevel) {
    1571       G4cout << tname << " table for " << part->GetParticleName()
    1572              << " is Retrieved from <" << filename << ">"
    1573              << G4endl;
    1574     }
    1575   } else {
    1576     if(mandatory) res = false;
    1577     if(mandatory || 1 < verboseLevel) {
    1578       G4cout << tname << " table for " << part->GetParticleName()
    1579              << " from file <"
    1580              << filename << "> is not Retrieved"
    1581              << G4endl;
    1582     }
    1583   }
    1584   return res;
    1585 }
    1586 
    1587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1588  
    1589 void G4VEnergyLossProcess::AddCollaborativeProcess(
    1590             G4VEnergyLossProcess* p)
    1591 {
    1592   G4bool add = true;
    1593   if(p->GetProcessName() != "eBrem") add = false;
    1594   if(add && nProcesses > 0) {
    1595     for(G4int i=0; i<nProcesses; i++) {
    1596       if(p == scProcesses[i]) {
    1597         add = false;
    1598         break;
    1599       }
    1600     }
    1601   }
    1602   if(add) {
    1603     scProcesses.push_back(p);
    1604     nProcesses++;
    1605     if (1 < verboseLevel) {
    1606       G4cout << "### The process " << p->GetProcessName()
    1607              << " is added to the list of collaborative processes of "
    1608              << GetProcessName() << G4endl;
    1609     }
    1610   }
    1611 }
    1612 
    1613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1614 
    1615 G4double G4VEnergyLossProcess::GetDEDXDispersion(
    1616                                   const G4MaterialCutsCouple *couple,
    1617                                   const G4DynamicParticle* dp,
    1618                                         G4double length)
    1619 {
    1620   DefineMaterial(couple);
    1621   G4double ekin = dp->GetKineticEnergy();
    1622   SelectModel(ekin*massRatio);
    1623   G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
    1624   tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
    1625   G4double d = 0.0;
    1626   G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
    1627   if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
    1628   return d;
    1629 }
    1630 
    1631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1632 
    1633 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
    1634 {}
    1635 
    1636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1637 

  • trunk/source/processes/electromagnetic/utils/src/G4VMscModel.cc

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.cc,v 1.4 2008/03/10 18:39:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VMscModel.cc,v 1.10 2009/02/24 09:56:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4949
    5050#include "G4VMscModel.hh"
     51#include "G4ParticleChangeForMSC.hh"
     52#include "G4TransportationManager.hh"
    5153
    5254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5557G4VMscModel::G4VMscModel(const G4String& nam):
    5658  G4VEmModel(nam),
     59  safetyHelper(0),
    5760  facrange(0.02),
    5861  facgeom(2.5),
    … …  
    6164  dtrl(0.05),
    6265  lambdalimit(mm),
     66  geommax(1.e50*mm),
    6367  steppingAlgorithm(fUseSafety),
    6468  samplez(false),
    … …  
    7276
    7377//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     78
     79void G4VMscModel::InitialiseSafetyHelper()
     80{
     81  if(!safetyHelper) {
     82    safetyHelper = G4TransportationManager::GetTransportationManager()
     83      ->GetSafetyHelper();
     84    safetyHelper->InitialiseHelper();
     85  }
     86}
     87
     88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     89
     90void G4VMscModel::ComputeDisplacement(G4ParticleChangeForMSC* fParticleChange, 
     91                                      const G4ThreeVector& dir,
     92                                      G4double displacement,
     93                                      G4double postsafety)
     94{
     95  const G4ThreeVector* pos = fParticleChange->GetProposedPosition();
     96  G4double r = displacement;
     97  if(r >  postsafety) {
     98    G4double newsafety = safetyHelper->ComputeSafety(*pos);
     99    if(r > newsafety) r = newsafety;
     100  }
     101  if(r > 0.) {
     102
     103    // compute new endpoint of the Step
     104    G4ThreeVector newPosition = *pos + r*dir;
     105
     106    // definitely not on boundary
     107    if(displacement == r) {
     108      safetyHelper->ReLocateWithinVolume(newPosition);
     109
     110    } else {
     111      // check safety after displacement
     112      G4double postsafety = safetyHelper->ComputeSafety(newPosition);
     113
     114      // displacement to boundary
     115      if(postsafety <= 0.0) {
     116        safetyHelper->Locate(newPosition,
     117                             *fParticleChange->GetProposedMomentumDirection());
     118
     119        // not on the boundary
     120      } else {
     121        safetyHelper->ReLocateWithinVolume(newPosition);
     122      }
     123    }
     124    fParticleChange->ProposePosition(newPosition);
     125  }
     126}
     127
     128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     129
     130void G4VMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
     131                                    const G4MaterialCutsCouple*,
     132                                    const G4DynamicParticle*,
     133                                    G4double, G4double)
     134{}
     135
     136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracChangeset for help on using the changeset viewer.