G4eIonisationCrossSectionHandler.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 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4eIonisationCrossSectionHandler
00034 //
00035 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
00036 // 
00037 // Creation date: 25 Sept 2001
00038 //
00039 // Modifications: 
00040 // 10 Oct 2001  M.G. Pia     Revision to improve code quality and consistency with design
00041 // 19 Jul 2002   VI          Create composite data set for material
00042 // 21 Jan 2003  V.Ivanchenko Cut per region
00043 // 28 Jan 2009  L.Pandola    Added public method to make a easier migration of 
00044 //                           G4LowEnergyIonisation to G4LivermoreIonisationModel
00045 // 15 Jul 2009   Nicolas A. Karakatsanis
00046 //
00047 //                           - BuildCrossSectionForMaterials method was revised in order to calculate the 
00048 //                             logarithmic values of the loaded data. 
00049 //                             It retrieves the data values from the G4EMLOW data files but, then, calculates the
00050 //                             respective log values and loads them to seperate data structures.
00051 //                             The EM data sets, initialized this way, contain both non-log and log values.
00052 //                             These initialized data sets can enhance the computing performance of data interpolation
00053 //                             operations
00054 //
00055 //
00056 //
00057 // -------------------------------------------------------------------
00058 
00059 #include "G4eIonisationCrossSectionHandler.hh"
00060 #include "G4SystemOfUnits.hh"
00061 #include "G4VEnergySpectrum.hh"
00062 #include "G4DataVector.hh"
00063 #include "G4CompositeEMDataSet.hh"
00064 #include "G4VDataSetAlgorithm.hh"
00065 #include "G4LinLogLogInterpolation.hh"
00066 #include "G4SemiLogInterpolation.hh"
00067 #include "G4VEMDataSet.hh"
00068 #include "G4EMDataSet.hh"
00069 #include "G4Material.hh"
00070 #include "G4ProductionCutsTable.hh"
00071 
00072 
00073 G4eIonisationCrossSectionHandler::G4eIonisationCrossSectionHandler(
00074     const G4VEnergySpectrum* spec, G4VDataSetAlgorithm* alg,
00075     G4double emin, G4double emax, G4int nbin)
00076  :  G4VCrossSectionHandler(),
00077     theParam(spec),verbose(0)
00078 {
00079   G4VCrossSectionHandler::Initialise(alg, emin, emax, nbin);
00080   interp = new G4LinLogLogInterpolation();
00081 }
00082 
00083 
00084 G4eIonisationCrossSectionHandler::~G4eIonisationCrossSectionHandler()
00085 {
00086   delete interp;
00087 }
00088 
00089 
00090 std::vector<G4VEMDataSet*>* G4eIonisationCrossSectionHandler::BuildCrossSectionsForMaterials(
00091                         const G4DataVector& energyVector,
00092                         const G4DataVector* energyCuts)
00093 {
00094   std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
00095 
00096   G4DataVector* energies;
00097   G4DataVector* cs;
00098 
00099   G4DataVector* log_energies;
00100   G4DataVector* log_cs;
00101 
00102   G4int nOfBins = energyVector.size();
00103 
00104   const G4ProductionCutsTable* theCoupleTable=
00105         G4ProductionCutsTable::GetProductionCutsTable();
00106   size_t numOfCouples = theCoupleTable->GetTableSize();
00107 
00108   for (size_t mLocal=0; mLocal<numOfCouples; mLocal++) {
00109 
00110     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(mLocal);
00111     const G4Material* material= couple->GetMaterial();
00112     const G4ElementVector* elementVector = material->GetElementVector();
00113     const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
00114     G4int nElements = material->GetNumberOfElements();
00115 
00116     if(verbose > 0) 
00117       {
00118         G4cout << "eIonisation CS for " << mLocal << "th material "
00119                << material->GetName()
00120                << "  eEl= " << nElements << G4endl;
00121       }
00122 
00123     G4double tcut  = (*energyCuts)[mLocal];
00124 
00125     G4VDataSetAlgorithm* algo = interp->Clone();
00126     G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
00127 
00128     for (G4int i=0; i<nElements; i++) {
00129 
00130       G4int Z = (G4int) (*elementVector)[i]->GetZ();
00131       G4int nShells = NumberOfComponents(Z);
00132 
00133       energies = new G4DataVector;
00134       cs       = new G4DataVector;
00135 
00136       log_energies = new G4DataVector;
00137       log_cs       = new G4DataVector;
00138 
00139       G4double density = nAtomsPerVolume[i];
00140 
00141       for (G4int bin=0; bin<nOfBins; bin++) {
00142 
00143         G4double e = energyVector[bin];
00144         energies->push_back(e);
00145         log_energies->push_back(std::log10(e));
00146         G4double value = 0.0;
00147         G4double log_value = -300;
00148 
00149         if(e > tcut) {
00150           for (G4int n=0; n<nShells; n++) {
00151             G4double cross = FindValue(Z, e, n);
00152             G4double p = theParam->Probability(Z, tcut, e, e, n);
00153             value += cross * p * density;
00154 
00155             if(verbose>0 && mLocal == 0 && e>=1. && e<=0.) 
00156             {
00157               G4cout << "G4eIonCrossSH: e(MeV)= " << e/MeV
00158                      << " n= " << n
00159                      << " cross= " << cross
00160                      << " p= " << p
00161                      << " value= " << value
00162                      << " tcut(MeV)= " << tcut/MeV
00163                      << " rho= " << density
00164                      << " Z= " << Z
00165                      << G4endl;
00166               }
00167 
00168           }
00169           if (value == 0.) value = 1e-300;
00170           log_value = std::log10(value);
00171         }
00172         cs->push_back(value);
00173         log_cs->push_back(log_value);
00174       }
00175       G4VDataSetAlgorithm* algoLocal = interp->Clone();
00176 
00177       //G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,algoLocal,1.,1.);
00178 
00179       G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,log_energies,log_cs,algoLocal,1.,1.);
00180 
00181       setForMat->AddComponent(elSet);
00182     }
00183     set->push_back(setForMat);
00184   }
00185 
00186   return set;
00187 }
00188 
00189 G4double G4eIonisationCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
00190                                                                                    G4double cutEnergy,
00191                                                                                    G4int Z)
00192 {
00193   G4int nShells = NumberOfComponents(Z);
00194   G4double value = 0.;
00195   if(energy > cutEnergy) 
00196     {
00197       for (G4int n=0; n<nShells; n++) {
00198         G4double cross = FindValue(Z, energy, n);
00199         G4double p = theParam->Probability(Z, cutEnergy, energy, energy, n);
00200         value += cross * p;
00201       }
00202     }
00203   return value;
00204 }

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