physics/src/G4RDeIonisationParameters.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 //
0030 // History:
0031 // -----------
0032 // 31 Jul 2001   MGP        Created, with dummy implementation
0033 // 12.09.01 V.Ivanchenko    Add param and interpolation of parameters
0034 // 04.10.01 V.Ivanchenko    Add BindingEnergy method
0035 // 25.10.01 MGP             Many bug fixes, mostly related to the
0036 //                          management of pointers
0037 // 29.11.01 V.Ivanchenko    New parametrisation + Excitation
0038 // 30.05.02 V.Ivanchenko    Format and names of the data files were
0039 //                          chenged to "ion-..."
0040 // 17.02.04 V.Ivanchenko    Increase buffer size
0041 //
0042 // -------------------------------------------------------------------
0043
0044 #include <fstream>
0045 #include <sstream>
0046
0047 #include "G4RDeIonisationParameters.hh"
0048 #include "G4RDVEMDataSet.hh"
0049 #include "G4RDShellEMDataSet.hh"
0050 #include "G4RDEMDataSet.hh"
0051 #include "G4RDCompositeEMDataSet.hh"
0052 #include "G4RDLinInterpolation.hh"
0053 #include "G4RDLogLogInterpolation.hh"
0054 #include "G4RDLinLogLogInterpolation.hh"
0055 #include "G4RDSemiLogInterpolation.hh"
0056 #include "G4SystemOfUnits.hh"
0057 #include "G4Material.hh"
0058 #include "G4DataVector.hh"
0059
0060 G4RDeIonisationParameters:: G4RDeIonisationParameters(G4int minZ, G4int maxZ)
0061   : zMin(minZ), zMax(maxZ),
0062   length(24)
0063 {
0064   LoadData();
0065 }
0066
0067
0068 G4RDeIonisationParameters::~G4RDeIonisationParameters()
0069 {
0070   // Reset the map of data sets: remove the data sets from the map
0071   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
0072
0073   for (pos = param.begin(); pos != param.end(); ++pos)
0074     {
0075       G4RDVEMDataSet* dataSet = (*pos).second;
0076       delete dataSet;
0077     }
0078
0079   for (pos = excit.begin(); pos != excit.end(); ++pos)
0080     {
0081       G4RDVEMDataSet* dataSet = (*pos).second;
0082       delete dataSet;
0083     }
0084
0085   activeZ.clear();
0086 }
0087
0088
0089 G4double G4RDeIonisationParameters::Parameter(G4int Z, G4int shellIndex,
0090                                                      G4int parameterIndex,
0091                                                      G4double e) const
0092 {
0093   G4double value = 0.;
0094   G4int id = Z*100 + parameterIndex;
0095   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0096
0097   pos = param.find(id);
0098   if (pos!= param.end()) {
0099     G4RDVEMDataSet* dataSet = (*pos).second;
0100     G4int nShells = dataSet->NumberOfComponents();
0101
0102     if(shellIndex < nShells) {
0103       const G4RDVEMDataSet* component = dataSet->GetComponent(shellIndex);
0104       const G4DataVector ener = component->GetEnergies(0);
0105       G4double ee = std::max(ener.front(),std::min(ener.back(),e));
0106       value = component->FindValue(ee);
0107     } else {
0108       G4cout << "WARNING: G4IonisationParameters::FindParameter "
0109              << "has no parameters for shell= " << shellIndex
0110              << "; Z= " << Z
0111              << G4endl;
0112     }
0113   } else {
0114     G4cout << "WARNING: G4IonisationParameters::Parameter "
0115            << "did not find ID = "
0116            << shellIndex << G4endl;
0117   }
0118
0119   return value;
0120 }
0121
0122 G4double G4RDeIonisationParameters::Excitation(G4int Z, G4double e) const
0123 {
0124   G4double value = 0.;
0125   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0126
0127   pos = excit.find(Z);
0128   if (pos!= excit.end()) {
0129     G4RDVEMDataSet* dataSet = (*pos).second;
0130
0131     const G4DataVector ener = dataSet->GetEnergies(0);
0132     G4double ee = std::max(ener.front(),std::min(ener.back(),e));
0133     value = dataSet->FindValue(ee);
0134   } else {
0135     G4cout << "WARNING: G4IonisationParameters::Excitation "
0136            << "did not find ID = "
0137            << Z << G4endl;
0138   }
0139
0140   return value;
0141 }
0142
0143
0144 void G4RDeIonisationParameters::LoadData()
0145 {
0146   // ---------------------------------------
0147   // Please document what are the parameters
0148   // ---------------------------------------
0149
0150   // define active elements
0151
0152   const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
0153   if (materialTable == 0)
0154      G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup",
0155                  FatalException, "No MaterialTable found!");
0156
0157   G4int nMaterials = G4Material::GetNumberOfMaterials();
0158
0159   for (G4int m=0; m<nMaterials; m++) {
0160
0161     const G4Material* material= (*materialTable)[m];
0162     const G4ElementVector* elementVector = material->GetElementVector();
0163     const size_t nElements = material->GetNumberOfElements();
0164
0165     for (size_t iEl=0; iEl<nElements; iEl++) {
0166       G4Element* element = (*elementVector)[iEl];
0167       G4double Z = element->GetZ();
0168       if (!(activeZ.contains(Z))) {
0169     activeZ.push_back(Z);
0170       }
0171     }
0172   }
0173
0174   const char* path = G4FindDataDir("G4LEDATA");
0175   if (!path)
0176     {
0177       G4String excep("G4LEDATA environment variable not set!");
0178       G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup",
0179                   FatalException, excep);
0180     }
0181
0182   G4String pathString(path);
0183   G4String path2("/ioni/ion-sp-");
0184   pathString += path2;
0185
0186   G4double energy, sum;
0187
0188   size_t nZ = activeZ.size();
0189
0190   for (size_t i=0; i<nZ; i++) {
0191
0192     G4int Z = (G4int)activeZ[i];
0193     std::ostringstream ost;
0194     ost << pathString << Z << ".dat";
0195     G4String name(ost.str());
0196
0197     std::ifstream file(name);
0198     std::filebuf* lsdp = file.rdbuf();
0199
0200     if (! (lsdp->is_open()) ) {
0201       G4String excep = G4String("Data file: ")
0202       + name + G4String(" not found. The version 1.# of G4LEDATA should be used");
0203       G4Exception("G4RDeIonisationParameters::LoadData()", "DataNotFound",
0204                   FatalException, excep);
0205     }
0206
0207     // The file is organized into:
0208     // For each shell there are two lines:
0209     //    1st line:
0210     // 1st column is the energy of incident e-,
0211     // 2d  column is the parameter of screan term;
0212     //    2d line:
0213     // 3 energy (MeV) subdividing different approximation area of the spectrum
0214     // 20 point on the spectrum
0215     // The shell terminates with the pattern: -1   -1
0216     // The file terminates with the pattern: -2   -2
0217
0218     std::vector<G4RDVEMDataSet*> p;
0219     for (size_t k=0; k<length; k++)
0220       {
0221     G4RDVDataSetAlgorithm* inter = new G4RDLinLogLogInterpolation();
0222     G4RDVEMDataSet* composite = new G4RDCompositeEMDataSet(inter,1.,1.);
0223     p.push_back(composite);
0224       }
0225
0226     G4int shell = 0;
0227     std::vector<G4DataVector*> a;
0228     for (size_t j=0; j<length; j++)
0229       {
0230     G4DataVector* aa = new G4DataVector();
0231     a.push_back(aa);
0232       }
0233     G4DataVector e;
0234     e.clear();
0235     do {
0236       file >> energy >> sum;
0237       if (energy == -2) break;
0238
0239       if (energy >  -1) {
0240         e.push_back(energy);
0241         a[0]->push_back(sum);
0242         for (size_t j=0; j<length-1; j++) {
0243       G4double qRead;
0244       file >> qRead;
0245       a[j + 1]->push_back(qRead);
0246     }
0247
0248       } else {
0249
0250         // End of set for a shell, fill the map
0251     for (size_t k=0; k<length; k++) {
0252
0253           G4RDVDataSetAlgorithm* interp;
0254           if(0 == k) interp  = new G4RDLinLogLogInterpolation();
0255           else       interp  = new G4RDLogLogInterpolation();
0256
0257       G4DataVector* eVector = new G4DataVector;
0258       size_t eSize = e.size();
0259       for (size_t s=0; s<eSize; s++) {
0260            eVector->push_back(e[s]);
0261       }
0262       G4RDVEMDataSet* set = new G4RDEMDataSet(shell,eVector,a[k],interp,1.,1.);
0263
0264       p[k]->AddComponent(set);
0265     }
0266
0267     // clear vectors
0268         for (size_t j2=0; j2<length; j2++) {
0269       a[j2] = new G4DataVector();
0270     }
0271         shell++;
0272         e.clear();
0273       }
0274     } while (energy > -2);
0275
0276     file.close();
0277
0278     for (size_t kk=0; kk<length; kk++)
0279       {
0280     G4int id = Z*100 + kk;
0281     param[id] = p[kk];
0282       }
0283   }
0284
0285   G4String pathString_a(path);
0286   G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
0287   std::ifstream file_a(name_a);
0288   std::filebuf* lsdp_a = file_a.rdbuf();
0289   G4String pathString_b(path);
0290   G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
0291   std::ifstream file_b(name_b);
0292   std::filebuf* lsdp_b = file_b.rdbuf();
0293
0294   if (! (lsdp_a->is_open()) ) {
0295      G4String excep = G4String("Cannot open file ")
0296                     + name_a;
0297      G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile",
0298                  FatalException, excep);
0299   }
0300   if (! (lsdp_b->is_open()) ) {
0301      G4String excep = G4String("Cannot open file ")
0302                     + name_b;
0303      G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile",
0304                  FatalException, excep);
0305   }
0306
0307   // The file is organized into two columns:
0308   // 1st column is the energy
0309   // 2nd column is the corresponding value
0310   // The file terminates with the pattern: -1   -1
0311   //                                       -2   -2
0312
0313   G4double ener, ener1, sig, sig1;
0314   G4int z    = 0;
0315
0316   G4DataVector e;
0317   e.clear();
0318   G4DataVector d;
0319   d.clear();
0320
0321   do {
0322     file_a >> ener >> sig;
0323     file_b >> ener1 >> sig1;
0324
0325     if(ener != ener1) {
0326       G4cout << "G4RDeIonisationParameters: problem in excitation data "
0327              << "ener= " << ener
0328              << " ener1= " << ener1
0329              << G4endl;
0330     }
0331
0332     // End of file
0333     if (ener == -2) {
0334       break;
0335
0336       // End of next element
0337     } else if (ener == -1) {
0338
0339       z++;
0340       G4double Z = (G4double)z;
0341
0342     // fill map if Z is used
0343       if (activeZ.contains(Z)) {
0344
0345     G4RDVDataSetAlgorithm* inter  = new G4RDLinInterpolation();
0346     G4DataVector* eVector = new G4DataVector;
0347     G4DataVector* dVector = new G4DataVector;
0348     size_t eSize = e.size();
0349     for (size_t s=0; s<eSize; s++) {
0350            eVector->push_back(e[s]);
0351            dVector->push_back(d[s]);
0352     }
0353     G4RDVEMDataSet* set = new G4RDEMDataSet(z,eVector,dVector,inter,1.,1.);
0354         excit[z] = set;
0355       }
0356       e.clear();
0357       d.clear();
0358
0359     } else {
0360
0361       e.push_back(ener*MeV);
0362       d.push_back(sig1*sig*barn*MeV);
0363     }
0364   } while (ener != -2);
0365
0366   file_a.close();
0367
0368 }
0369
0370
0371 void G4RDeIonisationParameters::PrintData() const
0372 {
0373   G4cout << G4endl;
0374   G4cout << "===== G4RDeIonisationParameters =====" << G4endl;
0375   G4cout << G4endl;
0376
0377   size_t nZ = activeZ.size();
0378   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0379
0380   for (size_t i=0; i<nZ; i++) {
0381     G4int Z = (G4int)activeZ[i];
0382
0383     for (size_t j=0; j<length; j++) {
0384
0385       G4int index = Z*100 + j;
0386
0387       pos = param.find(index);
0388       if (pos!= param.end()) {
0389         G4RDVEMDataSet* dataSet = (*pos).second;
0390         size_t nShells = dataSet->NumberOfComponents();
0391
0392         for (size_t k=0; k<nShells; k++) {
0393
0394           G4cout << "===== Z= " << Z << " shell= " << k
0395                  << " parameter[" << j << "]  ====="
0396                  << G4endl;
0397           const G4RDVEMDataSet* comp = dataSet->GetComponent(k);
0398           comp->PrintData();
0399     }
0400       }
0401     }
0402   }
0403   G4cout << "====================================" << G4endl;
0404 }
0405
0406