Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/electromagnetic/muons/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EnergyLossForExtrapolator.cc,v 1.13 2007/07/28 13:44:25 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4EnergyLossForExtrapolator.cc,v 1.18 2008/11/13 14:14:07 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929//---------------------------------------------------------------------------
    … …  
    6767#include "G4MuBremsstrahlungModel.hh"
    6868#include "G4ProductionCuts.hh"
     69#include "G4LossTableManager.hh"
     70#include "G4WentzelVIModel.hh"
    6971
    7072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    8890  delete invRangePositron;
    8991  delete invRangeProton;
     92  delete mscElectron;
    9093  delete cuts;
    9194}
    … …  
    100103  if(!isInitialised) Initialisation();
    101104  G4double kinEnergyFinal = kinEnergy;
    102   if(mat && part) {
    103     G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     105  if(SetupKinematics(part, mat, kinEnergy)) {
     106    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
    104107    G4double r  = ComputeRange(kinEnergy,part);
    105108    if(r <= step) {
    … …  
    125128  G4double kinEnergyFinal = kinEnergy;
    126129
    127   if(mat && part) {
    128     G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     130  if(SetupKinematics(part, mat, kinEnergy)) {
     131    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
    129132    G4double r  = ComputeRange(kinEnergy,part);
    130133
    … …  
    141144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    142145
    143 G4double G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat,
    144                                                       const G4ParticleDefinition* part,
    145                                                       G4double kinEnergy, G4double stepLength)
    146 {
     146G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
     147                                                     G4double stepLength,
     148                                                     const G4Material* mat,
     149                                                     const G4ParticleDefinition* part)
     150{
     151  G4double res = stepLength;
     152  if(!isInitialised) Initialisation();
     153  if(SetupKinematics(part, mat, kinEnergy)) {
     154    if(part == electron || part == positron) {
     155      G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
     156      if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
     157      else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
     158      else res = ComputeRange(kinEnergy,part);
     159   
     160    } else {
     161      res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     162    }
     163  }
     164  return res;
     165}
     166
     167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     168
     169G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
     170                                                    const G4Material* mat,
     171                                                    G4double kinEnergy)
     172{
     173  if(!part || !mat || kinEnergy < keV) return false;
    147174  if(!isInitialised) Initialisation();
    148175  G4bool flag = false;
    … …  
    182209    if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
    183210  }
    184   G4double theta = ComputeScatteringAngle(stepLength);
    185   return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
     211  return true;
    186212}
    187213
    … …  
    231257  invRangeMuon     = PrepareTable();
    232258  invRangeProton   = PrepareTable();
     259  mscElectron      = PrepareTable();
    233260
    234261  G4LossTableBuilder builder;
    … …  
    262289  builder.BuildInverseRangeTable(rangeProton, invRangeProton); 
    263290
     291  ComputeTrasportXS(electron, mscElectron);
    264292}
    265293
    … …  
    273301
    274302    G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
     303    v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
    275304    table->push_back(v);
    276305  }
    … …  
    488517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    489518
     519void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
     520                                                    G4PhysicsTable* table)
     521{
     522  G4DataVector v;
     523  G4WentzelVIModel* msc = new G4WentzelVIModel();
     524  msc->SetPolarAngleLimit(CLHEP::pi);
     525  msc->Initialise(part, v);
     526
     527  mass    = part->GetPDGMass();
     528  charge2 = 1.0;
     529  currentParticle = part;
     530
     531  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
     532
     533  if(0<verbose) {
     534    G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
     535           << G4endl;
     536  }
     537 
     538  for(G4int i=0; i<nmat; i++) { 
     539
     540    const G4Material* mat = (*mtable)[i];
     541    if(1<verbose)
     542      G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
     543    G4PhysicsVector* aVector = (*table)[i];
     544    for(G4int j=0; j<nbins; j++) {
     545       
     546       G4double e = aVector->GetLowEdgeEnergy(j);
     547       G4double xs = msc->CrossSectionPerVolume(mat,part,e);
     548       aVector->PutValue(j,xs);
     549       if(1<verbose) {
     550         G4cout << "j= " << j << "  e(MeV)= " << e/MeV 
     551                << " xs(1/mm)= " << xs*mm << G4endl;
     552       }
     553    }
     554  }
     555  delete msc;
     556}
     557
     558//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     559

  • trunk/source/processes/electromagnetic/muons/src/G4MuBetheBlochModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBetheBlochModel.cc,v 1.23 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBetheBlochModel.cc,v 1.25 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8888  theElectron = G4Electron::Electron();
    8989  corr = G4LossTableManager::Instance()->EmCorrections();
     90  fParticleChange = 0;
    9091
    9192  if(p) SetParticle(p);
    … …  
    99100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    100101
    101 void G4MuBetheBlochModel::SetParticle(const G4ParticleDefinition* p)
    102 {
    103   if(!particle) {
    104     particle = p;
    105     mass = particle->GetPDGMass();
    106     massSquare = mass*mass;
    107     ratio = electron_mass_c2/mass;
    108   }
    109 }
    110 
    111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    112 
    113102G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
    114103                                           const G4MaterialCutsCouple* couple)
    … …  
    117106}
    118107
     108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     109
     110G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     111                                                 G4double kinEnergy)
     112{
     113  G4double tau  = kinEnergy/mass;
     114  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
     115                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
     116  return tmax;
     117}
     118
    119119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    120120
    … …  
    124124  if(p) SetParticle(p);
    125125
    126   if(pParticleChange)
    127     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
    128                                                              (pParticleChange);
    129   else
    130     fParticleChange = new G4ParticleChangeForLoss();
     126  if(!fParticleChange) {
     127    if(pParticleChange) {
     128      fParticleChange =
     129        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
     130    } else {
     131      fParticleChange = new G4ParticleChangeForLoss();
     132    }
     133  }
    131134}
    132135
    … …  
    275278
    276279  //High order corrections
    277   dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
     280  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
    278281
    279282  return dedx;

  • trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlung.cc,v 1.38 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlung.cc,v 1.42 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8282    lowestKinEnergy(1.*GeV),
    8383    isInitialised(false)
    84 {}
     84{
     85  SetProcessSubType(fBremsstrahlung);
     86}
    8587
    8688//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    9294
    93 void G4MuBremsstrahlung::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
    94                                                      const G4ParticleDefinition*)
     95G4bool G4MuBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
     96{
     97  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     98}
     99
     100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     101
     102G4double G4MuBremsstrahlung::MinPrimaryEnergy(const G4ParticleDefinition*,
     103                                              const G4Material*,
     104                                              G4double)
     105{
     106  return lowestKinEnergy;
     107}
     108
     109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     110
     111void G4MuBremsstrahlung::InitialiseEnergyLossProcess(
     112                                 const G4ParticleDefinition* part,
     113                                 const G4ParticleDefinition*)
    95114{
    96115  if(!isInitialised) {
    … …  
    105124    em->SetLowestKineticEnergy(lowestKinEnergy);
    106125
    107     G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
    108     em->SetLowEnergyLimit(0.1*keV);
    109     em->SetHighEnergyLimit(100.0*TeV);
     126    G4VEmFluctuationModel* fm = 0;
     127    em->SetLowEnergyLimit(MinKinEnergy());
     128    em->SetHighEnergyLimit(MaxKinEnergy());
    110129    AddEmModel(1, em, fm);
    111130  }
    … …  
    115134
    116135void G4MuBremsstrahlung::PrintInfo()
    117 {
    118   G4cout << "      Parametrised model "
    119          << G4endl;
    120 }
     136{}
    121137
    122138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlungModel.cc,v 1.24 2007/11/08 11:48:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlungModel.cc,v 1.33 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646// 13-02-03 Add name (V.Ivanchenko)
    4747// 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
    48 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
    49 // 03-08-05 Angular correlations according to PRM (V.Ivantchenko)
     48// 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
     49// 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
    5050// 13-02-06 add ComputeCrossSectionPerAtom (mma)
    5151// 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
    5252// 07-11-07 Improve sampling of final state (A.Bogdanov)
     53// 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
    5354//
    5455
    … …  
    7475
    7576//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    76 
    77 // static members
    78 //
    79 G4double G4MuBremsstrahlungModel::zdat[]={1., 4., 13., 29., 92.};
    80 G4double G4MuBremsstrahlungModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
    81 G4double G4MuBremsstrahlungModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
    82                                           1.e9, 1.e10};
    83 
    8477//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    8578
    … …  
    9083  : G4VEmModel(nam),
    9184    particle(0),
     85    sqrte(sqrt(exp(1.))),
     86    bh(202.4),
     87    bh1(446.),
     88    btf(183.),
     89    btf1(1429.),
     90    fParticleChange(0),
    9291    lowestKinEnergy(1.0*GeV),
    93     minThreshold(1.0*keV),
    94     nzdat(5),
    95     ntdat(8),
    96     NBIN(1000),
    97     cutFixed(0.98*keV),
    98     ignoreCut(false),
    99     samplingTablesAreFilled(false)
     92    minThreshold(1.0*keV)
    10093{
    10194  theGamma = G4Gamma::Gamma();
     95  nist = G4NistManager::Instance();
    10296  if(p) SetParticle(p);
    10397}
    … …  
    125119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    126120
    127 void G4MuBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
    128 {
    129   if(!particle) {
    130     particle = p;
    131     mass = particle->GetPDGMass();
    132   }
    133 }
    134 
    135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    136 
    137121void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
    138122                                         const G4DataVector& cuts)
    … …  
    142126  highKinEnergy = HighEnergyLimit();
    143127
     128  // partial cross section is computed for fixed energy
    144129  G4double fixedEnergy = 0.5*highKinEnergy;
    145130
    … …  
    148133  if(theCoupleTable) {
    149134    G4int numOfCouples = theCoupleTable->GetTableSize();
    150    
     135
     136    // clear old data   
    151137    G4int nn = partialSumSigma.size();
    152138    G4int nc = cuts.size();
    153139    if(nn > 0) {
    154140      for (G4int ii=0; ii<nn; ii++){
    155         G4DataVector* a=partialSumSigma[ii];
     141        G4DataVector* a = partialSumSigma[ii];
    156142        if ( a )  delete a;   
    157143      }
    158144      partialSumSigma.clear();
    159145    }
     146    // fill new data
    160147    if (numOfCouples>0) {
    161148      for (G4int i=0; i<numOfCouples; i++) {
    162149        G4double cute = DBL_MAX;
     150
     151        // protection for usage with extrapolator
    163152        if(i < nc) cute = cuts[i];
    164         if(cute < cutFixed || ignoreCut) cute = cutFixed;
     153
    165154        const G4MaterialCutsCouple* couple =
    166155          theCoupleTable->GetMaterialCutsCouple(i);
    … …  
    171160    }
    172161  }
    173   if(!samplingTablesAreFilled) MakeSamplingTables();
    174   if(pParticleChange)
    175     fParticleChange =
    176       reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
    177   else
    178     fParticleChange = new G4ParticleChangeForLoss();
     162
     163  // define pointer to G4ParticleChange
     164  if(!fParticleChange) {
     165    if(pParticleChange)
     166      fParticleChange =
     167        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
     168    else
     169      fParticleChange = new G4ParticleChangeForLoss();
     170  }
    179171}
    180172
    … …  
    188180{
    189181  G4double dedx = 0.0;
    190   if (kineticEnergy <= lowestKinEnergy || ignoreCut) return dedx;
     182  if (kineticEnergy <= lowestKinEnergy) return dedx;
    191183
    192184  G4double tmax = kineticEnergy;
    193   G4double cut  = min(cutEnergy,tmax);
    194   if(cut < cutFixed) cut = cutFixed;
     185  G4double cut  = std::min(cutEnergy,tmax);
     186  if(cut < minThreshold) cut = minThreshold;
    195187
    196188  const G4ElementVector* theElementVector = material->GetElementVector();
    197189  const G4double* theAtomicNumDensityVector =
    198                                           material->GetAtomicNumDensityVector();
     190    material->GetAtomicNumDensityVector();
    199191
    200192  //  loop for elements in the material
    201193  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    202194
    203     G4double Z = (*theElementVector)[i]->GetZ();
    204     G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
    205 
    206     G4double loss = ComputMuBremLoss(Z, A, kineticEnergy, cut);
     195    G4double loss =
     196      ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
    207197
    208198    dedx += loss*theAtomicNumDensityVector[i];
    209199  }
     200  //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
    210201  if(dedx < 0.) dedx = 0.;
    211202  return dedx;
    … …  
    214205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    215206
    216 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z, G4double A,
     207G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
    217208                                                   G4double tkin, G4double cut)
    218209{
    … …  
    239230    {
    240231      G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
    241       loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
     232      loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
    242233    }
    243234    aa += hhh;
    … …  
    254245                                           G4double tkin,
    255246                                           G4double Z,
    256                                            G4double A,
    257247                                           G4double cut)
    258248{
    … …  
    272262  G4double bbb = log(vmax);
    273263  G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
    274   G4double hhh = (bbb-aaa)/float(kkk);
     264  G4double hhh = (bbb-aaa)/G4double(kkk);
    275265
    276266  G4double aa = aaa;
    … …  
    281271    {
    282272      G4double ep = exp(aa + xgi[i]*hhh)*totalEnergy;
    283       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
     273      cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
    284274    }
    285275    aa += hhh;
    … …  
    287277
    288278  cross *=hhh;
     279
     280  //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
    289281
    290282  return cross;
    … …  
    296288                                           G4double tkin,
    297289                                           G4double Z,
    298                                            G4double A,
    299290                                           G4double gammaEnergy)
    300291//  differential cross section
    301292{
    302   static const G4double sqrte=sqrt(exp(1.)) ;
    303   static const G4double bh=202.4,bh1=446.,btf=183.,btf1=1429. ;
    304   static const G4double rmass=mass/electron_mass_c2 ;
    305   static const G4double cc=classic_electr_radius/rmass ;
    306   static const G4double coeff= 16.*fine_structure_const*cc*cc/3. ;
    307 
    308293  G4double dxsection = 0.;
    309294
    … …  
    315300  G4double rab0=delta*sqrte ;
    316301
    317   G4double z13 = exp(-log(Z)/3.) ;
    318   G4double dn  = 1.54*exp(0.27*log(A)) ;
     302  G4int iz = G4int(Z);
     303  if(iz < 1) iz = 1;
     304
     305  G4double z13 = 1.0/nist->GetZ13(iz);
     306  G4double dn  = 1.54*nist->GetA27(iz);
    319307
    320308  G4double b,b1,dnstar ;
    321309
    322   if(Z<1.5)
     310  if(1 == iz)
    323311  {
    324     b=bh;
    325     b1=bh1;
    326     dnstar=dn ;
     312    b  = bh;
     313    b1 = bh1;
     314    dnstar = dn;
    327315  }
    328316  else
    329317  {
    330     b=btf;
    331     b1=btf1;
    332     dnstar = exp((1.-1./Z)*log(dn)) ;
     318    b  = btf;
     319    b1 = btf1;
     320    dnstar = dn/std::pow(dn, 1./Z);
    333321  }
    334322
    … …  
    359347                                           const G4ParticleDefinition*,
    360348                                                 G4double kineticEnergy,
    361                                                  G4double Z, G4double A,
     349                                                 G4double Z, G4double,
    362350                                                 G4double cutEnergy,
    363                                                  G4double)
    364 {
    365   G4double cut  = min(cutEnergy, kineticEnergy);
    366   if(cut < cutFixed || ignoreCut) cut = cutFixed;
    367   G4double cross =
    368     ComputeMicroscopicCrossSection (kineticEnergy, Z, A/(g/mole), cut);
    369   return cross;
    370 }
    371 
    372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    373 
    374 G4double G4MuBremsstrahlungModel::CrossSectionPerVolume(
    375                                                const G4Material* material,
    376                                                const G4ParticleDefinition*,
    377                                                      G4double kineticEnergy,
    378                                                      G4double cutEnergy,
    379                                                      G4double maxEnergy)
     351                                                 G4double maxEnergy)
    380352{
    381353  G4double cross = 0.0;
    382   if (cutEnergy >= maxEnergy || kineticEnergy <= lowestKinEnergy) return cross;
    383  
    384   G4double tmax = min(maxEnergy, kineticEnergy);
    385   G4double cut  = min(cutEnergy, tmax);
    386   if(cut < cutFixed || ignoreCut) cut = cutFixed;
    387 
    388   const G4ElementVector* theElementVector = material->GetElementVector();
    389   const G4double* theAtomNumDensityVector =
    390                                           material->GetAtomicNumDensityVector();
    391 
    392   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    393 
    394     G4double Z = (*theElementVector)[i]->GetZ();
    395     G4double A = (*theElementVector)[i]->GetA()/(g/mole);
    396 
    397     G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
    398 
    399     if(tmax < kineticEnergy) {
    400       cr -= ComputeMicroscopicCrossSection(kineticEnergy, Z, A, tmax);
    401     }
    402     cross += theAtomNumDensityVector[i] * cr;
    403   }
    404 
     354  if (kineticEnergy <= lowestKinEnergy) return cross;
     355  G4double tmax = std::min(maxEnergy, kineticEnergy);
     356  G4double cut  = std::min(cutEnergy, kineticEnergy);
     357  if(cut < minThreshold) cut = minThreshold;
     358  if (cut >= tmax) return cross;
     359
     360  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
     361  if(tmax < kineticEnergy) {
     362    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
     363  }
    405364  return cross;
    406365}
    … …  
    410369G4DataVector* G4MuBremsstrahlungModel::ComputePartialSumSigma(
    411370                                       const G4Material* material,
    412                                              G4double kineticEnergy,
    413                                              G4double cut)
    414 
    415 // Build the table of cross section per element. The table is built for MATERIAL
    416 // This table is used by DoIt to select randomly an element in the material.
     371                                       G4double kineticEnergy,
     372                                       G4double cut)
     373
     374// Build the table of cross section per element.
     375// The table is built for material
     376// This table is used to select randomly an element in the material.
    417377{
    418378  G4int nElements = material->GetNumberOfElements();
    419379  const G4ElementVector* theElementVector = material->GetElementVector();
    420380  const G4double* theAtomNumDensityVector =
    421                                           material->GetAtomicNumDensityVector();
     381    material->GetAtomicNumDensityVector();
    422382
    423383  G4DataVector* dv = new G4DataVector();
    … …  
    426386
    427387  for (G4int i=0; i<nElements; i++ ) {
    428 
    429     G4double Z = (*theElementVector)[i]->GetZ();
    430     G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
    431388    cross += theAtomNumDensityVector[i]
    432              * ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
     389      * ComputeMicroscopicCrossSection(kineticEnergy,
     390                                       (*theElementVector)[i]->GetZ(), cut);
    433391    dv->push_back(cross);
    434392  }
    … …  
    438396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    439397
    440 void G4MuBremsstrahlungModel::MakeSamplingTables()
    441 {
    442 
    443   G4double AtomicNumber,AtomicWeight,KineticEnergy,
    444            TotalEnergy,Maxep;
    445 
    446   for (G4int iz=0; iz<nzdat; iz++)
    447    {
    448      AtomicNumber = zdat[iz];
    449      AtomicWeight = adat[iz]*g/mole ;
    450 
    451      for (G4int it=0; it<ntdat; it++)
    452      {
    453        KineticEnergy = tdat[it];
    454        TotalEnergy = KineticEnergy + mass;
    455        Maxep = KineticEnergy ;
    456 
    457        G4double CrossSection = 0.0 ;
    458 
    459        // calculate the differential cross section
    460        // numerical integration in
    461        //  log ...............
    462        G4double c = log(Maxep/cutFixed) ;
    463        G4double ymin = -5. ;
    464        G4double ymax = 0. ;
    465        G4double dy = (ymax-ymin)/NBIN ;
    466 
    467        G4double y = ymin - 0.5*dy ;
    468        G4double yy = ymin - dy ;
    469        G4double x = exp(y);
    470        G4double fac = exp(dy);
    471        G4double dx = exp(yy)*(fac - 1.0);
    472 
    473        for (G4int i=0 ; i<NBIN; i++)
    474        {
    475          y += dy ;
    476          x *= fac;
    477          dx*= fac;
    478          G4double ep = cutFixed*exp(c*x) ;
    479 
    480          CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
    481                                                  KineticEnergy,AtomicNumber,
    482                                                  AtomicWeight,ep) ;
    483          ya[i]=y ;
    484          proba[iz][it][i] = CrossSection ;
    485 
    486        }
    487 
    488        proba[iz][it][NBIN] = CrossSection ;
    489        ya[NBIN] = 0. ;   //   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    490 
    491        if(CrossSection > 0.)
    492        {
    493          for(G4int ib=0; ib<=NBIN; ib++)
    494          {
    495            proba[iz][it][ib] /= CrossSection ;
    496          }
    497        }
    498      }
    499    }
    500   samplingTablesAreFilled = true;
    501 }
    502 
    503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    504 
    505 void G4MuBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
    506                                                 const G4MaterialCutsCouple* couple,
    507                                                 const G4DynamicParticle* dp,
    508                                                 G4double minEnergy,
    509                                                 G4double maxEnergy)
     398void G4MuBremsstrahlungModel::SampleSecondaries(
     399                              std::vector<G4DynamicParticle*>* vdp,
     400                              const G4MaterialCutsCouple* couple,
     401                              const G4DynamicParticle* dp,
     402                              G4double minEnergy,
     403                              G4double maxEnergy)
    510404{
    511405  G4double kineticEnergy = dp->GetKineticEnergy();
    512406  // check against insufficient energy
    513   G4double tmax = min(kineticEnergy, maxEnergy);
    514   G4double tmin = min(kineticEnergy, minEnergy);
    515   if(tmin < cutFixed || ignoreCut) tmin = cutFixed;
     407  G4double tmax = std::min(kineticEnergy, maxEnergy);
     408  G4double tmin = std::min(kineticEnergy, minEnergy);
     409  if(tmin < minThreshold) tmin = minThreshold;
    516410  if(tmin >= tmax) return;
    517411
    518   // ===== the begining of a new code  ======
    519412  // ===== sampling of energy transfer ======
    520413
    … …  
    523416  // select randomly one element constituing the material
    524417  const G4Element* anElement = SelectRandomAtom(couple);
     418  G4double Z = anElement->GetZ();
    525419
    526420  G4double totalEnergy   = kineticEnergy + mass;
    527421  G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
    528422
    529   G4double AtomicNumber = anElement->GetZ();
    530   G4double AtomicWeight = anElement->GetA()/(g/mole);
    531 
    532   G4double func1 = tmin*ComputeDMicroscopicCrossSection(
    533                                     kineticEnergy,AtomicNumber,
    534                                     AtomicWeight,tmin);
     423  G4double func1 = tmin*
     424    ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
    535425
    536426  G4double lnepksi, epksi;
    537427  G4double func2;
    538   G4double ksi2;
    539428
    540429  do {
    541430    lnepksi = log(tmin) + G4UniformRand()*log(kineticEnergy/tmin);
    542431    epksi   = exp(lnepksi);
    543     func2   = epksi*ComputeDMicroscopicCrossSection(
    544                                 kineticEnergy,AtomicNumber,
    545                                 AtomicWeight,epksi);
    546     ksi2 = G4UniformRand();
    547 
    548   } while(func2/func1 < ksi2);
    549 
    550   // ===== the end of a new code =====
    551 
    552   // create G4DynamicParticle object for the Gamma
     432    func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
     433
     434  } while(func2 < func1*G4UniformRand());
     435
    553436  G4double gEnergy = epksi;
    554437
    555   // sample angle
     438  // ===== sample angle =====
     439
    556440  G4double gam  = totalEnergy/mass;
    557   G4double rmax = gam*min(1.0, totalEnergy/gEnergy - 1.0);
    558   rmax *= rmax;
    559   G4double x = G4UniformRand()*rmax/(1.0 + rmax);
     441  G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
     442  G4double rmax2= rmax*rmax;
     443  G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
    560444
    561445  G4double theta = sqrt(x/(1.0 - x))/gam;
    … …  
    577461
    578462  // save secondary
    579   G4DynamicParticle* aGamma = new G4DynamicParticle(theGamma,gDirection,gEnergy);
     463  G4DynamicParticle* aGamma =
     464    new G4DynamicParticle(theGamma,gDirection,gEnergy);
    580465  vdp->push_back(aGamma);
    581466}

  • trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuIonisation.cc,v 1.54 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuIonisation.cc,v 1.59 2009/02/26 11:04:20 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8686#include "G4MuBetheBlochModel.hh"
    8787#include "G4UniversalFluctuation.hh"
     88#include "G4IonFluctuations.hh"
    8889#include "G4BohrFluctuations.hh"
    8990#include "G4UnitsTable.hh"
    … …  
    99100    isInitialised(false)
    100101{
    101   SetStepFunction(0.2, 1*mm);
    102   SetIntegral(true);
    103   SetVerboseLevel(1);
     102  //  SetStepFunction(0.2, 1*mm);
     103  //SetIntegral(true);
     104  //SetVerboseLevel(1);
     105  SetProcessSubType(fIonisation);
    104106}
    105107
    … …  
    108110G4MuIonisation::~G4MuIonisation()
    109111{}
     112
     113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     114
     115G4bool G4MuIonisation::IsApplicable(const G4ParticleDefinition& p)
     116{
     117  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     118}
     119
     120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     121
     122G4double G4MuIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
     123                                          const G4Material*,
     124                                          G4double cut)
     125{
     126  G4double x = 0.5*cut/electron_mass_c2;
     127  G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
     128  return mass*(g - 1.0);
     129}
    110130
    111131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    122142    SetSecondaryParticle(G4Electron::Electron());
    123143
    124     flucModel = new G4UniversalFluctuation();
     144    // Bragg peak model
     145    if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
     146    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
     147    EmModel(1)->SetHighEnergyLimit(0.2*MeV);
     148    AddEmModel(1, EmModel(1), new G4IonFluctuations());
    125149
    126     G4VEmModel* em = new G4BraggModel();
    127     em->SetLowEnergyLimit(0.1*keV);
    128     em->SetHighEnergyLimit(0.2*MeV);
    129     AddEmModel(1, em, flucModel);
    130     G4VEmModel* em1 = new G4BetheBlochModel();
    131     em1->SetLowEnergyLimit(0.2*MeV);
    132     em1->SetHighEnergyLimit(1.0*GeV);
    133     AddEmModel(2, em1, flucModel);
    134     G4VEmModel* em2 = new G4MuBetheBlochModel();
    135     em2->SetLowEnergyLimit(1.0*GeV);
    136     em2->SetHighEnergyLimit(100.0*TeV);
    137     AddEmModel(3, em2, flucModel);
     150    // high energy fluctuation model
     151    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
     152
     153    // moderate energy model
     154    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);
     155    EmModel(2)->SetLowEnergyLimit(0.2*MeV);
     156    EmModel(2)->SetHighEnergyLimit(1.0*GeV);
     157    AddEmModel(2, EmModel(2), FluctModel());
     158
     159    // high energy model
     160    if (!EmModel(3)) SetEmModel(new G4MuBetheBlochModel(),3);
     161    EmModel(3)->SetLowEnergyLimit(1.0*GeV);
     162    EmModel(3)->SetHighEnergyLimit(MaxKinEnergy());
     163    AddEmModel(3, EmModel(3), FluctModel());
    138164
    139165    ratio = electron_mass_c2/mass;
    … …  
    145171
    146172void G4MuIonisation::PrintInfo()
    147 {
    148   G4cout << "      Bether-Bloch model for E > 0.2 MeV, "
    149          << "parametrisation of Bragg peak below, "
    150          << G4endl;
    151   G4cout << "      radiative corrections for E > 1 GeV" << G4endl;
    152 }
     173{}
    153174
    154175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuMultipleScattering.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuMultipleScattering.cc,v 1.3 2007/11/09 19:48:10 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuMultipleScattering.cc,v 1.12 2008/10/16 13:37:04 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -----------------------------------------------------------------------------
    … …  
    4747
    4848#include "G4MuMultipleScattering.hh"
    49 #include "G4MuMscModel.hh"
     49#include "G4WentzelVIModel.hh"
    5050#include "G4MscStepLimitType.hh"
    5151
    … …  
    6161  samplez           = false ;
    6262  isInitialized     = false; 
    63   SetRangeFactor(0.04);
     63  SetRangeFactor(0.2);
    6464  SetLateralDisplasmentFlag(true);
    6565}
    … …  
    8484  if(isInitialized) {
    8585
    86     if (p->GetParticleType() != "nucleus") {
     86    if (p->GetParticleType() != "nucleus" && p->GetPDGMass() < GeV) {
    8787      mscModel->SetStepLimitType(StepLimitType());
    8888      mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
    89       //mscModel->SetThetaLimit(thetaLimit);
    9089      mscModel->SetRangeFactor(RangeFactor());
    9190    }
     91    mscModel->SetPolarAngleLimit(PolarAngleLimit());
    9292    return;
    9393  }
    9494
    95   if (p->GetParticleType() == "nucleus") {
     95  if (p->GetParticleType() == "nucleus" || p->GetPDGMass() > GeV) {
    9696    SetLateralDisplasmentFlag(false);
    9797    SetBuildLambdaTable(false);
    98     //    SetRangeFactor(0.2);
    9998  }
    10099
    101   // initialisation of parameters
    102   //  G4String part_name = p->GetParticleName();
    103   mscModel = new G4MuMscModel(RangeFactor(),thetaLimit);
     100  // initialisation of the model
     101
     102  mscModel = new G4WentzelVIModel();
     103  mscModel->SetStepLimitType(StepLimitType());
    104104  mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
     105  mscModel->SetRangeFactor(RangeFactor());
     106  mscModel->SetPolarAngleLimit(PolarAngleLimit());
     107  mscModel->SetLowEnergyLimit(MinKinEnergy());
     108  mscModel->SetHighEnergyLimit(MaxKinEnergy());
    105109
    106110  AddEmModel(1,mscModel);
    … …  
    112116void G4MuMultipleScattering::PrintInfo()
    113117{
    114   G4cout << "      Boundary/stepping algorithm is active with RangeFactor= "
    115          << RangeFactor()
    116          << "  Step limit type " << StepLimitType()
    117         << G4endl;
     118  G4cout << "      RangeFactor= " << RangeFactor()
     119         << ", step limit type: " << StepLimitType()
     120         << ", lateralDisplacement: " << LateralDisplasmentFlag()
     121        << G4endl;
    118122}
    119123

  • trunk/source/processes/electromagnetic/muons/src/G4MuPairProduction.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProduction.cc,v 1.48 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProduction.cc,v 1.52 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8989    lowestKinEnergy(1.*GeV),
    9090    isInitialised(false)
    91 {}
     91{
     92  SetProcessSubType(fPairProdByCharged);
     93}
    9294
    9395//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    98100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    99101
    100 void G4MuPairProduction::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
    101                                                      const G4ParticleDefinition*)
     102G4bool G4MuPairProduction::IsApplicable(const G4ParticleDefinition& p)
     103{
     104  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     105}
     106
     107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     108
     109G4double G4MuPairProduction::MinPrimaryEnergy(const G4ParticleDefinition*,
     110                                              const G4Material*,
     111                                              G4double)
     112{
     113  return lowestKinEnergy;
     114}
     115
     116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     117
     118void G4MuPairProduction::InitialiseEnergyLossProcess(
     119                         const G4ParticleDefinition* part,
     120                         const G4ParticleDefinition*)
    102121{
    103122  if (!isInitialised) {
    … …  
    110129    G4MuPairProductionModel* em = new G4MuPairProductionModel();
    111130    em->SetLowestKineticEnergy(lowestKinEnergy);
    112     G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
    113     em->SetLowEnergyLimit(0.1*keV);
    114     em->SetHighEnergyLimit(100.0*TeV);
     131    G4VEmFluctuationModel* fm = 0;
     132    em->SetLowEnergyLimit(MinKinEnergy());
     133    em->SetHighEnergyLimit(MaxKinEnergy());
    115134    AddEmModel(1, em, fm);
    116135  }
    … …  
    120139
    121140void G4MuPairProduction::PrintInfo()
    122 {
    123   G4cout << "      Parametrised model "
    124          << G4endl;
    125 }
     141{}
    126142
    127143//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProductionModel.cc,v 1.35 2007/10/11 13:52:04 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProductionModel.cc,v 1.40 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    104104                                                 const G4String& nam)
    105105  : G4VEmModel(nam),
    106   minPairEnergy(4.*electron_mass_c2),
    107   lowestKinEnergy(1.*GeV),
    108   factorForCross(4.*fine_structure_const*fine_structure_const
     106    particle(0),
     107    factorForCross(4.*fine_structure_const*fine_structure_const
    109108                   *classic_electr_radius*classic_electr_radius/(3.*pi)),
    110109    sqrte(sqrt(exp(1.))),
    111110    currentZ(0),
    112     particle(0),
     111    fParticleChange(0),
     112    minPairEnergy(4.*electron_mass_c2),
     113    lowestKinEnergy(1.*GeV),
    113114    nzdat(5),
    114115    ntdat(8),
    … …  
    118119    ymax(0.),
    119120    dy((ymax-ymin)/nbiny),
    120     ignoreCut(false),
    121121    samplingTablesAreFilled(false)
    122122{
    123123  SetLowEnergyLimit(minPairEnergy);
     124  nist = G4NistManager::Instance();
    124125
    125126  theElectron = G4Electron::Electron();
    … …  
    144145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    145146
    146 void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
    147 {
    148   if(!particle) {
    149     particle = p;
    150     particleMass = particle->GetPDGMass();
    151   }
     147G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     148                                                     G4double kineticEnergy)
     149{
     150  G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
     151  return maxPairEnergy;
    152152}
    153153
    … …  
    161161    MakeSamplingTables();
    162162  }
    163   if(pParticleChange) {
    164     if(ignoreCut) {
    165       gParticleChange =
    166         reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    167       fParticleChange = 0;
    168     } else {
     163  if(!fParticleChange) {
     164    if(pParticleChange)
    169165      fParticleChange =
    170166        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
    171       gParticleChange = 0;
    172     }
    173   } else {
    174     fParticleChange = new G4ParticleChangeForLoss();
    175     gParticleChange = 0;
     167    else
     168      fParticleChange = new G4ParticleChangeForLoss();
    176169  }
    177170}
    … …  
    186179{
    187180  G4double dedx = 0.0;
    188   if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy
    189       || ignoreCut)
     181  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
    190182    return dedx;
    191183
    … …  
    216208  G4double loss = 0.0;
    217209
    218   G4double cut  = min(cutEnergy,tmax);
     210  G4double cut = std::min(cutEnergy,tmax);
    219211  if(cut <= minPairEnergy) return loss;
    220212
    … …  
    251243                                           G4double Z,
    252244                                           G4double cut)
    253 
    254 {
    255   G4double cross = 0. ;
    256 
     245{
     246  G4double cross = 0.;
    257247  SetCurrentElement(Z);
    258248  G4double tmax = MaxSecondaryEnergy(particle, tkin);
    259 
    260249  if (tmax <= cut) return cross;
    261250
    … …  
    400389                                                 G4double Z, G4double,
    401390                                                 G4double cutEnergy,
    402                                                  G4double)
    403 {
    404   G4double cut  = max(minPairEnergy,cutEnergy);
    405   if(ignoreCut) cut = minPairEnergy;
    406   G4double cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
    407   return cross;
    408 }
    409 
    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    411 
    412 G4double G4MuPairProductionModel::CrossSectionPerVolume(
    413                                                const G4Material* material,
    414                                                const G4ParticleDefinition*,
    415                                                      G4double kineticEnergy,
    416                                                      G4double cutEnergy,
    417                                                      G4double maxEnergy)
     391                                                 G4double maxEnergy)
    418392{
    419393  G4double cross = 0.0;
    420394  if (kineticEnergy <= lowestKinEnergy) return cross;
    421395
    422   maxEnergy += particleMass;
    423 
    424   const G4ElementVector* theElementVector = material->GetElementVector();
    425   const G4double* theAtomNumDensityVector = material->
    426                                                     GetAtomicNumDensityVector();
    427 
    428   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    429     G4double Z = (*theElementVector)[i]->GetZ();
    430     SetCurrentElement(Z);
    431     G4double tmax = min(maxEnergy,MaxSecondaryEnergy(particle, kineticEnergy));
    432     G4double cut  = max(minPairEnergy,cutEnergy);
    433     if(ignoreCut) cut = minPairEnergy;
    434     if(cut < tmax) {
    435       G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut)
    436                   - ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
    437 
    438       cross += theAtomNumDensityVector[i] * cr;
    439     }
     396  SetCurrentElement(Z);
     397  G4double tmax = std::min(maxEnergy, kineticEnergy);
     398  G4double cut  = std::min(cutEnergy, kineticEnergy);
     399  if(cut < minPairEnergy) cut = minPairEnergy;
     400  if (cut >= tmax) return cross;
     401
     402  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
     403  if(tmax < kineticEnergy) {
     404    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
    440405  }
    441406  return cross;
    … …  
    451416    SetCurrentElement(Z);
    452417
    453     for (G4int it=0; it<ntdat; it++)
    454     {
     418    for (G4int it=0; it<ntdat; it++) {
     419
    455420      G4double kineticEnergy = tdat[it];
    456421      G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
    457 
     422      // G4cout << "Z= " << currentZ << " z13= " << z13
     423      //<< " mE= " << maxPairEnergy << G4endl;
    458424      G4double CrossSection = 0.0 ;
    459425
    460       G4double y = ymin - 0.5*dy ;
    461       G4double yy = ymin - dy ;
    462       G4double x = exp(y);
    463       G4double fac = exp(dy);
    464       G4double dx = exp(yy)*(fac - 1.0);
    465 
    466       G4double c = log(maxPairEnergy/minPairEnergy);
    467 
    468       for (G4int i=0 ; i<nbiny; i++)
    469       {
    470         y += dy ;
    471         if(c > 0.0) {
    472           x *= fac;
    473           dx*= fac;
    474           G4double ep = minPairEnergy*exp(c*x) ;
    475           CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
    476                                       kineticEnergy, Z, ep);
    477         }
    478         ya[i] = y;
    479         proba[iz][it][i] = CrossSection;
     426      if(maxPairEnergy > minPairEnergy) {
     427
     428        G4double y = ymin - 0.5*dy ;
     429        G4double yy = ymin - dy ;
     430        G4double x = exp(y);
     431        G4double fac = exp(dy);
     432        G4double dx = exp(yy)*(fac - 1.0);
     433
     434        G4double c = log(maxPairEnergy/minPairEnergy);
     435
     436        for (G4int i=0 ; i<nbiny; i++) {
     437          y += dy ;
     438          if(c > 0.0) {
     439            x *= fac;
     440            dx*= fac;
     441            G4double ep = minPairEnergy*exp(c*x) ;
     442            CrossSection +=
     443              ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
     444          }
     445          ya[i] = y;
     446          proba[iz][it][i] = CrossSection;
     447        }
     448       
     449      } else {
     450        for (G4int i=0 ; i<nbiny; i++) {
     451          proba[iz][it][i] = CrossSection;
     452        }
    480453      }
    481454
    482455      ya[nbiny]=ymax;
    483 
    484456      proba[iz][it][nbiny] = CrossSection;
    485457
    … …  
    515487  G4double maxEnergy     = std::min(tmax, maxPairEnergy);
    516488  G4double minEnergy     = std::max(tmin, minPairEnergy);
    517   if(ignoreCut)minEnergy = minPairEnergy;
     489
    518490  if(minEnergy >= maxEnergy) return;
    519491  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
    … …  
    618590  // primary change
    619591  kineticEnergy -= (ElectronEnergy + PositronEnergy);
    620   if(fParticleChange)
    621     fParticleChange->SetProposedKineticEnergy(kineticEnergy);
    622   else
    623     gParticleChange->SetProposedKineticEnergy(kineticEnergy);
     592  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
    624593
    625594  vdp->push_back(aParticle1);
    … …  
    655624    G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
    656625    G4double minEnergy     = std::max(tmin, minPairEnergy);
    657     if(ignoreCut)minEnergy = minPairEnergy;
    658626
    659627    G4int iz;
Note: See TracChangeset for help on using the changeset viewer.