Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (15 years ago)
Author:
garnier
Message:

update ti head

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

Legend:

Unmodified
Added
Removed

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ElectronIonPair.cc,v 1.2 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ElectronIonPair.cc,v 1.5 2010/10/25 17:23:01 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646
    4747#include "G4ElectronIonPair.hh"
    48 #include "G4Gamma.hh"
    4948#include "G4Material.hh"
    5049#include "G4MaterialTable.hh"
    … …  
    6362  curMeanEnergy = 0.0;
    6463  nMaterials = 0;
     64  FanoFactor = 0.2;
    6565  Initialise();
    6666}
    … …  
    9797        }
    9898      }
    99       if(curMeanEnergy > 0.0) nion = (edep - niel)/curMeanEnergy;
     99      if(curMeanEnergy > 0.0) { nion = (edep - niel)/curMeanEnergy; }
    100100    }
    101101  }
    … …  
    106106
    107107std::vector<G4ThreeVector>*
    108 G4ElectronIonPair::SampleIonsAlongStep(const G4ThreeVector& prePos,
    109                                        const G4ThreeVector& postPos,
    110                                        G4double meanion)
    111 {
    112   std::vector<G4ThreeVector>* v = new std::vector<G4ThreeVector>;
    113 
    114   G4double sig = 0.2*std::sqrt(meanion);
    115   G4int nion = G4int(G4RandGauss::shoot(meanion,sig) + 0.5);
     108G4ElectronIonPair::SampleIonsAlongStep(const G4Step* step)
     109{
     110  std::vector<G4ThreeVector>* v = 0;
     111
     112  G4int nion = SampleNumberOfIonsAlongStep(step);
    116113
    117114  // sample ionisation along step
    118115  if(nion > 0) {
    119116
    120     G4ThreeVector deltaPos = postPos - prePos; 
    121     for(G4int i=0; i<nion; i++) {
     117    v = new std::vector<G4ThreeVector>;
     118    G4ThreeVector prePos = step->GetPreStepPoint()->GetPosition();
     119    G4ThreeVector deltaPos = step->GetPostStepPoint()->GetPosition() - prePos; 
     120    for(G4int i=0; i<nion; ++i) {
    122121      v->push_back( prePos + deltaPos*G4UniformRand() );
    123122    }
    … …  
    138137  G4int nholes = 0;
    139138
    140   if(2 == subType || 12 == subType || 13 == subType) nholes = 1;
     139  if(2 == subType || 12 == subType || 13 == subType) { nholes = 1; }
    141140  return nholes;
    142141}
    … …  
    173172  if(nmat > 0) {
    174173    G4cout << "### G4ElectronIonPair: mean energy per ion pair avalable:" << G4endl;
    175     for(G4int i=0; i<nmat; i++) {
     174    for(G4int i=0; i<nmat; ++i) {
    176175      const G4Material* mat = (*mtable)[i];
    177176      G4double x = mat->GetIonisation()->GetMeanEnergyPerIonPair();
    … …  
    191190    G4cout << "### G4ElectronIonPair: mean energy per ion pair "
    192191           << " for Geant4 materials" << G4endl;
    193     for(G4int i=0; i<nMaterials; i++) {
     192    for(G4int i=0; i<nMaterials; ++i) {
    194193      G4cout << "   " << g4MatNames[i] << "    Epair= "
    195194             << g4MatData[i]/eV << " eV" << G4endl;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.cc,v 1.53 2010/04/13 10:58:03 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmCalculator.cc,v 1.57 2010/11/04 12:55:09 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9898  currentLambda      = 0;
    9999  currentModel       = 0;
     100  currentProcess     = 0;
    100101  loweModel          = 0;
    101102  chargeSquare       = 1.0;
    102103  massRatio          = 1.0;
     104  mass               = 0.0;
     105  currentCut         = 0.0;
    103106  currentParticleName= "";
    104107  currentMaterialName= "";

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.cc,v 1.8 2010/06/04 15:33:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmConfigurator.cc,v 1.9 2010/07/29 11:13:28 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    172172
    173173      G4VProcess* proc = 0;
    174       for(G4int i=0; i<np; i++) {
     174      for(G4int i=0; i<np; ++i) {
    175175        if(processName == (*plist)[i]->GetProcessName()) {
    176176          proc = (*plist)[i];
    … …  
    182182               << processName << "> for " << particleName << G4endl;
    183183        return;
    184 
    185184      }
    186185
    187186      if(mod) {
    188         if(!UpdateModelEnergyRange(mod, emin,emax)) { return; }
     187        if(!UpdateModelEnergyRange(mod, emin, emax)) { return; }
    189188      }
    190189      // classify process

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.cc,v 1.60 2010/04/12 18:28:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmModelManager.cc,v 1.63 2010/10/15 10:22:13 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    463463      }
    464464    }
    465     /*
    466     G4int nm = setOfRegionModels[reg]->NumberOfModels();
    467     for(G4int j=0; j<nm; ++j) {
    468 
    469       G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
    470       G4double tcutmin = model->MinEnergyCut(particle, couple);
    471 
    472       if(cut < tcutmin)    { cut = tcutmin; }
    473       if(subcut < tcutmin) { subcut = tcutmin; }
    474       if(1 < verboseLevel) {
    475             G4cout << "The model # " << j
    476                    << "; tcutmin(MeV)= " << tcutmin/MeV
    477                    << "; tcut(MeV)= " << cut/MeV
    478                    << "; tsubcut(MeV)= " << subcut/MeV
    479                    << " for " << particle->GetParticleName()
    480                    << G4endl;
    481       }
    482     }
    483     */
    484465    theCuts[i] = cut;
    485466    if(minSubRange < 1.0) { theSubCuts[i] = subcut; }
    … …  
    656637void G4EmModelManager::DumpModelList(G4int verb)
    657638{
    658   if(verb == 0) return;
     639  if(verb == 0) { return; }
    659640  for(G4int i=0; i<nRegions; ++i) {
    660641    G4RegionModels* r = setOfRegionModels[i];
    … …  
    665646             << " ======" << G4endl;;
    666647      for(G4int j=0; j<n; ++j) {
    667         const G4VEmModel* m = models[r->ModelIndex(j)];
     648        G4VEmModel* m = models[r->ModelIndex(j)];
    668649        G4cout << std::setw(20);
    669         G4cout << m->GetName() << " :     Emin= "
    670                << std::setw(10) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
    671                << "        Emax=   "
    672                << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy")
    673                << G4endl;
     650        G4cout << m->GetName() << " :   Emin= "
     651               << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
     652               << "      Emax=   "
     653               << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy");
     654        G4VEmAngularDistribution* an = m->GetAngularDistribution();
     655        if(an) { G4cout << "  " << an->GetName(); }
     656        G4cout << G4endl;
    674657      } 
    675658    }

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmMultiModel.cc,v 1.6 2007/05/22 17:31:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmMultiModel.cc,v 1.8 2010/08/17 17:36:59 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4040// Modifications:
    4141// 15-04-05 optimize internal interface (V.Ivanchenko)
    42 //
    43 
    44 // Class Description:
    45 //
    46 // Energy loss model using several G4VEmModels
    47 
    48 // -------------------------------------------------------------------
     42// 04-07-10 updated interfaces according to g4 9.4 (V.Ivanchenko)
    4943//
    5044
    … …  
    5953
    6054G4EmMultiModel::G4EmMultiModel(const G4String& nam)
    61   : G4VEmModel(nam),
    62   nModels(0)
     55  : G4VEmModel(nam), nModels(0)
    6356{
    6457  model.clear();
    65   tsecmin.clear();
    6658  cross_section.clear();
    6759}
    … …  
    7062
    7163G4EmMultiModel::~G4EmMultiModel()
     64{}
     65
     66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     67
     68void G4EmMultiModel::AddModel(G4VEmModel* p)
    7269{
    73   if(nModels) {
    74     for(G4int i=0; i<nModels; i++) {
    75       delete model[i];
    76     }
    77   }
     70  cross_section.push_back(0.0);
     71  model.push_back(p);
     72  ++nModels;
    7873}
    7974
    … …  
    8378                                const G4DataVector& cuts)
    8479{
    85   if(nModels) {
    86     for(G4int i=0; i<nModels; i++) {
     80  if(nModels > 0) {
     81    G4cout << "### Initialisation of EM MultiModel " << GetName()
     82           << " including following list of models:" << G4endl;
     83    for(G4int i=0; i<nModels; ++i) {
     84      G4cout << "    " << (model[i])->GetName();
     85      (model[i])->SetParticleChange(pParticleChange, GetModelOfFluctuations());
    8786      (model[i])->Initialise(p, cuts);
    8887    }
     88    G4cout << G4endl;
    8989  }
    90 }
    91 
    92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    93 
    94 G4double G4EmMultiModel::MinEnergyCut(const G4ParticleDefinition* p,
    95                                       const G4MaterialCutsCouple* couple)
    96 {
    97   G4double cut = DBL_MAX;
    98   if(nModels) {
    99     cut = (model[0])->MinEnergyCut(p, couple);
    100   }
    101   return cut;
    10290}
    10391
    … …  
    10694G4double G4EmMultiModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
    10795                                     const G4ParticleDefinition* p,
    108                                            G4double kineticEnergy,
    109                                            G4double cutEnergy)
     96                                     G4double kineticEnergy,
     97                                     G4double cutEnergy)
    11098{
    111   G4double dedx  = 0.0;
     99  SetCurrentCouple(couple);
     100  G4double dedx = 0.0;
    112101
    113   if(nModels) {
    114     dedx =  (model[0])->ComputeDEDX(couple, p, cutEnergy, kineticEnergy);
     102  if(nModels > 0) {
     103    for(G4int i=0; i<nModels; i++) {
     104      dedx += (model[i])->ComputeDEDX(couple, p, cutEnergy, kineticEnergy);
     105    }
    115106  }
    116107
    … …  
    120111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    121112
    122 G4double G4EmMultiModel::CrossSection(const G4MaterialCutsCouple* couple,
    123                                       const G4ParticleDefinition* p,
    124                                             G4double kineticEnergy,
    125                                             G4double cutEnergy,
    126                                             G4double maxKinEnergy)
     113G4double G4EmMultiModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
     114                                                    G4double kinEnergy,
     115                                                    G4double Z,
     116                                                    G4double A,
     117                                                    G4double cutEnergy,
     118                                                    G4double maxEnergy)
    127119{
    128   G4double cross     = 0.0;
    129   G4double t1        = cutEnergy;
    130   G4double t2        = cutEnergy;
    131   if(nModels) {
    132     for(G4int i=0; i<nModels; i++) {
    133       t1 = std::max(t2, tsecmin[i]);
    134       t2 = std::min(maxKinEnergy, tsecmin[i+1]);
    135       cross += (model[i])->CrossSection(couple, p, kineticEnergy, t1, t2);
     120  G4double cross = 0.0;
     121  if(nModels>0) {
     122    for(G4int i=0; i<nModels; ++i) {
     123      (model[i])->SetCurrentCouple(CurrentCouple());
     124      cross += (model[i])->ComputeCrossSectionPerAtom(p, kinEnergy, Z, A,
     125                                                      cutEnergy, maxEnergy);
    136126    }
    137127  }
    … …  
    144134                                       const G4MaterialCutsCouple* couple,
    145135                                       const G4DynamicParticle* dp,
    146                                        G4double tmin,
     136                                       G4double minEnergy,
    147137                                       G4double maxEnergy)
    148138{
     139  SetCurrentCouple(couple);
    149140  if(nModels > 0) {
    150141    G4int i;
    151142    G4double cross = 0.0;
    152     G4double t1    = tmin;
    153     G4double t2    = tmin;
    154     for(i=0; i<nModels; i++) {
    155       t1 = std::max(t2, tsecmin[i]);
    156       t2 = std::min(maxEnergy, tsecmin[i+1]);
    157       cross += (model[i])->CrossSection(couple, dp->GetDefinition(),
    158                                         dp->GetKineticEnergy(), t1, t2);
     143    for(i=0; i<nModels; ++i) {
     144      cross += (model[i])->CrossSection(couple, dp->GetParticleDefinition(),
     145                                        dp->GetKineticEnergy(), minEnergy, maxEnergy);
    159146      cross_section[i] = cross;
    160147    }
    161148
    162149    cross *= G4UniformRand();
    163     t2 = tmin;
    164150
    165     for(i=0; i<nModels; i++) {
    166       t1 = std::max(t2, tsecmin[i]);
    167       t2 = std::min(maxEnergy, tsecmin[i+1]);
     151    for(i=0; i<nModels; ++i) {
    168152      if(cross <= cross_section[i]) {
    169         (model[i])->SampleSecondaries(vdp, couple, dp, t1, t2);
    170         break;
     153        (model[i])->SampleSecondaries(vdp, couple, dp, minEnergy, maxEnergy);
     154        return;
    171155      }
    172156    }
    … …  
    176160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    177161
    178 G4double G4EmMultiModel:: MaxSecondaryEnergy(const G4ParticleDefinition*,
    179                                                    G4double kinEnergy)
    180 {
    181   G4cout << "Warning! G4EmMultiModel::"
    182          << "MaxSecondaryEnergy(const G4ParticleDefinition*,G4double kinEnergy)"
    183          << " should not be used!" << G4endl;
    184   return kinEnergy;
    185 }
    186 
    187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    188 
    189 void G4EmMultiModel::DefineForRegion(const G4Region* r)
    190 {
    191   if(nModels) {
    192     for(G4int i=0; i<nModels; i++) {(model[i])->DefineForRegion(r);}
    193   }
    194 }
    195 
    196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    197 
    198 void G4EmMultiModel::AddModel(G4VEmModel* p, G4double tmin, G4double tmax)
    199 {
    200   if(tmin < tmax && 0.0 < tmin) {
    201 
    202     if(nModels == 0) {
    203       tsecmin.push_back(tmin);
    204       tsecmin.push_back(tmax);
    205       cross_section.push_back(0.0);
    206       model.push_back(p);
    207       nModels++;
    208 
    209     } else {
    210       G4int i, j;
    211       G4bool increment = false;
    212       for(i=0; i<nModels; i++) {
    213 
    214         if(tmin < tsecmin[i]) {
    215           G4double t2 = std::min(tsecmin[i+1],tmax);
    216           if(tmin < t2) {
    217             tsecmin.push_back(0.0);
    218             cross_section.push_back(0.0);
    219             model.push_back(0);
    220             for(j=nModels; j>i; j--) {
    221               model[j] = model[j-1];
    222               tsecmin[j+1] = tsecmin[j];
    223             }
    224             model[i] = p;
    225             tsecmin[i+1] = t2;
    226             tsecmin[i]   = tmin;
    227             increment = true;
    228           }
    229         } else if(i == nModels-1) {
    230           G4double t1 = std::min(tsecmin[i+1],tmin);
    231           G4double t2 = std::max(tsecmin[i+1],tmax);
    232           if(t1 < t2) {
    233             tsecmin.push_back(t2);
    234             cross_section.push_back(0.0);
    235             model.push_back(p);
    236             increment = true;
    237           }
    238         }
    239       }
    240       if(increment) nModels++;
    241     }
    242   }
    243 }
    244 
    245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    246 

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmSaturation.cc,v 1.10 2009/09/25 09:16:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmSaturation.cc,v 1.11 2010/10/25 17:23:01 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5959  verbose = 1;
    6060  manager = 0;
     61
    6162  curMaterial = 0;
    62   curBirks = 0.0;
    63   curRatio = 1.0;
     63  curBirks    = 0.0;
     64  curRatio    = 1.0;
    6465  curChargeSq = 1.0;
    65   nMaterials = 0;
     66  nMaterials  = 0;
     67
    6668  electron = 0;
     69  proton   = 0;
     70  nist     = G4NistManager::Instance();
     71
    6772  Initialise();
    6873}
    … …  
    8287                                      G4double niel)
    8388{
    84   if(edep <= 0.0) return 0.0;
     89  if(edep <= 0.0) { return 0.0; }
    8590
    8691  G4double evis = edep;
    … …  
    109114
    110115      // continues energy loss
    111       if(eloss > 0.0) eloss /= (1.0 + bfactor*eloss/length);
     116      if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
    112117 
    113118      // non-ionizing energy loss
    114119      if(nloss > 0.0) {
    115         if(!proton) {proton = G4Proton::Proton();}
     120        if(!proton) { proton = G4Proton::Proton(); }
    116121        G4double escaled = nloss*curRatio;
    117122        G4double s = manager->GetRange(proton,escaled,couple)/curChargeSq;
    … …  
    133138  // is this material in the vector?
    134139 
    135   for(G4int j=0; j<nG4Birks; j++) {
     140  for(G4int j=0; j<nG4Birks; ++j) {
    136141    if(name == g4MatNames[j]) {
    137142      if(verbose > 0)
    … …  
    152157  if(!manager) {
    153158    manager = G4LossTableManager::Instance();
    154     nist    = G4NistManager::Instance();
    155159    electron= G4Electron::Electron();
    156     proton  = 0;
    157   }
    158 
    159   if(mat == curMaterial) return curBirks;
     160  }
     161
     162  if(mat == curMaterial) { return curBirks; }
    160163
    161164  curMaterial = mat;
    … …  
    165168
    166169  // seach in the run-time list
    167   for(G4int i=0; i<nMaterials; i++) {
     170  for(G4int i=0; i<nMaterials; ++i) {
    168171    if(mat == matPointers[i]) {
    169172      curBirks = mat->GetIonisation()->GetBirksConstant();
    … …  
    180183  // seach in the Geant4 list
    181184  if(curBirks == 0.0) {
    182     for(G4int j=0; j<nG4Birks; j++) {
     185    for(G4int j=0; j<nG4Birks; ++j) {
    183186      if(name == g4MatNames[j]) {
    184187        mat->GetIonisation()->SetBirksConstant(g4MatData[j]);
    … …  
    201204  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
    202205  size_t nelm = mat->GetNumberOfElements();
    203   for (size_t i=0; i<nelm; i++) {
     206  for (size_t i=0; i<nelm; ++i) {
    204207    const G4Element* elm = (*theElementVector)[i];
    205208    G4double Z = elm->GetZ();
    … …  
    231234  if(nMaterials > 0) {
    232235    G4cout << "### Birks coeffitients used in run time" << G4endl;
    233     for(G4int i=0; i<nMaterials; i++) {
     236    for(G4int i=0; i<nMaterials; ++i) {
    234237      G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
    235238      G4cout << "   " << matNames[i] << "     "
    … …  
    248251  if(nG4Birks > 0) {
    249252    G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
    250     for(G4int i=0; i<nG4Birks; i++) {
     253    for(G4int i=0; i<nG4Birks; ++i) {
    251254      G4cout << "   " << g4MatNames[i] << "   "
    252255             << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.cc,v 1.101 2010/06/04 15:33:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4LossTableManager.cc,v 1.105 2010/11/04 12:55:09 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9595#include "G4EmSaturation.hh"
    9696#include "G4EmConfigurator.hh"
     97#include "G4ElectronIonPair.hh"
    9798#include "G4EmTableType.hh"
    9899#include "G4LossTableBuilder.hh"
     100#include "G4VAtomDeexcitation.hh"
    99101#include "G4Region.hh"
    100102
    … …  
    140142  delete emCorrections;
    141143  delete emSaturation;
     144  delete emConfigurator;
     145  delete emElectronIonPair;
     146  delete atomDeexcitation;
    142147}
    143148
    … …  
    182187  emSaturation = new G4EmSaturation();
    183188  emConfigurator = new G4EmConfigurator(verbose);
     189  emElectronIonPair = new G4ElectronIonPair();
    184190  tableBuilder->SetSplineFlag(splineFlag);
     191  atomDeexcitation = 0;
    185192}
    186193
    … …  
    438445    emConfigurator->Clear();
    439446    firstParticle = aParticle;
     447  }
     448  if(startInitialisation && atomDeexcitation) {
     449    atomDeexcitation->InitialiseAtomicDeexcitation();
    440450  }
    441451  startInitialisation = false;
    … …  
    590600
    591601  G4int n_dedx = t_list.size();
    592   if (!n_dedx) {
     602  if (0 == n_dedx || !em) {
    593603    G4cout << "G4LossTableManager WARNING: no DEDX processes for "
    594604           << aParticle->GetParticleName() << G4endl;
    … …  
    896906  //emCorrections->SetVerbose(val);
    897907  emSaturation->SetVerbose(val);
     908  emElectronIonPair->SetVerbose(val);
     909  if(atomDeexcitation) { atomDeexcitation->SetVerboseLevel(val); }
    898910}
    899911
    … …  
    10441056}
    10451057
     1058//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     1059
     1060G4ElectronIonPair* G4LossTableManager::ElectronIonPair()
     1061{
     1062  return emElectronIonPair;
     1063}
     1064
    10461065//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1066
     1067G4VAtomDeexcitation* G4LossTableManager::AtomDeexcitation()
     1068{
     1069  return atomDeexcitation;
     1070}
     1071
     1072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1073 
     1074void G4LossTableManager::SetAtomDeexcitation(G4VAtomDeexcitation* p)
     1075{
     1076  atomDeexcitation = p;
     1077}
     1078//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.cc,v 1.35 2010/05/26 18:06:25 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmModel.cc,v 1.37 2010/10/14 16:27:35 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6161
    6262G4VEmModel::G4VEmModel(const G4String& nam):
    63   fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
     63  flucModel(0),anglModel(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
    6464  eMinActive(0.0),eMaxActive(DBL_MAX),
    6565  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
    66   pParticleChange(0),nuclearStopping(false),
     66  pParticleChange(0),/*nuclearStopping(false),*/
    6767  currentCouple(0),currentElement(0),
    6868  nsec(5),flagDeexcitation(false)
    … …  
    8484    }
    8585  }
     86  delete anglModel;
    8687}
    8788
    … …  
    246247G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
    247248{
    248   return GetChargeSquareRatio(track.GetDefinition(),
     249  return GetChargeSquareRatio(track.GetParticleDefinition(),
    249250                              track.GetMaterial(), track.GetKineticEnergy());
    250251}
    … …  
    296297
    297298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     299
     300void
     301G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
     302{
     303  if(p && pParticleChange != p) { pParticleChange = p; }
     304  flucModel = f;
     305}
     306
     307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.cc,v 1.87 2010/05/10 13:35:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmProcess.cc,v 1.88 2010/08/17 17:36:59 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    516516    for (G4int i=0; i<num; ++i) {
    517517      G4DynamicParticle* dp = secParticles[i];
    518       const G4ParticleDefinition* p = dp->GetDefinition();
     518      const G4ParticleDefinition* p = dp->GetParticleDefinition();
    519519      G4double e = dp->GetKineticEnergy();
    520520      G4bool good = true;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.167 2010/05/26 10:41:34 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.171 2010/10/15 10:22:13 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    108108// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
    109109// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
     110// 15-10-10 Fixed 4-momentum balance if deexcitation is active (L.Pandola)
    110111//
    111112// Class Description:
    … …  
    144145#include "G4TransportationManager.hh"
    145146#include "G4EmConfigurator.hh"
     147#include "G4VAtomDeexcitation.hh"
    146148
    147149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    222224  (G4LossTableManager::Instance())->Register(this);
    223225  fluctModel = 0;
     226  atomDeexcitation = 0;
    224227
    225228  scTracks.reserve(5);
    … …  
    570573
    571574  // Added tracking cut to avoid tracking artifacts
    572   if(isIonisation) { fParticleChange.SetLowEnergyLimit(lowestKinEnergy); }
     575  if(isIonisation) {
     576    fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
     577    atomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
     578    if(atomDeexcitation) { useDeexcitation = true; }
     579  }
    573580
    574581  if(1 < verboseLevel) {
    … …  
    852859    for (G4int i=0; i<nDERegions; ++i) {
    853860      if (reg == deRegions[i]) {
    854         if(!val) deRegions[i] = 0;
     861        if(!val) { deRegions[i] = 0; }
    855862        return;
    856863      }
    … …  
    911918  DefineMaterial(track.GetMaterialCutsCouple());
    912919
    913   const G4ParticleDefinition* currPart = track.GetDefinition();
     920  const G4ParticleDefinition* currPart = track.GetParticleDefinition();
    914921  if(theGenericIon == particle) {
    915922    massRatio = proton_mass_c2/currPart->GetPDGMass();
    … …  
    9991006  /* 
    10001007  if(-1 < verboseLevel) {
    1001     const G4ParticleDefinition* d = track.GetDefinition();
     1008    const G4ParticleDefinition* d = track.GetParticleDefinition();
    10021009    G4cout << "AlongStepDoIt for "
    10031010           << GetProcessName() << " and particle "
    … …  
    10191026    eloss = preStepKinEnergy;
    10201027    if (useDeexcitation) {
    1021       if(idxDERegions[currentMaterialIndex]) {
     1028      if(atomDeexcitation) {
     1029        atomDeexcitation->AlongStepDeexcitation(&fParticleChange, step,
     1030                                                eloss, currentMaterialIndex);
     1031      } else if(idxDERegions[currentMaterialIndex]) {
    10221032        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
    1023         if(eloss < 0.0) eloss = 0.0;
    1024       }
     1033      }
     1034      if(eloss < 0.0) { eloss = 0.0; }
    10251035    }
    10261036    fParticleChange.SetProposedKineticEnergy(0.0);
    … …  
    11281138            G4Track* t = scTracks[i];
    11291139            G4double e = t->GetKineticEnergy();
    1130             if (t->GetDefinition() == thePositron) { e += 2.0*electron_mass_c2; }
     1140            if (t->GetParticleDefinition() == thePositron) {
     1141              e += 2.0*electron_mass_c2;
     1142            }
    11311143            esec += e;
    11321144            pParticleChange->AddSecondary(t);
    … …  
    11691181  // deexcitation
    11701182  if (useDeexcitation) {
    1171     if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
     1183    G4double eloss_before = eloss;
     1184    if(atomDeexcitation) {
     1185      atomDeexcitation->AlongStepDeexcitation(&fParticleChange, step,
     1186                                              eloss, currentMaterialIndex);
     1187    } else if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
    11721188      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
    1173       if(eloss < 0.0) { eloss = 0.0; }
    1174     }
     1189    }
     1190    esec += eloss_before - eloss;
    11751191  }
    11761192
    … …  
    11781194  G4double finalT = preStepKinEnergy - eloss - esec;
    11791195  if (finalT <= lowestKinEnergy) {
    1180     eloss  = preStepKinEnergy - esec;
     1196    eloss += finalT;
    11811197    finalT = 0.0;
    11821198  } else if(isIon) {
    11831199    fParticleChange.SetProposedCharge(
    1184       currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
    1185   }
    1186 
     1200      currentModel->GetParticleCharge(track.GetParticleDefinition(),
     1201                                      currentMaterial,finalT));
     1202  }
     1203
     1204  if(eloss < 0.0) { eloss = 0.0; }
    11871205  fParticleChange.SetProposedKineticEnergy(finalT);
    11881206  fParticleChange.ProposeLocalEnergyDeposit(eloss);
    … …  
    12601278      /*
    12611279      // do not track very low-energy delta-electrons
    1262       if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
     1280      if(theSecondaryRangeTable && (*it)->GetParticleDefinition() == theElectron) {
    12631281        G4double ekin = (*it)->GetKineticEnergy();
    12641282        G4double rg = ((*theSecondaryRangeTable)[idx]->Value(ekin));
    … …  
    12781296        /*     
    12791297        if(-1 < verboseLevel)
    1280           G4cout << "New track " << t->GetDefinition()->GetParticleName()
     1298          G4cout << "New track " << t->GetParticleDefinition()->GetParticleName()
    12811299                 << " e(keV)= " << t->GetKineticEnergy()/keV
    12821300                 << " fragment= " << fragment
    … …  
    12931311                                                      const G4Step&)
    12941312{
     1313  // In all cases clear number of interaction lengths
     1314  theNumberOfInteractionLengthLeft = -1.0;
     1315
    12951316  fParticleChange.InitializeForPostStep(track);
    12961317  G4double finalT = track.GetKineticEnergy();
    1297   if(finalT <= lowestKinEnergy) return &fParticleChange;
     1318  if(finalT <= lowestKinEnergy) { return &fParticleChange; }
    12981319
    12991320  G4double postStepScaledEnergy = finalT*massRatio;
    … …  
    13211342    }
    13221343    */
    1323     if(preStepLambda*G4UniformRand() > lx) {
    1324       ClearNumberOfInteractionLengthLeft();
     1344    if(lx <= 0.0) {
     1345      return &fParticleChange;
     1346    } else if(preStepLambda*G4UniformRand() > lx) {
    13251347      return &fParticleChange;
    13261348    }
    … …  
    13601382  }
    13611383  */
    1362   ClearNumberOfInteractionLengthLeft();
    13631384  return &fParticleChange;
    13641385}

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.cc,v 1.16 2010/04/06 09:24:46 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VMscModel.cc,v 1.18 2010/09/07 16:05:33 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6464  dtrl(0.05),
    6565  lambdalimit(mm),
    66   geommax(1.e50*mm),
     66  geomMin(1.e-6*CLHEP::mm),
     67  geomMax(1.e50*CLHEP::mm),
    6768  steppingAlgorithm(fUseSafety),
    6869  samplez(false),
    … …  
    106107                                      G4double postsafety)
    107108{
     109  if(displacement <= geomMin) { return; }
    108110  const G4ThreeVector* pos = fParticleChange->GetProposedPosition();
    109   G4double r = displacement;
    110   if(r >  postsafety) {
    111     G4double newsafety = safetyHelper->ComputeSafety(*pos);
    112     if(r > newsafety) r = newsafety;
     111
     112  // displaced point is definitely within the volume
     113  if(displacement < postsafety) {
     114
     115    // compute new endpoint of the Step
     116    G4ThreeVector newPosition = *pos + displacement*dir;
     117    safetyHelper->ReLocateWithinVolume(newPosition);
     118    fParticleChange->ProposePosition(newPosition);
     119    return;
    113120  }
    114   if(r > 0.) {
     121
     122  // displaced point may be outside the volume
     123  G4double newsafety = safetyHelper->ComputeSafety(*pos);
     124
     125  // add a factor which ensure numerical stability
     126  G4double r = std::min(displacement, newsafety*0.99);
     127
     128  if(r > geomMin) {
    115129
    116130    // compute new endpoint of the Step
    117131    G4ThreeVector newPosition = *pos + r*dir;
    118 
    119     // definitely not on boundary
    120     if(displacement == r) {
    121       safetyHelper->ReLocateWithinVolume(newPosition);
    122 
    123     } else {
    124       // check safety after displacement
    125       G4double postsafety = safetyHelper->ComputeSafety(newPosition);
    126 
    127       // displacement to boundary
    128       if(postsafety <= 0.0) {
    129         safetyHelper->Locate(newPosition,
    130                              *fParticleChange->GetProposedMomentumDirection());
    131 
    132         // not on the boundary
    133       } else {
    134         safetyHelper->ReLocateWithinVolume(newPosition);
    135       }
    136     }
     132    safetyHelper->ReLocateWithinVolume(newPosition);
    137133    fParticleChange->ProposePosition(newPosition);
    138134  }

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.cc,v 1.82 2010/04/12 11:45:03 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VMultipleScattering.cc,v 1.86 2010/10/26 11:30:46 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    294294  if(verboseLevel>0 && ( num == "e-" || num == "mu+" || 
    295295                         num == "proton" || num == "pi-" ||
    296                          num == "GenericIon")) {
     296                         num == "kaon+" || num == "GenericIon")) {
    297297    PrintInfoDefinition();
    298298    if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
    … …  
    348348  DefineMaterial(track.GetMaterialCutsCouple());
    349349  G4double ekin = track.GetKineticEnergy();
    350   if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); }
     350  if(isIon) { ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass(); }
    351351  currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
    352352  if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
    353353    G4double tPathLength =
    354354      currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
    355     if (tPathLength < x) *selection = CandidateForSelection;
     355    if (tPathLength < x) { *selection = CandidateForSelection; }
    356356    x = currentModel->ComputeGeomPathLength(tPathLength);
    357357    //  G4cout << "tPathLength= " << tPathLength
    … …  
    404404                                       G4double& currentSafety)
    405405{
    406   G4GPILSelection* selection = 0;
    407   return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
    408                                                currentSafety, selection);
     406  G4GPILSelection selection = NotCandidateForSelection;
     407  G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize,
     408                                                     currentMinimalStep,
     409                                                     currentSafety, &selection);
     410  return x;
    409411}
    410412
Note: See TracChangeset for help on using the changeset viewer.