ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPFFFissionFS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPFFFissionFS.cc
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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 //
33 #include "G4ParticleHPManager.hh"
34 #include "G4SystemOfUnits.hh"
35 
37 {
38  std::map<G4int,std::map<G4double,std::map<G4int,G4double >* >* >::iterator it = FissionProductYieldData.begin();
39  while ( it != FissionProductYieldData.end() ) { // Loop checking, 11.05.2015, T. Koi
40  std::map<G4double,std::map<G4int,G4double>* >* firstLevel = it->second;
41  if ( firstLevel ) {
42  std::map<G4double,std::map<G4int,G4double>*>::iterator it2 = firstLevel->begin();
43  while ( it2 != firstLevel->end() ) { // Loop checking, 11.05.2015, T. Koi
44  delete it2->second;
45  it2->second = 0;
46  firstLevel->erase(it2);
47  it2=firstLevel->begin();
48  }
49  }
50  delete firstLevel;
51  it->second = 0;
52  FissionProductYieldData.erase(it);
53  it = FissionProductYieldData.begin();
54  }
55 
56  std::map< G4int , std::map< G4double , G4int >* >::iterator ii = mMTInterpolation.begin();
57  while ( ii != mMTInterpolation.end() ) { // Loop checking, 11.05.2015, T. Koi
58  delete ii->second;
59  mMTInterpolation.erase(ii);
60  ii = mMTInterpolation.begin();
61  }
62 }
63 
65 {
66  //G4cout << "G4ParticleHPFFFissionFS::Init" << G4endl;
67  G4String aString = "FF";
68 
69  G4String tString = dirName;
70  G4bool dbool;
71  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aString , dbool);
72  G4String filename = aFile.GetName();
73  theBaseA = aFile.GetA();
74  theBaseZ = aFile.GetZ();
75 
76 //3456
77  if ( !dbool || ( Z < 2.5 && ( std::abs(theBaseZ-Z)>0.0001 || std::abs(theBaseA-A)>0.0001) ) )
78  {
79  hasAnyData = false;
80  hasFSData = false;
81  hasXsec = false;
82  return; // no data for exactly this isotope.
83  }
84  //std::ifstream theData(filename, std::ios::in);
85  std::istringstream theData(std::ios::in);
87  G4double dummy;
88  if ( !theData )
89  {
90  //theData.close();
91  hasFSData = false;
92  hasXsec = false;
93  hasAnyData = false;
94  return; // no data for this FS for this isotope
95  }
96 
97 
98  hasFSData = true;
99  // MT Energy FPS Yield
100  //std::map< int , std::map< double , std::map< int , double >* >* > FisionProductYieldData;
101  while ( theData.good() ) // Loop checking, 11.05.2015, T. Koi
102  {
103  G4int iMT, iMF;
104  G4int imax;
105  //Reading the data
106  // MT MF AWR
107  theData >> iMT >> iMF >> dummy;
108  // nBlock
109  theData >> imax;
110  //if ( !theData.good() ) continue;
111  // Ei FPS Yield
112  std::map< G4double , std::map< G4int , G4double >* >* mEnergyFSPData = new std::map< G4double , std::map< G4int , G4double >* >;
113 
114  std::map< G4double , G4int >* mInterporation = new std::map< G4double , G4int >;
115  for ( G4int i = 0 ; i <= imax ; i++ )
116  {
117 
118  G4double YY=0.0;
119  G4double Ei;
120  G4int jmax;
121  G4int ip;
122  // energy of incidence neutron
123  theData >> Ei;
124  // Number of data set followings
125  theData >> jmax;
126  // interpolation scheme
127  theData >> ip;
128  mInterporation->insert( std::pair<G4double,G4int>(Ei*eV,ip) );
129  // nNumber nIP
130  std::map<G4int,G4double>* mFSPYieldData = new std::map<G4int,G4double>;
131  for ( G4int j = 0 ; j < jmax ; j++ )
132  {
133  G4int FSP;
134  G4int mFSP;
135  G4double Y;
136  theData >> FSP >> mFSP >> Y;
137  G4int k = FSP*100+mFSP;
138  YY = YY + Y;
139  //if ( iMT == 454 )G4cout << iMT << " " << i << " " << j << " " << k << " " << Y << " " << YY << G4endl;
140  mFSPYieldData->insert( std::pair<G4int,G4double>( k , YY ) );
141  }
142  mEnergyFSPData->insert( std::pair<G4double,std::map<G4int,G4double>*>(Ei*eV,mFSPYieldData) );
143  }
144 
145  FissionProductYieldData.insert( std::pair< G4int , std::map< G4double , std::map< G4int , G4double >* >* > (iMT,mEnergyFSPData));
146  mMTInterpolation.insert( std::pair<G4int,std::map<G4double,G4int>*> (iMT,mInterporation) );
147  }
148  //theData.close();
149 }
150 
152 {
153  G4DynamicParticleVector * aResult;
154 // G4cout <<"G4ParticleHPFFFissionFS::ApplyYourself +"<<G4endl;
155  aResult = G4ParticleHPFissionBaseFS::ApplyYourself(nNeutrons);
156  return aResult;
157 }
158 
160 {
161  //G4cout << "G4ParticleHPFFFissionFS::GetAFissionFragment " << G4endl;
162 
163  G4double rand =G4UniformRand();
164  //G4cout << rand << G4endl;
165 
166  std::map< G4double , std::map< G4int , G4double >* >* mEnergyFSPData = FissionProductYieldData.find( 454 )->second;
167 
168  //It is not clear that the treatment of the scheme 2 on two-dimensional interpolation.
169  //So, here just use the closest energy point array of yield data.
170  //TK120531
171  G4double key_energy = DBL_MAX;
172  if ( mEnergyFSPData->size() == 1 )
173  {
174  key_energy = mEnergyFSPData->begin()->first;
175  }
176  else
177  {
178  //Find closest energy point
179  G4double Dmin=DBL_MAX;
180  G4int i = 0;
181  for ( std::map< G4double , std::map< G4int , G4double >* >::iterator it = mEnergyFSPData->begin() ;
182  it != mEnergyFSPData->end() ; it++ )
183  {
184  G4double e = (it->first);
185  G4double d = std::fabs ( energy - e );
186  if ( d < Dmin )
187  {
188  Dmin = d;
189  key_energy = e;
190  }
191  i++;
192  }
193  }
194 
195  std::map<G4int,G4double>* mFSPYieldData = (*mEnergyFSPData)[key_energy];
196 
197  G4int ifrag=0;
198  G4double ceilling = mFSPYieldData->rbegin()->second; // Because of numerical accuracy, this is not always 2
199  for ( std::map<G4int,G4double>::iterator it = mFSPYieldData->begin() ; it != mFSPYieldData->end() ; it++ )
200  {
201  //if ( ( rand - it->second/ceilling ) < 1.0e-6 ) std::cout << rand - it->second/ceilling << std::endl;
202  if ( rand <= it->second/ceilling )
203  {
204  //G4cout << it->first << " " << it->second/ceilling << G4endl;
205  ifrag = it->first;
206  break;
207  }
208  }
209 
210  fragZ = ifrag/100000;
211  fragA = (ifrag%100000)/100;
212  fragM = (ifrag%100);
213 
214  //G4cout << fragZ << " " << fragA << " " << fragM << G4endl;
215 }