Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/abla/src/G4AblaDataFile.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // ABLAXX statistical de-excitation model
 27 // Jose Luis Rodriguez, UDC (translation from ABLA07 and contact person)
 28 // Pekka Kaitaniemi, HIP (initial translation of ablav3p)
 29 // Aleksandra Kelic, GSI (ABLA07 code)
 30 // Davide Mancusi, CEA (contact person INCL)
 31 // Aatos Heikkinen, HIP (project coordination)
 32 //
 33 
 34 #include "G4AblaDataFile.hh"
 35 #include "G4AblaDataDefs.hh"
 36 #include "globals.hh"
 37 
 38 #include <cmath>
 39 #include <cstdlib>
 40 #include <fstream>
 41 #include <iostream>
 42 
 43 G4AblaDataFile::G4AblaDataFile() { verboseLevel = 0; }
 44 
 45 /**
 46  * Read all data from files.
 47  */
 48 G4bool G4AblaDataFile::readData()
 49 {
 50     if (!G4FindDataDir("G4ABLADATA"))
 51     {
 52         G4ExceptionDescription ed;
 53         ed << " Data missing: set environment variable G4ABLADATA\n"
 54            << " to point to the directory containing data files needed\n"
 55            << " by the ABLA model" << G4endl;
 56         G4Exception("G4AblaDataFile::readData()", "ABLA_001", FatalException, ed);
 57     }
 58     G4String dataPath(G4FindDataDir("G4ABLADATA"));
 59 
 60     G4String flAlphaFile(dataPath + "/flalpha.dat");
 61     G4String frldmFile(dataPath + "/frldm.dat");
 62     G4String vgsldFile(dataPath + "/vgsld.dat");
 63     G4String rmsFile(dataPath + "/rms.dat");
 64     G4String defoFile(dataPath + "/defo.dat");
 65     G4String massFile(dataPath + "/mass2020.dat");
 66 
 67     if (verboseLevel > 1)
 68     {
 69         // G4cout <<"Data path   = " << dataPath    << G4endl;
 70         // G4cout <<"FlAlphaFile = " << flAlphaFile << G4endl;
 71         // G4cout <<"FrldmFile   = " << frldmFile   << G4endl;
 72         // G4cout <<"VgsldFile   = " << vgsldFile   << G4endl;
 73     }
 74 
 75     std::ifstream flalphain(flAlphaFile.c_str());
 76     std::ifstream frldmin(frldmFile.c_str());
 77     std::ifstream vgsldin(vgsldFile.c_str());
 78     std::ifstream rmsin(rmsFile.c_str());
 79     std::ifstream defoin(defoFile.c_str());
 80     std::ifstream massin(massFile.c_str());
 81 
 82     if (!massin.is_open())
 83     {
 84         massFile = dataPath + "/mass2016.dat";
 85         massin.close();
 86         massin.open(massFile.c_str());
 87         std::cout << "Mass evaluation file mass2020.dat not found, current file: " << massFile.c_str() << std::endl;
 88 
 89         if (!massin.is_open())
 90         {
 91             massFile = dataPath + "/mass2003.dat";
 92             massin.close();
 93             massin.open(massFile.c_str());
 94             std::cout << "Mass evaluation file mass2016.dat not found, current file: " << massFile.c_str() << std::endl;
 95         }
 96     }
 97 
 98     std::filebuf* buf1 = flalphain.rdbuf();
 99     std::filebuf* buf2 = frldmin.rdbuf();
100     std::filebuf* buf3 = vgsldin.rdbuf();
101     std::filebuf* buf4 = rmsin.rdbuf();
102     std::filebuf* buf5 = defoin.rdbuf();
103     std::filebuf* buf6 = massin.rdbuf();
104     if (!((buf1->is_open()) && (buf2->is_open()) && (buf3->is_open()) && (buf4->is_open()) && (buf5->is_open()) &&
105           (buf6->is_open())))
106     {
107         G4ExceptionDescription ed;
108         ed << "Data missing: could not find ABLA data file in " << dataPath
109            << "defined by environment variable G4ABLADATA" << G4endl;
110         G4Exception("G4AblaDataFile::readData()", "ABLA", FatalException, ed);
111     }
112 
113     G4double fflalpha, ffrldm, fvgsld, frms;
114     G4int fj = 0, fk = 0, a2, a3, a4;
115     G4double fbeta2, fbeta4;
116     G4double a7;
117     const G4int rows = 99;
118     const G4int cols = 154;
119     const G4int rowsbeta = 137;
120     const G4int colsbeta = 251;
121 
122     for (G4int i = 0; i < zcols; i++)
123     {
124         for (G4int j = 0; j < nrows; j++)
125         {
126             setAlpha(j, i, 0.0);
127             setEcnz(j, i, 0.0);
128             setVgsld(j, i, 0.0);
129             setRms(j, i, 0.0);
130         }
131     }
132 
133     for (G4int i = 0; i < rows; i++)
134     {
135         for (G4int j = 0; j < cols; j++)
136         {
137             flalphain >> fflalpha;
138             frldmin >> ffrldm;
139             vgsldin >> fvgsld;
140             rmsin >> frms;
141             setAlpha(j, i, fflalpha);
142             setEcnz(j, i, ffrldm);
143             setVgsld(j, i, fvgsld);
144             setRms(j, i, frms);
145         }
146     }
147 
148     for (G4int i = 0; i < rowsbeta; i++)
149     {
150         for (G4int j = 0; j < colsbeta; j++)
151         {
152             setBeta2(j, i, 0.0);
153             setBeta4(j, i, 0.0);
154         }
155     }
156 
157     defoin >> fj >> fk >> fbeta2 >> fbeta4;
158     while (!defoin.eof())
159     {
160         setBeta2(fk, fj, fbeta2);
161         setBeta4(fk, fj, fbeta4);
162         defoin >> fj >> fk >> fbeta2 >> fbeta4;
163     }
164 
165     for (G4int i = 0; i < zcols; i++)
166     {
167         for (G4int j = 0; j < nrows; j++)
168         {
169             setMexp(j, i, 0.0);
170             setMexpID(j, i, 0);
171         }
172     }
173     massin >> a2 >> a3 >> a4 >> a7;
174     while (!massin.eof())
175     {
176         //
177         if (a3 < lpcols)
178         {
179             setMexpID(a2, a3, 1);
180             setMexp(a2, a3, 938.7829835 * a3 + 939.5653301 * a2 - 1. * a4 * a7 / 1000.);
181         }
182         massin >> a2 >> a3 >> a4 >> a7;
183     }
184 
185     flalphain.close();
186     frldmin.close();
187     vgsldin.close();
188     rmsin.close();
189     defoin.close();
190     massin.close();
191 
192     return true;
193 }
194