G4EnergyLossForExtrapolator.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 //---------------------------------------------------------------------------
00029 //
00030 // ClassName:    G4EnergyLossForExtrapolator
00031 //  
00032 // Description:  This class provide calculation of energy loss, fluctuation, 
00033 //               and msc angle
00034 //
00035 // Author:       09.12.04 V.Ivanchenko 
00036 //
00037 // Modification: 
00038 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
00039 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
00040 // 21-03-06 Add verbosity defined in the constructor and Initialisation
00041 //          start only when first public method is called (V.Ivanchenko)
00042 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
00043 // 12-05-06 SEt linLossLimit=0.001 (VI)
00044 //
00045 //----------------------------------------------------------------------------
00046 //
00047  
00048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00049 
00050 #include "G4EnergyLossForExtrapolator.hh"
00051 #include "G4PhysicalConstants.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4PhysicsLogVector.hh"
00054 #include "G4ParticleDefinition.hh"
00055 #include "G4Material.hh"
00056 #include "G4MaterialCutsCouple.hh"
00057 #include "G4Electron.hh"
00058 #include "G4Positron.hh"
00059 #include "G4Proton.hh"
00060 #include "G4MuonPlus.hh"
00061 #include "G4MuonMinus.hh"
00062 #include "G4ParticleTable.hh"
00063 #include "G4LossTableBuilder.hh"
00064 #include "G4MollerBhabhaModel.hh"
00065 #include "G4BetheBlochModel.hh"
00066 #include "G4eBremsstrahlungModel.hh"
00067 #include "G4MuPairProductionModel.hh"
00068 #include "G4MuBremsstrahlungModel.hh"
00069 #include "G4ProductionCuts.hh"
00070 #include "G4LossTableManager.hh"
00071 #include "G4WentzelVIModel.hh"
00072 
00073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00074 
00075 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
00076   :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false)
00077 {
00078   currentParticle = 0;
00079   currentMaterial = 0;
00080 
00081   linLossLimit = 0.01;
00082   emin         = 1.*MeV;
00083   emax         = 10.*TeV;
00084   nbins        = 70;
00085 
00086   nmat = index = 0;
00087   cuts = 0;
00088 
00089   mass = charge2 = electronDensity = radLength = bg2 = beta2 
00090     = kineticEnergy = tmax = 0;
00091   gam = 1.0;
00092 
00093   dedxElectron = dedxPositron = dedxProton = rangeElectron 
00094     = rangePositron = rangeProton = invRangeElectron = invRangePositron 
00095     = invRangeProton = mscElectron = dedxMuon = rangeMuon = invRangeMuon = 0;
00096   cuts = 0;
00097   electron = positron = proton = muonPlus = muonMinus = 0;
00098 }
00099 
00100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00101 
00102 G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator()
00103 {
00104   for(G4int i=0; i<nmat; i++) {delete couples[i];}
00105   delete dedxElectron;
00106   delete dedxPositron;
00107   delete dedxProton;
00108   delete dedxMuon;
00109   delete rangeElectron;
00110   delete rangePositron;
00111   delete rangeProton;
00112   delete rangeMuon;
00113   delete invRangeElectron;
00114   delete invRangePositron;
00115   delete invRangeProton;
00116   delete invRangeMuon;
00117   delete mscElectron;
00118   delete cuts;
00119 }
00120 
00121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00122 
00123 G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy, 
00124                                                       G4double stepLength, 
00125                                                       const G4Material* mat, 
00126                                                       const G4ParticleDefinition* part)
00127 {
00128   if(!isInitialised) Initialisation();
00129   G4double kinEnergyFinal = kinEnergy;
00130   if(SetupKinematics(part, mat, kinEnergy)) {
00131     G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
00132     G4double r  = ComputeRange(kinEnergy,part);
00133     if(r <= step) {
00134       kinEnergyFinal = 0.0;
00135     } else if(step < linLossLimit*r) {
00136       kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
00137     } else {  
00138       G4double r1 = r - step;
00139       kinEnergyFinal = ComputeEnergy(r1,part);
00140     }
00141   }
00142   return kinEnergyFinal;
00143 }
00144 
00145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00146 
00147 G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy, 
00148                                                        G4double stepLength, 
00149                                                        const G4Material* mat, 
00150                                                        const G4ParticleDefinition* part)
00151 {
00152   if(!isInitialised) Initialisation();
00153   G4double kinEnergyFinal = kinEnergy;
00154 
00155   if(SetupKinematics(part, mat, kinEnergy)) {
00156     G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
00157     G4double r  = ComputeRange(kinEnergy,part);
00158 
00159     if(step < linLossLimit*r) {
00160       kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
00161     } else {  
00162       G4double r1 = r + step;
00163       kinEnergyFinal = ComputeEnergy(r1,part);
00164     }
00165   }
00166   return kinEnergyFinal;
00167 }
00168 
00169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00170 
00171 G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy, 
00172                                                      G4double stepLength,
00173                                                      const G4Material* mat, 
00174                                                      const G4ParticleDefinition* part)
00175 {
00176   G4double res = stepLength;
00177   if(!isInitialised) Initialisation();
00178   if(SetupKinematics(part, mat, kinEnergy)) {
00179     if(part == electron || part == positron) {
00180       G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
00181       if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
00182       else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
00183       else res = ComputeRange(kinEnergy,part);
00184     
00185     } else {
00186       res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
00187     }
00188   }
00189   return res;
00190 }
00191 
00192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00193 
00194 G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part, 
00195                                                     const G4Material* mat, 
00196                                                     G4double kinEnergy)
00197 {
00198   if(!part || !mat || kinEnergy < keV) return false;
00199   if(!isInitialised) Initialisation();
00200   G4bool flag = false;
00201   if(part != currentParticle) {
00202     flag = true;
00203     currentParticle = part;
00204     mass = part->GetPDGMass();
00205     G4double q = part->GetPDGCharge()/eplus;
00206     charge2 = q*q;
00207   }
00208   if(mat != currentMaterial) {
00209     G4int i = mat->GetIndex();
00210     if(i >= nmat) {
00211       G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= " 
00212              << i << " is out of table - NO extrapolation" << G4endl;
00213     } else {
00214       flag = true;
00215       currentMaterial = mat;
00216       electronDensity = mat->GetElectronDensity();
00217       radLength       = mat->GetRadlen();
00218       index           = i;
00219     }
00220   }
00221   if(flag || kinEnergy != kineticEnergy) {
00222     kineticEnergy = kinEnergy;
00223     G4double tau  = kinEnergy/mass;
00224 
00225     gam   = tau + 1.0;
00226     bg2   = tau * (tau + 2.0);
00227     beta2 = bg2/(gam*gam);
00228     tmax  = kinEnergy;
00229     if(part == electron) tmax *= 0.5;
00230     else if(part != positron) {
00231       G4double r = electron_mass_c2/mass;
00232       tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
00233     }
00234     if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
00235   }
00236   return true;
00237 } 
00238 
00239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00240 
00241 void G4EnergyLossForExtrapolator::Initialisation()
00242 {
00243   isInitialised = true;
00244   if(verbose>1) {
00245     G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
00246   }
00247   currentParticle = 0;
00248   currentMaterial = 0;
00249   kineticEnergy   = 0.0;
00250 
00251   electron = G4Electron::Electron();
00252   positron = G4Positron::Positron();
00253   proton   = G4Proton::Proton();
00254   muonPlus = G4MuonPlus::MuonPlus();
00255   muonMinus= G4MuonMinus::MuonMinus();
00256 
00257   currentParticleName = "";
00258 
00259   nmat = G4Material::GetNumberOfMaterials();
00260   const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00261   cuts = new G4ProductionCuts();
00262 
00263   const G4MaterialCutsCouple* couple;
00264   for(G4int i=0; i<nmat; i++) {
00265     couple = new G4MaterialCutsCouple((*mtable)[i],cuts);  
00266     couples.push_back(couple);
00267   }
00268 
00269   dedxElectron     = PrepareTable();
00270   dedxPositron     = PrepareTable();
00271   dedxMuon         = PrepareTable();
00272   dedxProton       = PrepareTable();
00273   rangeElectron    = PrepareTable();
00274   rangePositron    = PrepareTable();
00275   rangeMuon        = PrepareTable();
00276   rangeProton      = PrepareTable();
00277   invRangeElectron = PrepareTable();
00278   invRangePositron = PrepareTable();
00279   invRangeMuon     = PrepareTable();
00280   invRangeProton   = PrepareTable();
00281   mscElectron      = PrepareTable();
00282 
00283   G4LossTableBuilder builder; 
00284 
00285   if(verbose>1) 
00286     G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl;
00287 
00288   ComputeElectronDEDX(electron, dedxElectron);
00289   builder.BuildRangeTable(dedxElectron,rangeElectron);  
00290   builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);  
00291 
00292   if(verbose>1) 
00293     G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl;
00294 
00295   ComputeElectronDEDX(positron, dedxPositron);
00296   builder.BuildRangeTable(dedxPositron, rangePositron);  
00297   builder.BuildInverseRangeTable(rangePositron, invRangePositron);  
00298 
00299   if(verbose>1) 
00300     G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl;
00301 
00302   ComputeMuonDEDX(muonPlus, dedxMuon);
00303   builder.BuildRangeTable(dedxMuon, rangeMuon);  
00304   builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);  
00305 
00306   if(verbose>1) 
00307     G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl;
00308 
00309   ComputeProtonDEDX(proton, dedxProton);
00310   builder.BuildRangeTable(dedxProton, rangeProton);  
00311   builder.BuildInverseRangeTable(rangeProton, invRangeProton);  
00312 
00313   ComputeTrasportXS(electron, mscElectron);
00314 }
00315 
00316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00317 
00318 G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable()
00319 {
00320   G4PhysicsTable* table = new G4PhysicsTable();
00321 
00322   for(G4int i=0; i<nmat; i++) {  
00323 
00324     G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
00325     v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
00326     table->push_back(v);
00327   }
00328   return table;
00329 }
00330 
00331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00332 
00333 const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
00334 {
00335   const G4ParticleDefinition* p = 0;
00336   if(name != currentParticleName) {
00337     p = G4ParticleTable::GetParticleTable()->FindParticle(name);
00338     if(!p) {
00339       G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find " 
00340              << name << G4endl;
00341     }
00342   } else {
00343     p = currentParticle;
00344   }
00345   return p;
00346 }
00347 
00348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00349 
00350 G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy, 
00351                                                   const G4ParticleDefinition* part)
00352 {
00353   G4double x = 0.0;
00354   if(part == electron)      x = ComputeValue(kinEnergy, dedxElectron);
00355   else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
00356   else if(part == muonPlus || part == muonMinus) {
00357     x = ComputeValue(kinEnergy, dedxMuon);
00358   } else {
00359     G4double e = kinEnergy*proton_mass_c2/mass;
00360     x = ComputeValue(e, dedxProton)*charge2;
00361   }
00362   return x;
00363 }
00364   
00365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00366 
00367 G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy, 
00368                                                    const G4ParticleDefinition* part)
00369 {
00370   G4double x = 0.0;
00371   if(part == electron)      x = ComputeValue(kinEnergy, rangeElectron);
00372   else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
00373   else if(part == muonPlus || part == muonMinus) 
00374     x = ComputeValue(kinEnergy, rangeMuon);
00375   else {
00376     G4double massratio = proton_mass_c2/mass;
00377     G4double e = kinEnergy*massratio;
00378     x = ComputeValue(e, rangeProton)/(charge2*massratio);
00379   }
00380   return x;
00381 }
00382 
00383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00384 
00385 G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range, 
00386                                                     const G4ParticleDefinition* part)
00387 {
00388   G4double x = 0.0;
00389   if(part == electron)      x = ComputeValue(range, invRangeElectron);
00390   else if(part == positron) x = ComputeValue(range, invRangePositron);
00391   else if(part == muonPlus || part == muonMinus) 
00392     x = ComputeValue(range, invRangeMuon);
00393   else {
00394     G4double massratio = proton_mass_c2/mass;
00395     G4double r = range*massratio*charge2;
00396     x = ComputeValue(r, invRangeProton)/massratio;
00397   }
00398   return x;
00399 }
00400 
00401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00402 
00403 void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part, 
00404                                                       G4PhysicsTable* table) 
00405 {
00406   G4DataVector v;
00407   G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
00408   G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
00409   ioni->Initialise(part, v);
00410   brem->Initialise(part, v);
00411 
00412   mass    = electron_mass_c2;
00413   charge2 = 1.0;
00414   currentParticle = part;
00415 
00416   const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00417   if(0<verbose) {
00418     G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for " 
00419            << part->GetParticleName() 
00420            << G4endl;
00421   }
00422   for(G4int i=0; i<nmat; i++) {  
00423 
00424     const G4Material* mat = (*mtable)[i];
00425     if(1<verbose)
00426       G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
00427     const G4MaterialCutsCouple* couple = couples[i];
00428     G4PhysicsVector* aVector = (*table)[i];
00429 
00430     for(G4int j=0; j<=nbins; j++) {
00431         
00432        G4double e = aVector->Energy(j);
00433        G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e);
00434        if(1<verbose) {
00435          G4cout << "j= " << j
00436                 << "  e(MeV)= " << e/MeV  
00437                 << " dedx(Mev/cm)= " << dedx*cm/MeV
00438                 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
00439        }
00440        aVector->PutValue(j,dedx);
00441     }
00442   }
00443   delete ioni;
00444   delete brem;
00445 } 
00446 
00447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00448 
00449 void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part, 
00450                                                   G4PhysicsTable* table)
00451 {
00452   G4DataVector v;
00453   G4BetheBlochModel* ioni = new G4BetheBlochModel();
00454   G4MuPairProductionModel* pair = new G4MuPairProductionModel();
00455   G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
00456   ioni->Initialise(part, v);
00457   pair->Initialise(part, v);
00458   brem->Initialise(part, v);
00459 
00460   mass    = part->GetPDGMass();
00461   charge2 = 1.0;
00462   currentParticle = part;
00463 
00464   const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00465 
00466   if(0<verbose) {
00467     G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName() 
00468            << G4endl;
00469   }
00470  
00471   for(G4int i=0; i<nmat; i++) {  
00472 
00473     const G4Material* mat = (*mtable)[i];
00474     if(1<verbose)
00475       G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
00476     const G4MaterialCutsCouple* couple = couples[i];
00477     G4PhysicsVector* aVector = (*table)[i];
00478     for(G4int j=0; j<=nbins; j++) {
00479         
00480        G4double e = aVector->Energy(j);
00481        G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
00482                        pair->ComputeDEDX(couple,part,e,e) +
00483                        brem->ComputeDEDX(couple,part,e,e);
00484        aVector->PutValue(j,dedx);
00485        if(1<verbose) {
00486          G4cout << "j= " << j
00487                 << "  e(MeV)= " << e/MeV  
00488                 << " dedx(Mev/cm)= " << dedx*cm/MeV
00489                 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
00490        }
00491     }
00492   }
00493   delete ioni;
00494 }
00495 
00496 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00497 
00498 void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part, 
00499                                                     G4PhysicsTable* table)
00500 {
00501   G4DataVector v;
00502   G4BetheBlochModel* ioni = new G4BetheBlochModel();
00503   ioni->Initialise(part, v);
00504 
00505   mass    = part->GetPDGMass();
00506   charge2 = 1.0;
00507   currentParticle = part;
00508 
00509   const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00510 
00511   if(0<verbose) {
00512     G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() 
00513            << G4endl;
00514   }
00515  
00516   for(G4int i=0; i<nmat; i++) {  
00517 
00518     const G4Material* mat = (*mtable)[i];
00519     if(1<verbose)
00520       G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
00521     const G4MaterialCutsCouple* couple = couples[i];
00522     G4PhysicsVector* aVector = (*table)[i];
00523     for(G4int j=0; j<=nbins; j++) {
00524         
00525        G4double e = aVector->Energy(j);
00526        G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
00527        aVector->PutValue(j,dedx);
00528        if(1<verbose) {
00529          G4cout << "j= " << j
00530                 << "  e(MeV)= " << e/MeV  
00531                 << " dedx(Mev/cm)= " << dedx*cm/MeV
00532                 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl;
00533        }
00534     }
00535   }
00536   delete ioni;
00537 }
00538 
00539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00540 
00541 void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part, 
00542                                                     G4PhysicsTable* table)
00543 {
00544   G4DataVector v;
00545   G4WentzelVIModel* msc = new G4WentzelVIModel();
00546   msc->SetPolarAngleLimit(CLHEP::pi);
00547   msc->Initialise(part, v);
00548 
00549   mass    = part->GetPDGMass();
00550   charge2 = 1.0;
00551   currentParticle = part;
00552 
00553   const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00554 
00555   if(0<verbose) {
00556     G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() 
00557            << G4endl;
00558   }
00559  
00560   for(G4int i=0; i<nmat; i++) {  
00561 
00562     const G4Material* mat = (*mtable)[i];
00563     msc->SetCurrentCouple(couples[i]);
00564     if(1<verbose)
00565       G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
00566     G4PhysicsVector* aVector = (*table)[i];
00567     for(G4int j=0; j<=nbins; j++) {
00568         
00569        G4double e = aVector->Energy(j);
00570        G4double xs = msc->CrossSectionPerVolume(mat,part,e);
00571        aVector->PutValue(j,xs);
00572        if(1<verbose) {
00573          G4cout << "j= " << j << "  e(MeV)= " << e/MeV  
00574                 << " xs(1/mm)= " << xs*mm << G4endl;
00575        }
00576     }
00577   }
00578   delete msc;
00579 }
00580 
00581 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00582 

Generated on Mon May 27 17:48:09 2013 for Geant4 by  doxygen 1.4.7