physics/src/G4RDBremsstrahlungParameters.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 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0029 //         V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
0030 //
0031 // History:
0032 // -----------
0033 // 31 Jul 2001   MGP        Created
0034 // 12.09.01 V.Ivanchenko    Add activeZ and paramA
0035 // 25.09.01 V.Ivanchenko    Add parameter C and change interface to B
0036 // 29.11.01 V.Ivanchenko    Update parametrisation
0037 // 18.11.02 V.Ivanchenko    Fix problem of load
0038 // 21.02.03 V.Ivanchenko    Number of parameters is defined in the constructor
0039 // 28.02.03 V.Ivanchenko    Filename is defined in the constructor
0040 //
0041 // -------------------------------------------------------------------
0042
0043 #include "G4RDBremsstrahlungParameters.hh"
0044 #include "G4RDVEMDataSet.hh"
0045 #include "G4RDEMDataSet.hh"
0046 #include "G4RDLogLogInterpolation.hh"
0047 #include "G4Material.hh"
0048 #include <fstream>
0049
0050
0051 G4RDBremsstrahlungParameters:: G4RDBremsstrahlungParameters(const G4String& name,
0052     size_t num, G4int minZ, G4int maxZ)
0053   : zMin(minZ),
0054     zMax(maxZ),
0055     length(num)
0056 {
0057   LoadData(name);
0058 }
0059
0060
0061 G4RDBremsstrahlungParameters::~G4RDBremsstrahlungParameters()
0062 {
0063   // Reset the map of data sets: remove the data sets from the map
0064   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
0065
0066   for (pos = param.begin(); pos != param.end(); ++pos)
0067     {
0068       G4RDVEMDataSet* dataSet = (*pos).second;
0069       delete dataSet;
0070     }
0071
0072   activeZ.clear();
0073   paramC.clear();
0074 }
0075
0076
0077 G4double G4RDBremsstrahlungParameters::Parameter(G4int parameterIndex,
0078                                                  G4int Z,
0079                                                  G4double energy) const
0080 {
0081   G4double value = 0.;
0082   G4int id = Z*length + parameterIndex;
0083   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0084
0085   pos = param.find(id);
0086   if (pos!= param.end()) {
0087
0088     G4RDVEMDataSet* dataSet = (*pos).second;
0089     const G4DataVector ener = dataSet->GetEnergies(0);
0090     G4double ee = std::max(ener.front(),std::min(ener.back(),energy));
0091     value = dataSet->FindValue(ee);
0092
0093   } else {
0094     G4cout << "WARNING: G4RDBremsstrahlungParameters::FindValue "
0095            << "did not find ID = "
0096            << id << G4endl;
0097   }
0098
0099   return value;
0100 }
0101
0102 void G4RDBremsstrahlungParameters::LoadData(const G4String& name)
0103 {
0104   // Build the complete string identifying the file with the data set
0105
0106   // define active elements
0107
0108   const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
0109   if (materialTable == 0)
0110      G4Exception("G4RDBremsstrahlungParameters::LoadData()",
0111                  "DataNotFound", FatalException, "No MaterialTable found!");
0112
0113   G4int nMaterials = G4Material::GetNumberOfMaterials();
0114
0115   G4double x = 1.e-9;
0116   for (G4int mm=0; mm<100; mm++) {
0117     paramC.push_back(x);
0118   }
0119
0120   for (G4int m=0; m<nMaterials; m++) {
0121
0122     const G4Material* material= (*materialTable)[m];
0123     const G4ElementVector* elementVector = material->GetElementVector();
0124     const G4int nElements = material->GetNumberOfElements();
0125
0126     for (G4int iEl=0; iEl<nElements; iEl++) {
0127       G4Element* element = (*elementVector)[iEl];
0128       G4double Z = element->GetZ();
0129       G4int iz = (G4int)Z;
0130       if(iz < 100)
0131             paramC[iz] = 0.217635e-33*(material->GetTotNbOfElectPerVolume());
0132       if (!(activeZ.contains(Z))) {
0133      activeZ.push_back(Z);
0134       }
0135     }
0136   }
0137
0138   // Read parameters
0139
0140   const char* path = G4FindDataDir("G4LEDATA");
0141   if (path == 0)
0142     {
0143       G4String excep("G4LEDATA environment variable not set!");
0144       G4Exception("G4RDBremsstrahlungParameters::LoadData()",
0145                   "InvalidSetup", FatalException, excep);
0146     }
0147
0148   G4String pathString_a(path);
0149   G4String name_a = pathString_a + name;
0150   std::ifstream file_a(name_a);
0151   std::filebuf* lsdp_a = file_a.rdbuf();
0152
0153   if (! (lsdp_a->is_open()) )
0154     {
0155       G4String stringConversion2("Cannot open file ");
0156       G4String excep = stringConversion2 + name_a;
0157       G4Exception("G4RDBremsstrahlungParameters::LoadData()",
0158                   "FileNotFound", FatalException, excep);
0159   }
0160
0161   // The file is organized into two columns:
0162   // 1st column is the energy
0163   // 2nd column is the corresponding value
0164   // The file terminates with the pattern: -1   -1
0165   //                                       -2   -2
0166
0167   G4double ener = 0.0;
0168   G4double sum = 0.0;
0169   G4int z    = 0;
0170
0171   std::vector<G4DataVector*> a;
0172   for (size_t j=0; j<length; j++) {
0173     G4DataVector* aa = new G4DataVector();
0174     a.push_back(aa);
0175   }
0176   G4DataVector e;
0177   e.clear();
0178
0179   do {
0180     file_a >> ener >> sum;
0181
0182     // End of file
0183     if (ener == (G4double)(-2)) {
0184       break;
0185
0186       // End of next element
0187     } else if (ener == (G4double)(-1)) {
0188
0189       z++;
0190       G4double Z = (G4double)z;
0191
0192     // fill map if Z is used
0193       if (activeZ.contains(Z)) {
0194
0195     for (size_t k=0; k<length; k++) {
0196
0197         G4int id = z*length + k;
0198         G4RDVDataSetAlgorithm* inter  = new G4RDLogLogInterpolation();
0199         G4DataVector* eVector = new G4DataVector;
0200         size_t eSize = e.size();
0201         for (size_t s=0; s<eSize; s++) {
0202            eVector->push_back(e[s]);
0203         }
0204         G4RDVEMDataSet* set = new G4RDEMDataSet(id,eVector,a[k],inter,1.,1.);
0205             param[id] = set;
0206     }
0207         a.clear();
0208         for (size_t j=0; j<length; j++) {
0209             G4DataVector* aa = new G4DataVector();
0210             a.push_back(aa);
0211         }
0212       } else {
0213         for (size_t j=0; j<length; j++) {
0214           a[j]->clear();
0215         }
0216       }
0217       e.clear();
0218
0219     } else {
0220
0221       if(ener > 1000.) ener = 1000.;
0222       e.push_back(ener);
0223       a[length-1]->push_back(sum);
0224
0225       for (size_t j=0; j<length-1; j++) {
0226     G4double qRead;
0227     file_a >> qRead;
0228     a[j]->push_back(qRead);
0229       }
0230
0231     }
0232   } while (ener != (G4double)(-2));
0233
0234   file_a.close();
0235
0236 }
0237
0238
0239 G4double G4RDBremsstrahlungParameters::ParameterC(G4int id) const
0240 {
0241   G4int n = paramC.size();
0242   if (id < 0 || id >= n)
0243     {
0244       G4String stringConversion1("Wrong id = ");
0245       G4String stringConversion2(id);
0246       G4String ex = stringConversion1 + stringConversion2;
0247       G4Exception("G4RDBremsstrahlungParameters::ParameterC()",
0248                   "InvalidSetup", FatalException, ex);
0249     }
0250
0251   return paramC[id];
0252 }
0253
0254
0255 void G4RDBremsstrahlungParameters::PrintData() const
0256 {
0257
0258   G4cout << G4endl;
0259   G4cout << "===== G4RDBremsstrahlungParameters =====" << G4endl;
0260   G4cout << G4endl;
0261   G4cout << "===== Parameters =====" << G4endl;
0262   G4cout << G4endl;
0263
0264   size_t nZ = activeZ.size();
0265   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0266
0267   for (size_t j=0; j<nZ; j++) {
0268     G4int Z = (G4int)activeZ[j];
0269
0270     for (size_t i=0; i<length; i++) {
0271
0272       pos = param.find(Z*length + i);
0273       if (pos!= param.end()) {
0274
0275         G4cout << "===== Z= " << Z
0276                  << " parameter[" << i << "]  ====="
0277                  << G4endl;
0278         G4RDVEMDataSet* dataSet = (*pos).second;
0279         dataSet->PrintData();
0280       }
0281     }
0282   }
0283
0284   G4cout << "==========================================" << G4endl;
0285 }
0286