physics/src/G4RDAugerData.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 // Author: Alfonso Mmantero (Alfonso.Mantero@ge.infn.it)
0030 //
0031 // History:
0032 // -----------
0033 // Based on G4RDFluoData by Elena Guardincerri
0034 //
0035 // Modified: 30.07.02 VI Add select active Z + clean up against pedantic compiler
0036 //
0037 // -------------------------------------------------------------------
0038
0039 #include <fstream>
0040 #include <sstream>
0041
0042 #include "G4RDAugerData.hh"
0043 #include "G4PhysicalConstants.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4DataVector.hh"
0046 #include "G4Material.hh"
0047 #include "G4Element.hh"
0048 #include "G4ElementVector.hh"
0049
0050 G4RDAugerData::G4RDAugerData()
0051 {
0052
0053   G4int n = 0;
0054   G4int pos = 0;
0055
0056     for (pos = 0 ; pos < 100; pos++)
0057       {
0058     numberOfVacancies.push_back(n);
0059       }
0060
0061   BuildAugerTransitionTable();
0062
0063
0064 }
0065
0066 G4RDAugerData::~G4RDAugerData()
0067 {
0068   /*
0069   std::map<G4int,std::vector<G4RDAugerTransition>,std::less<G4int> >::iterator pos;
0070
0071   for (pos = augerTransitionTable.begin(); pos != augerTransitionTable.end(); pos++)
0072     {
0073       std::vector<G4RDAugerTransition> dataSet = (*pos).second;
0074       delete dataSet;
0075     }
0076   for (pos = energyMap.begin(); pos != energyMap.end(); pos++)
0077     {
0078       std::map<G4Int,G4DataVector*,std::less<G4int>>* dataMap = (*pos).second;
0079       for (pos2 = newIdProbabilityMap.begin(); pos2 != idMap.end(); pos2++)
0080     {
0081       G4DataVector* dataSet = (*pos2).second;
0082       delete dataSet;
0083     }
0084     }
0085   for (pos = probabilityMap.begin(); pos != probabilityMap.end(); pos++)
0086     {
0087       std::map<G4Int,G4DataVector*,std::less<G4int>>* dataMap = (*pos).second;
0088       for (pos2 = newIdProbabilityMap.begin(); pos2 != idMap.end(); pos2++)
0089     {
0090       G4DataVector* dataSet = (*pos2).second;
0091       delete dataSet;
0092     }
0093     }
0094   for (pos2 = newIdMap.begin(); pos2 != idMap.end(); pos2++)
0095     {
0096       G4DataVector* dataSet = (*pos2).second;
0097       delete dataSet;
0098     }
0099   for (pos2 = newIdEnergyMap.begin(); pos2 != idMap.end(); pos2++)
0100     {
0101       G4DataVector* dataSet = (*pos2).second;
0102       delete dataSet;
0103     }
0104   for (pos2 = newIdProbabilityMap.begin(); pos2 != idMap.end(); pos2++)
0105     {
0106       G4DataVector* dataSet = (*pos2).second;
0107       delete dataSet;
0108     }
0109   */
0110
0111 }
0112
0113 size_t G4RDAugerData::NumberOfVacancies(G4int Z) const
0114 {
0115   return numberOfVacancies[Z];
0116 }
0117
0118 G4int G4RDAugerData::VacancyId(G4int Z, G4int vacancyIndex) const
0119 {
0120
0121   G4int n = 0;
0122   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0123     {G4Exception("G4RDAugerData::VacancyId()", "OutOfRange",
0124                  FatalException, "VacancyIndex outside boundaries!");}
0125   else {
0126     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0127     if (element == augerTransitionTable.end())
0128       {G4Exception("G4RDAugerData::VacancyId()", "NoDataFound",
0129                    FatalException, "Data not loaded!");}
0130     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0131     n = (G4int) dataSet[vacancyIndex].FinalShellId();
0132   }
0133
0134   return n;
0135 }
0136
0137
0138
0139 // Attenzione: questa funzione vuole l'indice della vacanza, non l'Id
0140
0141
0142 size_t G4RDAugerData::NumberOfTransitions(G4int Z, G4int vacancyIndex) const
0143 {
0144   G4int n = 0;
0145   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0146     {G4Exception("G4RDAugerData::NumberOfTransitions()", "OutOfRange",
0147                  FatalException, "VacancyIndex outside boundaries!");}
0148   else {
0149     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0150     if (element == augerTransitionTable.end())
0151       {G4Exception("G4RDAugerData::NumberOfTransitions()", "NoDataFound",
0152                    FatalException, "Data not loaded!");}
0153     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0154     n = (G4int)dataSet[vacancyIndex].TransitionOriginatingShellIds()->size();
0155   }
0156  return  n;
0157 }
0158
0159
0160
0161 size_t G4RDAugerData::NumberOfAuger(G4int Z, G4int initIndex, G4int vacancyId) const
0162 {
0163   size_t n = 0;
0164   if (initIndex<0 || initIndex>=numberOfVacancies[Z])
0165     {G4Exception("G4RDAugerData::NumberOfAuger()", "OutOfRange",
0166                  FatalException, "VacancyIndex outside boundaries!");}
0167   else {
0168     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0169     if (element == augerTransitionTable.end())
0170       {G4Exception("G4RDAugerData::NumberOfAuger()", "NoDataFound",
0171                    FatalException, "Data not loaded!");}
0172     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0173     const std::vector<G4int>* temp =  dataSet[initIndex].AugerOriginatingShellIds(vacancyId);
0174     n = temp->size();
0175   }
0176   return n;
0177 }
0178
0179 size_t G4RDAugerData::AugerShellId(G4int Z, G4int vacancyIndex, G4int transId, G4int augerIndex) const
0180 {
0181   size_t n = 0;
0182   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0183     {G4Exception("G4RDAugerData::AugerShellId()", "OutOfRange",
0184                  FatalException, "VacancyIndex outside boundaries!");}
0185   else {
0186     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0187     if (element == augerTransitionTable.end())
0188       {G4Exception("G4RDAugerData::AugerShellId()", "NoDataFound",
0189                    FatalException, "Data not loaded!");}
0190     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0191     n = dataSet[vacancyIndex].AugerOriginatingShellId(augerIndex,transId);
0192   }
0193   return n;
0194 }
0195
0196 G4int G4RDAugerData::StartShellId(G4int Z, G4int vacancyIndex, G4int transitionShellIndex) const
0197 {
0198   G4int n = 0;
0199
0200   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0201     {G4Exception("G4RDAugerData::StartShellId()", "OutOfRange",
0202                  FatalException, "VacancyIndex outside boundaries!");}
0203   else {
0204     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0205     if (element == augerTransitionTable.end())
0206       {G4Exception("G4RDAugerData::StartShellId()", "NoDataFound",
0207                    FatalException, "Data not loaded!");}
0208     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0209      n = dataSet[vacancyIndex].TransitionOriginatingShellId(transitionShellIndex);
0210   }
0211
0212
0213  return n;
0214 }
0215
0216 G4double G4RDAugerData::StartShellEnergy(G4int Z, G4int vacancyIndex, G4int transitionId, G4int augerIndex) const
0217 {
0218   G4double energy = 0;
0219
0220   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0221     {G4Exception("G4RDAugerData::StartShellEnergy()", "OutOfRange",
0222                  FatalException, "VacancyIndex outside boundaries!");}
0223   else {
0224     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0225     if (element == augerTransitionTable.end())
0226       {G4Exception("G4RDAugerData::StartShellEnergy()", "NoDataFound",
0227                    FatalException, "Data not loaded!");}
0228     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0229     energy = dataSet[vacancyIndex].AugerTransitionEnergy(augerIndex,transitionId);
0230
0231   }
0232   return energy;
0233 }
0234
0235
0236 G4double G4RDAugerData::StartShellProb(G4int Z, G4int vacancyIndex,G4int transitionId,G4int augerIndex) const
0237 {
0238   G4double prob = 0;
0239
0240   if (vacancyIndex<0 || vacancyIndex>=numberOfVacancies[Z])
0241     {G4Exception("G4RDAugerData::StartShellProb()", "OutOfRange",
0242                  FatalException, "VacancyIndex outside boundaries!");}
0243   else {
0244     trans_Table::const_iterator element = augerTransitionTable.find(Z);
0245     if (element == augerTransitionTable.end())
0246       {G4Exception("G4RDAugerData::StartShellProb()", "NoDataFound",
0247                    FatalException, "Data not loaded!");}
0248     std::vector<G4RDAugerTransition> dataSet = (*element).second;
0249     prob = dataSet[vacancyIndex].AugerTransitionProbability(augerIndex, transitionId);
0250
0251
0252
0253   }
0254      return prob;
0255 }
0256
0257 std::vector<G4RDAugerTransition> G4RDAugerData::LoadData(G4int Z)
0258 {
0259   // Build the complete string identifying the file with the data set
0260
0261     std::ostringstream ost;
0262     if(Z != 0){
0263       ost << "au-tr-pr-"<< Z << ".dat";
0264     }
0265     else{
0266       ost << "au-tr-pr-"<<".dat";
0267     }
0268     G4String name(ost.str());
0269
0270     const char* path = G4FindDataDir("G4LEDATA");
0271     if (!path)
0272       {
0273     G4String excep = "G4LEDATA environment variable not set";
0274     G4Exception("G4RDAugerData::LoadData()", "InvalidSetup",
0275                     FatalException, excep);
0276       }
0277
0278     G4String pathString(path);
0279     G4String dirFile = pathString + "/auger/" + name;
0280     std::ifstream file(dirFile);
0281     std::filebuf* lsdp = file.rdbuf();
0282
0283     if (! (lsdp->is_open()) )
0284       {
0285     G4String excep = "Data file: " + dirFile + " not found!";
0286     G4Exception("G4RDAugerData::LoadData()", "DataNotFound",
0287                     FatalException, excep);
0288       }
0289
0290
0291     G4double a = 0;
0292     G4int k = 1;
0293     G4int s = 0;
0294
0295     G4int vacId = 0;
0296     std::vector<G4int>* initIds = new std::vector<G4int>;
0297     std::vector<G4int>* newIds = new std::vector<G4int>;
0298     G4DataVector* transEnergies = new G4DataVector;
0299     G4DataVector* transProbabilities = new G4DataVector;
0300     std::vector<G4RDAugerTransition> augerTransitionVector;
0301     std::map<G4int,std::vector<G4int>,std::less<G4int> >* newIdMap =
0302       new std::map<G4int,std::vector<G4int>,std::less<G4int> >;
0303     std::map<G4int,G4DataVector,std::less<G4int> >* newEnergyMap =
0304       new std::map<G4int,G4DataVector,std::less<G4int> >;
0305     std::map<G4int,G4DataVector,std::less<G4int> >* newProbabilityMap =
0306       new std::map<G4int,G4DataVector,std::less<G4int> >;
0307
0308
0309     do {
0310       file >> a;
0311
0312
0313       G4int nColumns = 4;
0314
0315       if (a == -1)
0316     {
0317
0318
0319
0320       if (s == 0)
0321         {
0322           // End of a shell data set
0323
0324
0325
0326           std::vector<G4int>::iterator vectorIndex = initIds->begin();
0327
0328           vacId = *vectorIndex;
0329
0330           //initIds->erase(vectorIndex);
0331
0332
0333
0334           std::vector<G4int> identifiers;
0335           for (vectorIndex = initIds->begin()+1 ; vectorIndex != initIds->end(); ++vectorIndex){
0336
0337         identifiers.push_back(*vectorIndex);
0338           }
0339
0340           vectorIndex = (initIds->end())-1;
0341
0342           G4int augerShellId = *(vectorIndex);
0343
0344
0345           (*newIdMap)[augerShellId] = *newIds;
0346           (*newEnergyMap)[augerShellId] = *transEnergies;
0347           (*newProbabilityMap)[augerShellId] = *transProbabilities;
0348
0349           augerTransitionVector.push_back(G4RDAugerTransition(vacId, identifiers, newIdMap, newEnergyMap, newProbabilityMap));
0350
0351           // Now deleting all the variables I used, and creating new ones for the next shell
0352
0353           delete newIdMap;
0354           delete newEnergyMap;
0355           delete newProbabilityMap;
0356
0357           G4int n = initIds->size();
0358           nInitShells.push_back(n);
0359           numberOfVacancies[Z]++;
0360           delete initIds;
0361           delete newIds;
0362           delete transEnergies;
0363           delete transProbabilities;
0364           initIds = new std::vector<G4int>;
0365           newIds = new std::vector<G4int>;
0366           transEnergies = new G4DataVector;
0367           transProbabilities = new G4DataVector;
0368           newIdMap = new std::map<G4int,std::vector<G4int>,std::less<G4int> >;
0369           newEnergyMap = new std::map<G4int,G4DataVector,std::less<G4int> >;
0370           newProbabilityMap = new std::map<G4int,G4DataVector,std::less<G4int> >;
0371
0372
0373
0374         }
0375       s++;
0376       if (s == nColumns)
0377         {
0378           s = 0;
0379         }
0380     }
0381       else if (a == -2)
0382     {
0383       // End of file; delete the empty vectors created
0384       //when encountering the last -1 -1 row
0385       delete initIds;
0386       delete newIds;
0387       delete transEnergies;
0388       delete transProbabilities;
0389       delete newIdMap ;
0390       delete newEnergyMap;
0391       delete newProbabilityMap;
0392     }
0393       else
0394     {
0395
0396       if (k%nColumns == 3){
0397         // 3rd column is the transition probabilities
0398         transProbabilities->push_back(a);
0399
0400         k++;}
0401       else if(k%nColumns == 2){
0402         // 2nd column is new auger vacancy
0403
0404         // 2nd column is new auger vacancy
0405
0406         G4int l = (G4int)a;
0407         newIds->push_back(l);
0408
0409
0410         k++;
0411       }
0412       else if (k%nColumns == 1)
0413         {
0414           // 1st column is shell id
0415
0416           if(initIds->size() == 0) {
0417
0418         // if this is the first data of the shell, all the colums are equal
0419         // to the shell Id; so we skip the next colums ang go to the next row
0420
0421         initIds->push_back((G4int)a);
0422         // first line of initIds is the original shell of the vacancy
0423         file >> a;
0424         file >> a;
0425         file >> a;
0426         k = k+3;
0427           }
0428           else {
0429
0430         //        std::vector<G4int>::iterator vectorIndex = (initIds->end())-1;
0431         if((G4int)a != initIds->back()){
0432
0433
0434           if((initIds->size()) == 1) {
0435             initIds->push_back((G4int)a);
0436           }
0437           else {
0438
0439
0440             G4int augerShellId = 0;
0441             augerShellId = initIds->back();
0442
0443             (*newIdMap)[augerShellId] = *newIds;
0444             (*newEnergyMap)[augerShellId] = *transEnergies;
0445             (*newProbabilityMap)[augerShellId] = *transProbabilities;
0446             delete newIds;
0447             delete transEnergies;
0448             delete transProbabilities;
0449             newIds = new std::vector<G4int>;
0450             transEnergies = new G4DataVector;
0451             transProbabilities = new G4DataVector;
0452             initIds->push_back((G4int)a);
0453           }
0454         }
0455           }
0456
0457           k++;
0458
0459         }
0460       else if (k%nColumns == 0)
0461
0462         {//fourth column is transition energies
0463           G4double e = a * MeV;
0464
0465           transEnergies->push_back(e);
0466           k=1;
0467
0468         }
0469     }
0470     }
0471
0472
0473     while (a != -2); // end of file
0474     file.close();
0475     return augerTransitionVector;
0476
0477 }
0478
0479 void G4RDAugerData::BuildAugerTransitionTable()
0480 {
0481
0482   //  trans_Table::iterator pos = augerTransitionTable.begin();
0483
0484   const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
0485
0486   G4int nMaterials = G4Material::GetNumberOfMaterials();
0487
0488   G4DataVector activeZ;
0489   activeZ.clear();
0490
0491   for (G4int m=0; m<nMaterials; m++) {
0492
0493     const G4Material* material= (*materialTable)[m];
0494     const G4ElementVector* elementVector = material->GetElementVector();
0495     const size_t nElements = material->GetNumberOfElements();
0496
0497     for (size_t iEl=0; iEl<nElements; iEl++) {
0498       G4Element* element = (*elementVector)[iEl];
0499       G4double Z = element->GetZ();
0500       if (!(activeZ.contains(Z))) {
0501     activeZ.push_back(Z);
0502       }
0503     }
0504   }
0505
0506
0507   for (G4int element = 6; element < 101; element++)
0508     {
0509       if(nMaterials == 0 || activeZ.contains(element)) {
0510         augerTransitionTable.insert(trans_Table::value_type(element,LoadData(element)));
0511
0512         G4cout << "G4RDAugerData for Element no. " << element << " are loaded" << G4endl;
0513       //      PrintData(element);
0514       }
0515     }
0516
0517   G4cout << "AugerTransitionTable complete"<< G4endl;
0518 }
0519
0520 void G4RDAugerData::PrintData(G4int Z)
0521 {
0522
0523   for (G4int i = 0; i < numberOfVacancies[Z]; i++)
0524     {
0525       G4cout << "---- TransitionData for the vacancy nb "
0526          <<i
0527          <<" of the atomic number elemnt "
0528          << Z
0529          <<"----- "
0530          <<G4endl;
0531
0532       for (size_t k = 0; k<=NumberOfTransitions(Z,i); k++)
0533     {
0534       G4int id = StartShellId(Z,i,k);
0535
0536       for (size_t a = 0; a <= NumberOfAuger(Z,i,id); a++) {
0537
0538         G4double e = StartShellEnergy(Z,i,id,a) /MeV;
0539         G4double p = StartShellProb(Z,i,id,a);
0540         G4int augerId = AugerShellId(Z, i, id, a);
0541         G4cout << k <<") Shell id: " << id <<G4endl;
0542         G4cout << "    Auger Originatig Shell Id :"<< augerId <<G4endl;
0543         G4cout << " - Transition energy = " << e << " MeV "<<G4endl;
0544         G4cout   << " - Transition probability = " << p <<G4endl;
0545       }
0546     }
0547       G4cout << "-------------------------------------------------"
0548          << G4endl;
0549     }
0550 }
0551 G4RDAugerTransition* G4RDAugerData::GetAugerTransition(G4int Z,G4int vacancyShellIndex)
0552     {
0553       std::vector<G4RDAugerTransition>* dataSet = &augerTransitionTable[Z];
0554       std::vector<G4RDAugerTransition>::iterator vectorIndex = dataSet->begin() + vacancyShellIndex;
0555
0556       G4RDAugerTransition* augerTransition = &(*vectorIndex);
0557       return augerTransition;
0558     }
0559
0560 std::vector<G4RDAugerTransition>* G4RDAugerData::GetAugerTransitions(G4int Z)
0561   {
0562     std::vector<G4RDAugerTransition>* dataSet = &augerTransitionTable[Z];
0563     return dataSet;
0564   }
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585