ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPLabAngularEnergy.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPLabAngularEnergy.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 // 080808 Bug fix in serching mu bin and index for theBuff2b by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 #include "G4Gamma.hh"
39 #include "G4Electron.hh"
40 #include "G4Positron.hh"
41 #include "G4Neutron.hh"
42 #include "G4Proton.hh"
43 #include "G4Deuteron.hh"
44 #include "G4Triton.hh"
45 #include "G4He3.hh"
46 #include "G4Alpha.hh"
47 
48 void G4ParticleHPLabAngularEnergy::Init(std::istream & aDataFile)
49 {
50  aDataFile >> nEnergies;
51  theManager.Init(aDataFile);
53  nCosTh = new G4int[nEnergies];
56  for(G4int i=0; i<nEnergies; i++)
57  {
58  aDataFile >> theEnergies[i];
59  theEnergies[i]*=eV;
60  aDataFile >> nCosTh[i];
61  theSecondManager[i].Init(aDataFile);
62  theData[i] = new G4ParticleHPVector[nCosTh[i]];
63  G4double label;
64  for(G4int ii=0; ii<nCosTh[i]; ii++)
65  {
66  aDataFile >> label;
67  theData[i][ii].SetLabel(label);
68  theData[i][ii].Init(aDataFile, eV);
69  }
70  }
71 }
72 
74 {
75  G4ReactionProduct * result = new G4ReactionProduct;
76  G4int Z = static_cast<G4int>(massCode/1000);
77  G4int A = static_cast<G4int>(massCode-1000*Z);
78 
79  if(massCode==0)
80  {
81  result->SetDefinition(G4Gamma::Gamma());
82  }
83  else if(A==0)
84  {
86  if(Z==1) result->SetDefinition(G4Positron::Positron());
87  }
88  else if(A==1)
89  {
91  if(Z==1) result->SetDefinition(G4Proton::Proton());
92  }
93  else if(A==2)
94  {
96  }
97  else if(A==3)
98  {
99  result->SetDefinition(G4Triton::Triton());
100  if(Z==2) result->SetDefinition(G4He3::He3());
101  }
102  else if(A==4)
103  {
104  result->SetDefinition(G4Alpha::Alpha());
105  if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
106  }
107  else
108  {
109  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPLabAngularEnergy: Unknown ion case 2");
110  }
111 
112  // get theta, E
113  G4double cosTh, secEnergy;
114  G4int i, it(0);
115  // find the energy bin
116  for(i=0; i<nEnergies; i++)
117  {
118  it = i;
119  if ( anEnergy < theEnergies[i] ) break;
120  }
121  //080808
122  //if ( it == 0 || it == nEnergies-1 ) // it marks the energy bin
123  if ( it == 0 ) // it marks the energy bin
124  {
125  G4cout << "080808 Something unexpected is happen in G4ParticleHPLabAngularEnergy " << G4endl;
126  // integrate the prob for each costh, and select theta.
127  G4double * running = new G4double [nCosTh[it]];
128  running[0]=0;
129  for(i=0;i<nCosTh[it]; i++)
130  {
131  if(i!=0) running[i] = running[i-1];
132  running[i]+=theData[it][i].GetIntegral(); // Does interpolated integral.
133  }
134  G4double random = running[nCosTh[it]-1]*G4UniformRand();
135  G4int ith(0);
136  for(i=0;i<nCosTh[it]; i++)
137  {
138  ith = i;
139  if(random<running[i]) break;
140  }
141  //080807
142  //if ( ith == 0 || ith == nCosTh[it]-1 ) //ith marks the angluar bin
143  if ( ith == 0 ) //ith marks the angluar bin
144  {
145  cosTh = theData[it][ith].GetLabel();
146  secEnergy = theData[it][ith].Sample();
148  }
149  else
150  {
151  //080808
152  //G4double x1 = theData[it][ith-1].GetIntegral();
153  //G4double x2 = theData[it][ith].GetIntegral();
154  G4double x1 = running [ ith-1 ];
155  G4double x2 = running [ ith ];
156  G4double x = random;
157  G4double y1 = theData[it][ith-1].GetLabel();
158  G4double y2 = theData[it][ith].GetLabel();
159  cosTh = theInt.Interpolate(theSecondManager[it].GetInverseScheme(ith),
160  x, x1, x2, y1, y2);
161  G4ParticleHPVector theBuff1;
162  theBuff1.SetInterpolationManager(theData[it][ith-1].GetInterpolationManager());
163  G4ParticleHPVector theBuff2;
164  theBuff2.SetInterpolationManager(theData[it][ith].GetInterpolationManager());
165  x1=y1;
166  x2=y2;
167  G4double y, mu;
168  for(i=0;i<theData[it][ith-1].GetVectorLength(); i++)
169  {
170  mu = theData[it][ith-1].GetX(i);
171  y1 = theData[it][ith-1].GetY(i);
172  y2 = theData[it][ith].GetY(mu);
173 
174  y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
175  cosTh, x1,x2,y1,y2);
176  theBuff1.SetData(i, mu, y);
177  }
178  for(i=0;i<theData[it][ith].GetVectorLength(); i++)
179  {
180  mu = theData[it][ith].GetX(i);
181  y1 = theData[it][ith-1].GetY(mu);
182  y2 = theData[it][ith].GetY(i);
183  y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
184  cosTh, x1,x2,y1,y2);
185  theBuff2.SetData(i, mu, y);
186  }
187  G4ParticleHPVector theStore;
188  theStore.Merge(&theBuff1, &theBuff2);
189  secEnergy = theStore.Sample();
190  currentMeanEnergy = theStore.GetMeanX();
191  }
192  delete [] running;
193  }
194  else // this is the small big else.
195  {
196  G4double x, x1, x2, y1, y2, y, tmp, E;
197  // integrate the prob for each costh, and select theta.
198  G4ParticleHPVector run1;
199  run1.SetY(0, 0.);
200  for(i=0;i<nCosTh[it-1]; i++)
201  {
202  if(i!=0) run1.SetY(i, run1.GetY(i-1));
203  run1.SetX(i, theData[it-1][i].GetLabel());
204  run1.SetY(i, run1.GetY(i)+theData[it-1][i].GetIntegral());
205  }
206  G4ParticleHPVector run2;
207  run2.SetY(0, 0.);
208  for(i=0;i<nCosTh[it]; i++)
209  {
210  if(i!=0) run2.SetY(i, run2.GetY(i-1));
211  run2.SetX(i, theData[it][i].GetLabel());
212  run2.SetY(i, run2.GetY(i)+theData[it][i].GetIntegral());
213  }
214  // get the distributions for the correct neutron energy
215  x = anEnergy;
216  x1 = theEnergies[it-1];
217  x2 = theEnergies[it];
218  G4ParticleHPVector thBuff1; // to be interpolated as run1.
220  for(i=0; i<run1.GetVectorLength(); i++)
221  {
222  tmp = run1.GetX(i); //theta
223  y1 = run1.GetY(i); // integral
224  y2 = run2.GetY(tmp);
226  thBuff1.SetData(i, tmp, y);
227  }
228  G4ParticleHPVector thBuff2;
230  for(i=0; i<run2.GetVectorLength(); i++)
231  {
232  tmp = run2.GetX(i); //theta
233  y1 = run1.GetY(tmp); // integral
234  y2 = run2.GetY(i);
235  y = theInt.Lin(x, x1,x2,y1,y2);
236  thBuff2.SetData(i, tmp, y);
237  }
238  G4ParticleHPVector theThVec;
239  theThVec.Merge(&thBuff1 ,&thBuff2); // takes care of interpolation
240  G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
241  -theThVec.GetY(0)) *G4UniformRand();
242  G4int ith(0);
243  for(i=1;i<theThVec.GetVectorLength(); i++)
244  {
245  ith = i;
246  if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
247  }
248  {
249  // calculate theta
250  G4double xx, xx1, xx2, yy1, yy2;
251  xx = random;
252  xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
253  xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
254  yy1 = theThVec.GetX(ith-1); // std::cos(theta)
255  yy2 = theThVec.GetX(ith);
256  cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
257  xx, xx1,xx2,yy1,yy2);
258  }
259  G4int i1(0), i2(0);
260  // get the indixes of the vectors close to theta for low energy
261  // first it-1 !!!! i.e. low in energy
262  for(i=0; i<nCosTh[it-1]; i++)
263  {
264  i1 = i;
265  if(cosTh<theData[it-1][i].GetLabel()) break;
266  }
267  // now get the prob at this energy for the right theta value
268  x = cosTh;
269  x1 = theData[it-1][i1-1].GetLabel();
270  x2 = theData[it-1][i1].GetLabel();
271  G4ParticleHPVector theBuff1a;
272  theBuff1a.SetInterpolationManager(theData[it-1][i1-1].GetInterpolationManager());
273  for(i=0;i<theData[it-1][i1-1].GetVectorLength(); i++)
274  {
275  E = theData[it-1][i1-1].GetX(i);
276  y1 = theData[it-1][i1-1].GetY(i);
277  y2 = theData[it-1][i1].GetY(E);
278  y = theInt.Lin(x, x1,x2,y1,y2);
279  theBuff1a.SetData(i, E, y); // wrong E, right theta.
280  }
281  G4ParticleHPVector theBuff2a;
282  theBuff2a.SetInterpolationManager(theData[it-1][i1].GetInterpolationManager());
283  for(i=0;i<theData[it-1][i1].GetVectorLength(); i++)
284  {
285  E = theData[it-1][i1].GetX(i);
286  y1 = theData[it-1][i1-1].GetY(E);
287  y2 = theData[it-1][i1].GetY(i);
288  y = theInt.Lin(x, x1,x2,y1,y2);
289  theBuff2a.SetData(i, E, y); // wrong E, right theta.
290  }
291  G4ParticleHPVector theStore1;
292  theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
293 
294  // get the indixes of the vectors close to theta for high energy
295  // then it !!!! i.e. high in energy
296  for(i=0; i<nCosTh[it]; i++)
297  {
298  i2 = i;
299  if(cosTh<theData[it][i2].GetLabel()) break;
300  } // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
301  x1 = theData[it][i2-1].GetLabel();
302  x2 = theData[it][i2].GetLabel();
303  G4ParticleHPVector theBuff1b;
304  theBuff1b.SetInterpolationManager(theData[it][i2-1].GetInterpolationManager());
305  for(i=0;i<theData[it][i2-1].GetVectorLength(); i++)
306  {
307  E = theData[it][i2-1].GetX(i);
308  y1 = theData[it][i2-1].GetY(i);
309  y2 = theData[it][i2].GetY(E);
310  y = theInt.Lin(x, x1,x2,y1,y2);
311  theBuff1b.SetData(i, E, y); // wrong E, right theta.
312  }
313  G4ParticleHPVector theBuff2b;
314  theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
315  //080808 i1 -> i2
316  //for(i=0;i<theData[it][i1].GetVectorLength(); i++)
317  for(i=0;i<theData[it][i2].GetVectorLength(); i++)
318  {
319  //E = theData[it][i1].GetX(i);
320  //y1 = theData[it][i1-1].GetY(E);
321  //y2 = theData[it][i1].GetY(i);
322  E = theData[it][i2].GetX(i);
323  y1 = theData[it][i2-1].GetY(E);
324  y2 = theData[it][i2].GetY(i);
325  y = theInt.Lin(x, x1,x2,y1,y2);
326  theBuff2b.SetData(i, E, y); // wrong E, right theta.
327  }
328  G4ParticleHPVector theStore2;
329  theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
330  // now get to the right energy.
331 
332  x = anEnergy;
333  x1 = theEnergies[it-1];
334  x2 = theEnergies[it];
335  G4ParticleHPVector theOne1;
337  for(i=0; i<theStore1.GetVectorLength(); i++)
338  {
339  E = theStore1.GetX(i);
340  y1 = theStore1.GetY(i);
341  y2 = theStore2.GetY(E);
343  theOne1.SetData(i, E, y); // both correct
344  }
345  G4ParticleHPVector theOne2;
347  for(i=0; i<theStore2.GetVectorLength(); i++)
348  {
349  E = theStore2.GetX(i);
350  y1 = theStore1.GetY(E);
351  y2 = theStore2.GetY(i);
353  theOne2.SetData(i, E, y); // both correct
354  }
355  G4ParticleHPVector theOne;
356  theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
357 
358  secEnergy = theOne.Sample();
359  currentMeanEnergy = theOne.GetMeanX();
360  }
361 
362 // now do random direction in phi, and fill the result.
363 
364  result->SetKineticEnergy(secEnergy);
365 
367  G4double theta = std::acos(cosTh);
368  G4double sinth = std::sin(theta);
369  G4double mtot = result->GetTotalMomentum();
370  G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
371  result->SetMomentum(tempVector);
372 
373  return result;
374 }