physics/src/G4RDeIonisationCrossSectionHandler.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // -------------------------------------------------------------------
0028 //
0029 // GEANT4 Class file
0030 //
0031 //
0032 // File name:     G4RDeIonisationCrossSectionHandler
0033 //
0034 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
0035 //
0036 // Creation date: 25 Sept 2001
0037 //
0038 // Modifications:
0039 // 10 Oct 2001  M.G. Pia     Revision to improve code quality and consistency with design
0040 // 19 Jul 2002   VI          Create composite data set for material
0041 // 21 Jan 2003  V.Ivanchenko Cut per region
0042 //
0043 // -------------------------------------------------------------------
0044
0045 #include "G4RDeIonisationCrossSectionHandler.hh"
0046 #include "G4PhysicalConstants.hh"
0047 #include "G4SystemOfUnits.hh"
0048 #include "G4RDVEnergySpectrum.hh"
0049 #include "G4DataVector.hh"
0050 #include "G4RDCompositeEMDataSet.hh"
0051 #include "G4RDVDataSetAlgorithm.hh"
0052 #include "G4RDLinLogLogInterpolation.hh"
0053 #include "G4RDSemiLogInterpolation.hh"
0054 #include "G4RDVEMDataSet.hh"
0055 #include "G4RDEMDataSet.hh"
0056 #include "G4Material.hh"
0057 #include "G4ProductionCutsTable.hh"
0058
0059
0060 G4RDeIonisationCrossSectionHandler::G4RDeIonisationCrossSectionHandler(
0061     const G4RDVEnergySpectrum* spec, G4RDVDataSetAlgorithm* alg,
0062     G4double emin, G4double emax, G4int nbin)
0063  :  G4RDVCrossSectionHandler(),
0064     theParam(spec)
0065 {
0066   G4RDVCrossSectionHandler::Initialise(alg, emin, emax, nbin);
0067   interp = new G4RDLinLogLogInterpolation();
0068 }
0069
0070
0071 G4RDeIonisationCrossSectionHandler::~G4RDeIonisationCrossSectionHandler()
0072 {
0073   delete interp;
0074 }
0075
0076
0077 std::vector<G4RDVEMDataSet*>* G4RDeIonisationCrossSectionHandler::BuildCrossSectionsForMaterials(
0078                         const G4DataVector& energyVector,
0079                         const G4DataVector* energyCuts)
0080 {
0081   G4int verbose = 0;
0082   std::vector<G4RDVEMDataSet*>* set = new std::vector<G4RDVEMDataSet*>;
0083
0084   G4DataVector* energies;
0085   G4DataVector* cs;
0086   G4int nOfBins = energyVector.size();
0087
0088   const G4ProductionCutsTable* theCoupleTable=
0089         G4ProductionCutsTable::GetProductionCutsTable();
0090   size_t numOfCouples = theCoupleTable->GetTableSize();
0091
0092   for (size_t m=0; m<numOfCouples; m++) {
0093
0094     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
0095     const G4Material* material= couple->GetMaterial();
0096     const G4ElementVector* elementVector = material->GetElementVector();
0097     const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
0098     G4int nElements = material->GetNumberOfElements();
0099
0100     if(verbose > 0) {
0101       G4cout << "eIonisation CS for " << m << "th material "
0102              << material->GetName()
0103              << "  eEl= " << nElements << G4endl;
0104     }
0105
0106     G4double tcut  = (*energyCuts)[m];
0107
0108     G4RDVDataSetAlgorithm* algo = interp->Clone();
0109     G4RDVEMDataSet* setForMat = new G4RDCompositeEMDataSet(algo,1.,1.);
0110
0111     for (G4int i=0; i<nElements; i++) {
0112
0113       G4int Z = (G4int) (*elementVector)[i]->GetZ();
0114       G4int nShells = NumberOfComponents(Z);
0115       energies = new G4DataVector;
0116       cs       = new G4DataVector;
0117       G4double density = nAtomsPerVolume[i];
0118
0119       for (G4int bin=0; bin<nOfBins; bin++) {
0120
0121         G4double e = energyVector[bin];
0122         energies->push_back(e);
0123         G4double value = 0.0;
0124
0125         if(e > tcut) {
0126           for (G4int n=0; n<nShells; n++) {
0127             G4double cross = FindValue(Z, e, n);
0128             G4double p = theParam->Probability(Z, tcut, e, e, n);
0129             value += cross * p * density;
0130
0131             if(verbose>0 && m == 0 && e>=1. && e<=0.) {
0132               G4cout << "G4eIonCrossSH: e(MeV)= " << e/MeV
0133                      << " n= " << n
0134                      << " cross= " << cross
0135                      << " p= " << p
0136                      << " value= " << value
0137                      << " tcut(MeV)= " << tcut/MeV
0138                      << " rho= " << density
0139                      << " Z= " << Z
0140                      << G4endl;
0141         }
0142
0143       }
0144     }
0145         cs->push_back(value);
0146       }
0147       G4RDVDataSetAlgorithm* algo = interp->Clone();
0148       G4RDVEMDataSet* elSet = new G4RDEMDataSet(i,energies,cs,algo,1.,1.);
0149       setForMat->AddComponent(elSet);
0150     }
0151     set->push_back(setForMat);
0152   }
0153
0154   return set;
0155 }
0156
0157