ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPFastLegendre.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPFastLegendre.hh
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 //
27 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
28 //
29 #ifndef G4ParticleHPFastLegendre_h
30 #define G4ParticleHPFastLegendre_h 1
31 
32 #include "globals.hh"
33 
34 class G4ParticleHPFastLegendre
35 {
36  public:
37 
38  G4ParticleHPFastLegendre()
39  {
40  value = new const G4double * [31];
41  value[0] = l0;
42  value[1] = l1;
43  value[2] = l2;
44  value[3] = l3;
45  value[4] = l4;
46  value[5] = l5;
47  value[6] = l6;
48  value[7] = l7;
49  value[8] = l8;
50  value[9] = l9;
51  value[10] = l10;
52  value[11] = l11;
53  value[12] = l12;
54  value[13] = l13;
55  value[14] = l14;
56  value[15] = l15;
57  value[16] = l16;
58  value[17] = l17;
59  value[18] = l18;
60  value[19] = l19;
61  value[20] = l20;
62  value[21] = l21;
63  value[22] = l22;
64  value[23] = l23;
65  value[24] = l24;
66  value[25] = l25;
67  value[26] = l26;
68  value[27] = l27;
69  value[28] = l28;
70  value[29] = l29;
71  value[30] = l30;
72  integral = new const G4double * [31];
73  integral[0] = i0;
74  integral[1] = i1;
75  integral[2] = i2;
76  integral[3] = i3;
77  integral[4] = i4;
78  integral[5] = i5;
79  integral[6] = i6;
80  integral[7] = i7;
81  integral[8] = i8;
82  integral[9] = i9;
83  integral[10] = i10;
84  integral[11] = i11;
85  integral[12] = i12;
86  integral[13] = i13;
87  integral[14] = i14;
88  integral[15] = i15;
89  integral[16] = i16;
90  integral[17] = i17;
91  integral[18] = i18;
92  integral[19] = i19;
93  integral[20] = i20;
94  integral[21] = i21;
95  integral[22] = i22;
96  integral[23] = i23;
97  integral[24] = i24;
98  integral[25] = i25;
99  integral[26] = i26;
100  integral[27] = i27;
101  integral[28] = i28;
102  integral[29] = i29;
103  integral[30] = i30;
104 
105  G4int i;
106  for(i=0;i<31;i++) theNbin[i]=1+200*(i+1);
107  }
108 
109  ~G4ParticleHPFastLegendre()
110  {
111  delete [] value;
112  delete [] integral;
113  }
114 
115  G4double Integrate(G4int l, G4double costh)
116  {
117  if(l>30) return regularIntegrate(l,costh);
118  G4int bin = GetBin(l, costh);
119  G4double y1, y2;
120 // G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate "<<l<<" "<<bin<<G4endl;
121  y1 = integral[l][bin];
122  y2 = integral[l][bin+1];
123 // G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate exit"<<G4endl;
124  return Interpolate(bin, l, y1, y2, costh);
125  }
126 
127  inline G4double Evaluate(G4int l, G4double costh)
128  {
129  if(l>30) return regularEvaluate(l,costh);
130  G4double result;
131  G4int bin = GetBin(l, costh);
132  if(bin != theNbin[l]-1)
133  {
134  G4double y1, y2;
135  y1 = value[l][bin];
136  y2 = value[l][bin+1];
137  result = Interpolate(bin, l, y1, y2, costh);
138  }
139  else
140  {
141  result = value[l][bin];
142  }
143  return result;
144  }
145 
146  private:
147 
148  G4double regularEvaluate( int l , double x );
149  G4double regularIntegrate( int l , double x );
150 
151 
152  inline G4int GetBin(G4int l, G4double costh)
153  {
154  G4int bin=0;
155  bin = G4int( (theNbin[l]-1)*(costh+1)/2. );
156  if(bin == theNbin[l]-1) bin--;
157  return bin;
158  }
159 
160  inline G4double Interpolate(G4int bin, G4int l, G4double y1, G4double y2, G4double x)
161  {
162  G4double slope = 0, off = 0, x2=0, x1mx2;
163  G4int half = (theNbin[l]-1)/2;
164 // x1 = (bin-half)/G4double(half);
165  x2 = (bin+1-half)/G4double(half);
166  x1mx2 = 1./G4double( (theNbin[l]-1)/2 );
167 // slope = (y2-y1)/(x2-x1);
168  slope = (y2-y1)/x1mx2;
169  off = y2-x2*slope;
170  return x*slope+off;
171  }
172 
173  const G4double ** value;
174  const G4double ** integral;
175  G4int theNbin[31];
176  static const G4double l0[201];
177  static const G4double i0[201];
178  static const G4double l1[401];
179  static const G4double i1[401];
180  static const G4double l2[601];
181  static const G4double i2[601];
182  static const G4double l3[801];
183  static const G4double i3[801];
184  static const G4double l4[1001];
185  static const G4double i4[1001];
186  static const G4double l5[1201];
187  static const G4double i5[1201];
188  static const G4double l6[1401];
189  static const G4double i6[1401];
190  static const G4double l7[1601];
191  static const G4double i7[1601];
192  static const G4double l8[1801];
193  static const G4double i8[1801];
194  static const G4double l9[2001];
195  static const G4double i9[2001];
196  static const G4double l10[2201];
197  static const G4double i10[2201];
198  static const G4double l11[2401];
199  static const G4double i11[2401];
200  static const G4double l12[2601];
201  static const G4double i12[2601];
202  static const G4double l13[2801];
203  static const G4double i13[2801];
204  static const G4double l14[3001];
205  static const G4double i14[3001];
206  static const G4double l15[3201];
207  static const G4double i15[3201];
208  static const G4double l16[3401];
209  static const G4double i16[3401];
210  static const G4double l17[3601];
211  static const G4double i17[3601];
212  static const G4double l18[3801];
213  static const G4double i18[3801];
214  static const G4double l19[4001];
215  static const G4double i19[4001];
216  static const G4double l20[4201];
217  static const G4double i20[4201];
218  static const G4double l21[4401];
219  static const G4double i21[4401];
220  static const G4double l22[4601];
221  static const G4double i22[4601];
222  static const G4double l23[4801];
223  static const G4double i23[4801];
224  static const G4double l24[5001];
225  static const G4double i24[5001];
226  static const G4double l25[5201];
227  static const G4double i25[5201];
228  static const G4double l26[5401];
229  static const G4double i26[5401];
230  static const G4double l27[5601];
231  static const G4double i27[5601];
232  static const G4double l28[5801];
233  static const G4double i28[5801];
234  static const G4double l29[6001];
235  static const G4double i29[6001];
236  static const G4double l30[6201];
237  static const G4double i30[6201];
238 };
239 #endif