ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPChannel.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 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 071031 bug fix T. Koi on behalf of A. Howard
32 // 081203 bug fix in Register method by T. Koi
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 //
36 // June-2019 - E. Mendoza --> Modification to allow using an incomplete data library if the G4NEUTRONHP_SKIP_MISSING_ISOTOPES environmental flag is defined. The missing XS are set to 0.
37 
38 
39 #include <stdlib.h>
40 
41 #include "G4ParticleHPChannel.hh"
42 #include "globals.hh"
43 #include "G4SystemOfUnits.hh"
45 #include "G4HadTmpUtil.hh"
46 
47 #include "G4ParticleHPManager.hh"
49 
51  {
52  return std::max(0., theChannelData->GetXsec(energy));
53  }
54 
56  {
57  return theIsotopeWiseData[isoNumber].GetXsec(energy);
58  }
59 
61  {
62  return theFinalStates[isoNumber]->GetXsec(energy);
63  }
64 
66  Init(G4Element * anElement, const G4String dirName, const G4String aFSType)
67  {
68  theFSType = aFSType;
69  Init(anElement, dirName);
70  }
71 
72  void G4ParticleHPChannel::Init(G4Element * anElement, const G4String dirName)
73  {
74  theDir = dirName;
75  theElement = anElement;
76  }
77 
79  {
80  registerCount++;
82 
83  Z = Z-registerCount;
84  if ( registerCount > 5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material"); // for Elastic, Capture, Fission case
85  if ( Z < 1 ) return false;
86 /*
87  if(registerCount<5)
88  {
89  Z = Z-registerCount;
90  }
91 */
92  //if(Z=theElement->GetZ()-5) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
93  // Bug fix by TK on behalf of AH
94  //if ( Z <=theElement->GetZ()-5 ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this material");
95  G4int count = 0;
96  if(registerCount==0) count = theElement->GetNumberOfIsotopes();
97  if(count == 0||registerCount!=0) count +=
99  niso = count;
100  delete [] theIsotopeWiseData;
102  delete [] active;
103  active = new G4bool[niso];
104 
105  delete [] theFinalStates;
107  delete theChannelData;
109  for(G4int i=0; i<niso; i++)
110  {
111  theFinalStates[i] = theFS->New();
113  }
114  count = 0;
115  G4int nIsos = niso;
116  if(theElement->GetNumberOfIsotopes()!=0&&registerCount==0)
117  {
118  for (G4int i1=0; i1<nIsos; i1++)
119  {
120  // G4cout <<" Init: normal case"<<G4endl;
121  G4int A = theElement->GetIsotope(i1)->GetN();
122  G4int M = theElement->GetIsotope(i1)->Getm();
124  //theFinalStates[i1]->SetA_Z(A, Z);
125  //UpdateData(A, Z, count++, frac);
126  theFinalStates[i1]->SetA_Z(A, Z, M);
127  UpdateData(A, Z, M, count++, frac, theProjectile);
128  }
129  } else {
130  //G4cout <<" Init: mean case: "
131  // <<theStableOnes.GetNumberOfIsotopes(Z)<<" "
132  // <<Z<<" "<<theElement
133  // << G4endl;
135  for(G4int i1=0;
137  i1++)
138  {
140  G4double frac = theStableOnes.GetAbundance(first+i1);
141  theFinalStates[i1]->SetA_Z(A, Z);
142  UpdateData(A, Z, count++, frac, theProjectile);
143  }
144  }
145  G4bool result = HasDataInAnyFinalState();
146 
147  //To avoid issuing hash by worker threads
148  if ( result ) theChannelData->Hash();
149 
150  return result;
151  }
152 
153  //void G4ParticleHPChannel::UpdateData(G4int A, G4int Z, G4int index, G4double abundance)
155  {
156  // Initialze the G4FissionFragment generator for this isomer if needed
158  {
160  }
161 
162  theFinalStates[index]->Init(A, Z, M, theDir, theFSType, projectile);
163  if(!theFinalStates[index]->HasAnyData()) return; // nothing there for exactly this isotope.
164 
165  // the above has put the X-sec into the FS
166  theBuffer = 0;
167  if(theFinalStates[index]->HasXsec())
168  {
169  theBuffer = theFinalStates[index]->GetXsec();
170  theBuffer->Times(abundance/100.);
172  }
173  else // get data from CrossSection directory
174  {
175  G4String tString = "/CrossSection";
176  //active[index] = theIsotopeWiseData[index].Init(A, Z, abundance, theDir, tString);
177  active[index] = theIsotopeWiseData[index].Init(A, Z, M, abundance, theDir, tString);
178  if(active[index]) theBuffer = theIsotopeWiseData[index].MakeChannelData();
179  }
181  }
182 
184  {
185  G4int s_tmp = 0, n=0, m_tmp=0;
186  G4ParticleHPVector * theMerge = new G4ParticleHPVector;
187  G4ParticleHPVector * anActive = theStore;
188  G4ParticleHPVector * aPassive = theNew;
190  G4int a = s_tmp, p = n, t;
191  while (a<anActive->GetVectorLength()&&p<aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
192  {
193  if(anActive->GetEnergy(a) <= aPassive->GetEnergy(p))
194  {
195  G4double xa = anActive->GetEnergy(a);
196  theMerge->SetData(m_tmp, xa, anActive->GetXsec(a)+std::max(0., aPassive->GetXsec(xa)) );
197  m_tmp++;
198  a++;
199  G4double xp = aPassive->GetEnergy(p);
200  if( std::abs(std::abs(xp-xa)/xa)<0.001 )
201  {
202  p++;
203  }
204  } else {
205  tmp = anActive; t=a;
206  anActive = aPassive; a=p;
207  aPassive = tmp; p=t;
208  }
209  }
210  while (a!=anActive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
211  {
212  theMerge->SetData(m_tmp++, anActive->GetEnergy(a), anActive->GetXsec(a));
213  a++;
214  }
215  while (p!=aPassive->GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
216  {
217  if(std::abs(theMerge->GetEnergy(std::max(0,m_tmp-1))-aPassive->GetEnergy(p))/aPassive->GetEnergy(p)>0.001)
218  theMerge->SetData(m_tmp++, aPassive->GetEnergy(p), aPassive->GetXsec(p));
219  p++;
220  }
221  delete theStore;
222  theStore = theMerge;
223  }
224 
226 
228  ApplyYourself(const G4HadProjectile & theTrack, G4int anIsotope)
229  {
230 // G4cout << "G4ParticleHPChannel::ApplyYourself+"<<niso<<G4endl;
231  if ( anIsotope != -1 && anIsotope != -2 )
232  {
233  //Inelastic Case
234  //G4cout << "G4ParticleHPChannel Inelastic Case"
235  //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
238  return theFinalStates[anIsotope]->ApplyYourself(theTrack);
239  }
240  G4double sum=0;
241  G4int it=0;
242  G4double * xsec = new G4double[niso];
243  G4ParticleHPThermalBoost aThermalE;
244  for (G4int i=0; i<niso; i++)
245  {
246  if(theFinalStates[i]->HasAnyData())
247  {
248  xsec[i] = theIsotopeWiseData[i].GetXsec(aThermalE.GetThermalEnergy(theTrack,
249  theFinalStates[i]->GetN(),
250  theFinalStates[i]->GetZ(),
251  theTrack.GetMaterial()->GetTemperature()));
252  sum += xsec[i];
253  }
254  else
255  {
256  xsec[i]=0;
257  }
258  }
259  if(sum == 0)
260  {
261 // G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize+"<<G4endl;
262 // G4cout << "G4ParticleHPChannel::ApplyYourself theFinalState->Initialize-"<<G4endl;
263  it = static_cast<G4int>(niso*G4UniformRand());
264  }
265  else
266  {
267 // G4cout << "Are we still here? "<<sum<<G4endl;
268 // G4cout << "TESTHP 23 NISO="<<niso<<G4endl;
269  G4double random = G4UniformRand();
270  G4double running=0;
271 // G4cout << "G4ParticleHPChannel::ApplyYourself Done the sum"<<niso<<G4endl;
272 // G4cout << "TESTHP 24 NISO="<<niso<<G4endl;
273  for (G4int ix=0; ix<niso; ix++)
274  {
275  running += xsec[ix];
276  //if(random<=running/sum)
277  if( sum == 0 || random <= running/sum )
278  {
279  it = ix;
280  break;
281  }
282  }
283  if(it==niso) it--;
284  }
285  delete [] xsec;
286  G4HadFinalState * theFinalState=0;
287  const G4int A = (G4int)this->GetN(it);
288  const G4int Z = (G4int)this->GetZ(it);
289  const G4int M = (G4int)this->GetM(it);
290 
291  //-2:Marker for Fission
292  if(wendtFissionGenerator&&anIsotope==-2)
293  {
294  theFinalState = wendtFissionGenerator->ApplyYourself(theTrack, Z, A);
295  }
296 
297  // Use the standard procedure if the G4FissionFragmentGenerator model fails
298  if (!theFinalState)
299  {
300 
301  G4int icounter=0;
302  G4int icounter_max=1024;
303  while(theFinalState==0) // Loop checking, 11.05.2015, T. Koi
304  {
305  icounter++;
306  if ( icounter > icounter_max ) {
307  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
308  break;
309  }
310 // G4cout << "TESTHP 24 it="<<it<<G4endl;
311  theFinalState = theFinalStates[it]->ApplyYourself(theTrack);
312  }
313  }
314 
315  //G4cout <<"THE IMPORTANT RETURN"<<G4endl;
316  //G4cout << "TK G4ParticleHPChannel Elastic, Capture and Fission Cases "
317  //<< " Z= " << this->GetZ(it) << " A = " << this->GetN(it) << G4endl;
321 
322  return theFinalState;
323  }
324 
325 
327 
328  G4cout<<" Element: "<<theElement->GetName()<<G4endl;
329  G4cout<<" Directory name: "<<theDir<<G4endl;
330  G4cout<<" FS name: "<<theFSType<<G4endl;
331  G4cout<<" Number of Isotopes: "<<niso<<G4endl;
332  G4cout<<" Have cross sections: "<<G4endl;
333  for(int i=0;i<niso;i++){
334  G4cout<<theFinalStates[i]->HasXsec()<<" ";
335  }
336  G4cout<<G4endl;
337  if(theChannelData){
338  G4cout<<" Cross Section (total for this channel):"<<G4endl;
340  G4cout<<np<<G4endl;
341  for(int i=0;i<np;i++){
343  }
344  }
345 
346 }
347 
348 
349