ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPElementData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPElementData.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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 02-08-06 Modified Harmonise to reslove cross section trouble at high-end. T. KOI
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 #include "G4ParticleHPElementData.hh"
35 
36  G4ParticleHPElementData::G4ParticleHPElementData()
37  {
38  precision = 0.02;
39  theFissionData = new G4ParticleHPVector;
40  theCaptureData = new G4ParticleHPVector;
41  theElasticData = new G4ParticleHPVector;
42  theInelasticData = new G4ParticleHPVector;
43  theIsotopeWiseData = 0;
44  theBuffer = NULL;
45  }
46 
47  G4ParticleHPElementData::~G4ParticleHPElementData()
48  {
49  delete theFissionData;
50  delete theCaptureData;
51  delete theElasticData;
52  delete theInelasticData;
53  delete [] theIsotopeWiseData;
54  }
55 
56  void G4ParticleHPElementData::Init(G4Element * theElement, G4ParticleDefinition* projectile, const char* dataDirVariable )
57  {
58  G4int count = theElement->GetNumberOfIsotopes();
59  if(count == 0) count +=
60  theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()));
61  theIsotopeWiseData = new G4ParticleHPIsoData[count];
62  // filename = ein data-set je isotope.
63  count = 0;
64  G4int nIso = theElement->GetNumberOfIsotopes();
65  G4int Z = static_cast<G4int> (theElement->GetZ());
66  // G4int i1;
67  if(nIso!=0)
68  {
69  for (G4int i1=0; i1<nIso; i1++)
70  {
71 // G4cout <<" Init: normal case"<<G4endl;
72  G4int A = theElement->GetIsotope(i1)->GetN();
73  G4int M = theElement->GetIsotope(i1)->Getm();
74  G4double frac = theElement->GetRelativeAbundanceVector()[i1]/CLHEP::perCent;
75  //UpdateData(A, Z, count++, frac);
76  UpdateData(A, Z, M, count++, frac, projectile, dataDirVariable);
77  }
78  }else{
79 // G4cout <<" Init: theStableOnes case: Z="<<Z<<G4endl;
80  G4int first = theStableOnes.GetFirstIsotope(Z);
81 // G4cout <<"first="<<first<<" "<<theStableOnes.GetNumberOfIsotopes(theElement->GetZ())<<G4endl;
82  for(G4int i1=0;
83  i1<theStableOnes.GetNumberOfIsotopes(static_cast<G4int>(theElement->GetZ()) );
84  i1++)
85  {
86 // G4cout <<" Init: theStableOnes in the loop"<<G4endl;
87  G4int A = theStableOnes.GetIsotopeNucleonCount(first+i1);
88  G4double frac = theStableOnes.GetAbundance(first+i1);
89 // G4cout <<" Init: theStableOnes in the loop: "<<A<<G4endl;
90  UpdateData(A, Z, count++, frac, projectile, dataDirVariable);
91  }
92  }
93  theElasticData->ThinOut(precision);
94  if( projectile == G4Neutron::Neutron() ) theInelasticData->ThinOut(precision);
95 
96  theCaptureData->ThinOut(precision);
97  theFissionData->ThinOut(precision);
98  }
99 
100  //void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
101  void G4ParticleHPElementData::UpdateData(G4int A, G4int Z, G4int M, G4int index, G4double abundance, G4ParticleDefinition* projectile, const char* dataDirVariable )
102  {
103  //Reads in the Data, using G4ParticleHPIsoData[], and its Init
104 // G4cout << "entered: ElementWiseData::UpdateData"<<G4endl;
105  //theIsotopeWiseData[index].Init(A, Z, abundance);
106  theIsotopeWiseData[index].Init(A, Z, M, abundance,projectile, dataDirVariable);
107 // G4cout << "ElementWiseData::UpdateData Init finished"<<G4endl;
108 
109  theBuffer = theIsotopeWiseData[index].MakeElasticData();
110 // G4cout << "ElementWiseData::UpdateData MakeElasticData finished: "
111 // <<theBuffer->GetVectorLength()<<G4endl;
112  Harmonise(theElasticData, theBuffer);
113 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
114 // <<theElasticData->GetVectorLength()<<G4endl;
115  delete theBuffer;
116 
117  theBuffer = theIsotopeWiseData[index].MakeInelasticData();
118  // G4cout << "ElementWiseData::UpdateData MakeInelasticData finished: "
119 // <<theBuffer->GetVectorLength()<<G4endl;
120  Harmonise(theInelasticData, theBuffer);
121 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
122 // <<theInelasticData->GetVectorLength()<<G4endl;
123  delete theBuffer;
124 
125  theBuffer = theIsotopeWiseData[index].MakeCaptureData();
126 // G4cout << "ElementWiseData::UpdateData MakeCaptureData finished: "
127 // <<theBuffer->GetVectorLength()<<G4endl;
128  Harmonise(theCaptureData, theBuffer);
129 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
130 // <<theCaptureData->GetVectorLength()<<G4endl;
131  delete theBuffer;
132 
133  theBuffer = theIsotopeWiseData[index].MakeFissionData();
134 // G4cout << "ElementWiseData::UpdateData MakeFissionData finished: "
135 // <<theBuffer->GetVectorLength()<<G4endl;
136  Harmonise(theFissionData, theBuffer);
137 // G4cout << "ElementWiseData::UpdateData Harmonise finished: "
138 // <<theFissionData->GetVectorLength()<<G4endl;
139  delete theBuffer;
140 
141 // G4cout << "ElementWiseData::UpdateData finished"<endl;
142  }
143 
144  void G4ParticleHPElementData::Harmonise(G4ParticleHPVector *& theStore, G4ParticleHPVector * theNew)
145  {
146  if(theNew == 0) { return; }
147  G4int s_tmp = 0, n=0, m_tmp=0;
148  G4ParticleHPVector * theMerge = new G4ParticleHPVector(theStore->GetVectorLength());
149 // G4cout << "Harmonise 1: "<<theStore->GetEnergy(s_tmp)<<" "<<theNew->GetEnergy(0)<<G4endl;
150  while ( theStore->GetEnergy(s_tmp)<theNew->GetEnergy(0)&&s_tmp<theStore->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi
151  {
152  theMerge->SetData(m_tmp++, theStore->GetEnergy(s_tmp), theStore->GetXsec(s_tmp));
153  s_tmp++;
154  }
155  G4ParticleHPVector *active = theStore;
156  G4ParticleHPVector * passive = theNew;
157  G4ParticleHPVector * tmp;
158  G4int a = s_tmp, p = n, t;
159 // G4cout << "Harmonise 2: "<<active->GetVectorLength()<<" "<<passive->GetVectorLength()<<G4endl;
160  while (a<active->GetVectorLength()&&p<passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
161  {
162  if(active->GetEnergy(a) <= passive->GetEnergy(p))
163  {
164  theMerge->SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
165  G4double x = theMerge->GetEnergy(m_tmp);
166  G4double y = std::max(0., passive->GetXsec(x));
167  theMerge->SetData(m_tmp, x, theMerge->GetXsec(m_tmp)+y);
168  m_tmp++;
169  a++;
170  } else {
171 // G4cout << "swapping in Harmonise"<<G4endl;
172  tmp = active; t=a;
173  active = passive; a=p;
174  passive = tmp; p=t;
175  }
176  }
177 // G4cout << "Harmonise 3: "<< a <<" "<<active->GetVectorLength()<<" "<<m<<G4endl;
178  while (a!=active->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
179  {
180  theMerge->SetData(m_tmp++, active->GetEnergy(a), active->GetXsec(a));
181  a++;
182  }
183 // G4cout << "Harmonise 4: "<< p <<" "<<passive->GetVectorLength()<<" "<<m<<G4endl;
184  while (p!=passive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
185  {
186  // Modified by T. KOI
187  //theMerge->SetData(m++, passive->GetEnergy(p), passive->GetXsec(p));
188  G4double x = passive->GetEnergy(p);
189  G4double y = std::max(0., active->GetXsec(x));
190  theMerge->SetData(m_tmp++, x, passive->GetXsec(p)+y);
191  p++;
192  }
193 // G4cout <<"Harmonise 5: "<< theMerge->GetVectorLength() << " " << m << G4endl;
194  delete theStore;
195  theStore = theMerge;
196 // G4cout <<"Harmonise 6: "<< theStore->GetVectorLength() << " " << m << G4endl;
197  }
198 
199  G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
200  G4ParticleDefinition * projectile,
201  G4ParticleHPFissionData* theSet,
202  char* dataDirVariable)
203  {
204  if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
205  Init ( theElement, projectile, dataDirVariable );
206  return GetData(theSet);
207  }
208  G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
209  G4ParticleDefinition * projectile,
210  G4ParticleHPCaptureData * theSet,
211  char* dataDirVariable)
212  {
213  if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
214  Init ( theElement, projectile, dataDirVariable );
215  return GetData(theSet);
216  }
217  G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
218  G4ParticleDefinition * projectile,
219  G4ParticleHPElasticData * theSet,
220  char* dataDirVariable)
221  {
222  if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
223  Init ( theElement, projectile, dataDirVariable );
224  return GetData(theSet);
225  }
226  G4ParticleHPVector * G4ParticleHPElementData::MakePhysicsVector(G4Element * theElement,
227  G4ParticleDefinition * projectile,
228  G4ParticleHPInelasticData * theSet,
229  char* dataDirVariable)
230  {
231  if(projectile != G4Neutron::Neutron()) throw G4HadronicException(__FILE__, __LINE__, "not a neutron");
232  Init ( theElement, projectile, dataDirVariable );
233  return GetData(theSet);
234  }