Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:08 

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
0033 //
0034 // -------------------------------------------------------------------
0035 
0036 #include "G4RDShellData.hh"
0037 #include "G4DataVector.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include <fstream>
0040 #include <sstream>
0041 #include <numeric>
0042 #include <algorithm>
0043 #include <valarray>
0044 #include <functional>
0045 #include "Randomize.hh"
0046 
0047 // The following deprecated header is included because <functional> seems not to be found on MGP's laptop
0048 //#include "function.h"
0049 
0050 // Constructor
0051 
0052 G4RDShellData::G4RDShellData(G4int minZ, G4int maxZ, G4bool isOccupancy)
0053   : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
0054 {  }
0055 
0056 // Destructor
0057 G4RDShellData::~G4RDShellData()
0058 {
0059   std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos;
0060   for (pos = idMap.begin(); pos != idMap.end(); ++pos)
0061     {
0062       std::vector<G4double>* dataSet = (*pos).second;
0063       delete dataSet;
0064     }
0065 
0066   std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos2;
0067   for (pos2 = bindingMap.begin(); pos2 != bindingMap.end(); ++pos2)
0068     {
0069       G4DataVector* dataSet = (*pos2).second;
0070       delete dataSet;
0071     }
0072 
0073   if (occupancyData)
0074     {
0075       std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos3;
0076       for (pos3 = occupancyPdfMap.begin(); pos3 != occupancyPdfMap.end(); ++pos3)
0077     {
0078       std::vector<G4double>* dataSet = (*pos3).second;
0079       delete dataSet;
0080     }
0081     }
0082 }
0083 
0084 
0085 size_t G4RDShellData::NumberOfShells(G4int Z) const
0086 {
0087   G4int z = Z - 1;
0088   G4int n = 0;
0089 
0090   if (Z>= zMin && Z <= zMax)
0091     {
0092       n = nShells[z];
0093     }
0094   return n;
0095 }
0096 
0097 
0098 const std::vector<G4double>& G4RDShellData::ShellIdVector(G4int Z) const
0099 {
0100   std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
0101   if (Z < zMin || Z > zMax)
0102     G4Exception("G4RDShellData::ShellIdVector()", "OutOfRange",
0103                 FatalException, "Z outside boundaries!");
0104   pos = idMap.find(Z);
0105   std::vector<G4double>* dataSet = (*pos).second;
0106   return *dataSet;
0107 }
0108 
0109 
0110 const std::vector<G4double>& G4RDShellData::ShellVector(G4int Z) const
0111 {
0112   std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
0113   if (Z < zMin || Z > zMax)
0114     G4Exception("G4RDShellData::ShellVector()", "OutOfRange",
0115                 FatalException, "Z outside boundaries!");
0116   pos = occupancyPdfMap.find(Z);
0117   std::vector<G4double>* dataSet = (*pos).second;
0118   return *dataSet;
0119 }
0120 
0121 
0122 G4int G4RDShellData::ShellId(G4int Z, G4int shellIndex) const
0123 {
0124   G4int n = -1;
0125 
0126   if (Z >= zMin && Z <= zMax)
0127     {
0128       std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
0129       pos = idMap.find(Z);
0130       if (pos!= idMap.end())
0131     {
0132       std::vector<G4double> dataSet = *((*pos).second);
0133       G4int nData = dataSet.size();
0134       if (shellIndex >= 0 && shellIndex < nData)
0135         {
0136           n = (G4int) dataSet[shellIndex];
0137         }
0138     }
0139     }
0140   return n;
0141 }
0142 
0143 
0144 G4double G4RDShellData::ShellOccupancyProbability(G4int Z, G4int shellIndex) const
0145 {
0146   G4double prob = -1.;
0147 
0148   if (Z >= zMin && Z <= zMax)
0149     {
0150       std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
0151       pos = idMap.find(Z);
0152       if (pos!= idMap.end())
0153     {
0154       std::vector<G4double> dataSet = *((*pos).second);
0155       G4int nData = dataSet.size();
0156       if (shellIndex >= 0 && shellIndex < nData)
0157         {
0158           prob = dataSet[shellIndex];
0159         }
0160     }
0161     }
0162   return prob;
0163 }
0164 
0165 
0166 
0167 G4double G4RDShellData::BindingEnergy(G4int Z, G4int shellIndex)  const
0168 {
0169   G4double value = 0.;
0170 
0171   if (Z >= zMin && Z <= zMax)
0172     {
0173       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
0174       pos = bindingMap.find(Z);
0175       if (pos!= bindingMap.end())
0176     {
0177       G4DataVector dataSet = *((*pos).second);
0178       G4int nData = dataSet.size();
0179       if (shellIndex >= 0 && shellIndex < nData)
0180         {
0181           value = dataSet[shellIndex];
0182         }
0183     }
0184     }
0185   return value;
0186 }
0187 
0188 void G4RDShellData::PrintData() const
0189 {
0190   for (G4int Z = zMin; Z <= zMax; Z++)
0191     {
0192       G4cout << "---- Shell data for Z = "
0193          << Z
0194          << " ---- "
0195          << G4endl;
0196       G4int nSh = nShells[Z-1];
0197       std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posId;
0198       posId = idMap.find(Z);
0199       std::vector<G4double>* ids = (*posId).second;
0200       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posE;
0201       posE = bindingMap.find(Z);
0202       G4DataVector* energies = (*posE).second;
0203       for (G4int i=0; i<nSh; i++)
0204     {
0205       G4int id = (G4int) (*ids)[i];
0206       G4double e = (*energies)[i] / keV;
0207       G4cout << i << ") ";
0208 
0209       if (occupancyData) 
0210         {
0211           G4cout << " Occupancy: ";
0212         }
0213       else 
0214         {
0215           G4cout << " Shell id: ";
0216         }
0217       G4cout << id << " - Binding energy = "
0218          << e << " keV ";
0219         if (occupancyData)
0220           {
0221         std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posOcc;
0222         posOcc = occupancyPdfMap.find(Z);
0223                 std::vector<G4double> probs = *((*posOcc).second);
0224                 G4double prob = probs[i];
0225         G4cout << "- Probability = " << prob;
0226           }
0227         G4cout << G4endl;
0228     }
0229       G4cout << "-------------------------------------------------" 
0230          << G4endl;
0231     }
0232 }
0233 
0234 
0235 void G4RDShellData::LoadData(const G4String& fileName)
0236 { 
0237   // Build the complete string identifying the file with the data set
0238   
0239   std::ostringstream ost;
0240   
0241   ost << fileName << ".dat";
0242   
0243   G4String name(ost.str());
0244   
0245   const char* path = G4FindDataDir("G4LEDATA");
0246   if (!path)
0247     { 
0248       G4String excep("G4LEDATA environment variable not set!");
0249       G4Exception("G4RDShellData::LoadData()", "InvalidSetup",
0250                   FatalException, excep);
0251     }
0252   
0253   G4String pathString(path);
0254   G4String dirFile = pathString + name;
0255   std::ifstream file(dirFile);
0256   std::filebuf* lsdp = file.rdbuf();
0257 
0258   if (! (lsdp->is_open()) )
0259     {
0260       G4String s1("Data file: ");
0261       G4String s2(" not found");
0262       G4String excep = s1 + dirFile + s2;
0263       G4Exception("G4RDShellData::LoadData()", "DataNotFound",
0264                   FatalException, excep);
0265     }
0266 
0267   G4double a = 0;
0268   G4int k = 1;
0269   G4int s = 0;
0270   
0271   G4int Z = 1;
0272   G4DataVector* energies = new G4DataVector;
0273   std::vector<G4double>* ids = new std::vector<G4double>;
0274 
0275   do {
0276     file >> a;
0277     G4int nColumns = 2;
0278     if (a == -1)
0279       {
0280     if (s == 0)
0281       {
0282         // End of a shell data set
0283         idMap[Z] = ids;
0284             bindingMap[Z] = energies;
0285             G4int n = ids->size();
0286         nShells.push_back(n);
0287         // Start of new shell data set
0288         ids = new std::vector<G4double>;
0289             energies = new G4DataVector;
0290             Z++;        
0291       }      
0292     s++;
0293     if (s == nColumns)
0294     {
0295       s = 0;
0296     }
0297       }
0298     else if (a == -2)
0299       {
0300     // End of file; delete the empty vectors created when encountering the last -1 -1 row
0301     delete energies;
0302     delete ids;
0303     //nComponents = components.size();
0304       }
0305     else
0306       {
0307     // 1st column is shell id
0308     if(k%nColumns != 0)
0309       {     
0310         ids->push_back(a);
0311         k++;
0312       }
0313     else if (k%nColumns == 0)
0314       {
0315         // 2nd column is binding energy
0316         G4double e = a * MeV;
0317         energies->push_back(e);
0318         k = 1;
0319       }
0320       }
0321   } while (a != -2); // end of file
0322   file.close();    
0323 
0324   // For Doppler broadening: the data set contains shell occupancy and binding energy for each shell
0325   // Build additional map with probability for each shell based on its occupancy
0326 
0327   if (occupancyData)
0328     {
0329       // Build cumulative from raw shell occupancy
0330 
0331       for (G4int Z=zMin; Z <= zMax; Z++)
0332     {
0333       std::vector<G4double> occupancy = ShellIdVector(Z);
0334 
0335       std::vector<G4double>* prob = new std::vector<G4double>;
0336       G4double scale = 1. / G4double(Z);
0337 
0338       prob->push_back(occupancy[0] * scale);
0339       for (size_t i=1; i<occupancy.size(); i++)
0340         {
0341           prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
0342         }
0343       occupancyPdfMap[Z] = prob;
0344 
0345       /*
0346         G4double scale = 1. / G4double(Z);
0347         //      transform((*prob).begin(),(*prob).end(),(*prob).begin(),bind2nd(multiplies<G4double>(),scale));
0348 
0349         for (size_t i=0; i<occupancy.size(); i++)
0350         {
0351         (*prob)[i] *= scale;
0352         }
0353       */
0354     }
0355     }
0356 }
0357 
0358 
0359 G4int G4RDShellData::SelectRandomShell(G4int Z) const
0360 {
0361   if (Z < zMin || Z > zMax)
0362     G4Exception("G4RDShellData::SelectRandomShell()", "OutOfRange",
0363                 FatalException, "Z outside boundaries!");
0364 
0365   G4int shellIndex = 0;    
0366   std::vector<G4double> prob = ShellVector(Z);
0367   G4double random = G4UniformRand();
0368 
0369   // std::vector<G4double>::const_iterator pos;
0370   // pos = lower_bound(prob.begin(),prob.end(),random);
0371 
0372   // Binary search the shell with probability less or equal random
0373 
0374   G4int nShells = NumberOfShells(Z);
0375   G4int upperBound = nShells;
0376 
0377   while (shellIndex <= upperBound) 
0378     {
0379       G4int midShell = (shellIndex + upperBound) / 2;
0380       if ( random < prob[midShell] ) 
0381     upperBound = midShell - 1;
0382       else 
0383     shellIndex = midShell + 1;
0384     }  
0385   if (shellIndex >= nShells) shellIndex = nShells - 1;
0386 
0387   return shellIndex;
0388 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!