include/Geant4/G4ParticleHPFastLegendre.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
0028 //
0029 #ifndef G4ParticleHPFastLegendre_h
0030 #define G4ParticleHPFastLegendre_h 1
0031
0032 #include "globals.hh"
0033
0034 class G4ParticleHPFastLegendre
0035 {
0036   public:
0037     G4ParticleHPFastLegendre()
0038     {
0039       value = new const G4double*[31];
0040       value[0] = l0;
0041       value[1] = l1;
0042       value[2] = l2;
0043       value[3] = l3;
0044       value[4] = l4;
0045       value[5] = l5;
0046       value[6] = l6;
0047       value[7] = l7;
0048       value[8] = l8;
0049       value[9] = l9;
0050       value[10] = l10;
0051       value[11] = l11;
0052       value[12] = l12;
0053       value[13] = l13;
0054       value[14] = l14;
0055       value[15] = l15;
0056       value[16] = l16;
0057       value[17] = l17;
0058       value[18] = l18;
0059       value[19] = l19;
0060       value[20] = l20;
0061       value[21] = l21;
0062       value[22] = l22;
0063       value[23] = l23;
0064       value[24] = l24;
0065       value[25] = l25;
0066       value[26] = l26;
0067       value[27] = l27;
0068       value[28] = l28;
0069       value[29] = l29;
0070       value[30] = l30;
0071       integral = new const G4double*[31];
0072       integral[0] = i0;
0073       integral[1] = i1;
0074       integral[2] = i2;
0075       integral[3] = i3;
0076       integral[4] = i4;
0077       integral[5] = i5;
0078       integral[6] = i6;
0079       integral[7] = i7;
0080       integral[8] = i8;
0081       integral[9] = i9;
0082       integral[10] = i10;
0083       integral[11] = i11;
0084       integral[12] = i12;
0085       integral[13] = i13;
0086       integral[14] = i14;
0087       integral[15] = i15;
0088       integral[16] = i16;
0089       integral[17] = i17;
0090       integral[18] = i18;
0091       integral[19] = i19;
0092       integral[20] = i20;
0093       integral[21] = i21;
0094       integral[22] = i22;
0095       integral[23] = i23;
0096       integral[24] = i24;
0097       integral[25] = i25;
0098       integral[26] = i26;
0099       integral[27] = i27;
0100       integral[28] = i28;
0101       integral[29] = i29;
0102       integral[30] = i30;
0103
0104       G4int i;
0105       for (i = 0; i < 31; i++)
0106         theNbin[i] = 1 + 200 * (i + 1);
0107     }
0108
0109     ~G4ParticleHPFastLegendre()
0110     {
0111       delete[] value;
0112       delete[] integral;
0113     }
0114
0115     G4double Integrate(G4int l, G4double costh)
0116     {
0117       if (l > 30) return regularIntegrate(l, costh);
0118       G4int bin = GetBin(l, costh);
0119       G4double y1, y2;
0120       //    G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate "<<l<<" "<<bin<<G4endl;
0121       y1 = integral[l][bin];
0122       y2 = integral[l][bin + 1];
0123       //    G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate exit"<<G4endl;
0124       return Interpolate(bin, l, y1, y2, costh);
0125     }
0126
0127     inline G4double Evaluate(G4int l, G4double costh)
0128     {
0129       if (l > 30) return regularEvaluate(l, costh);
0130       G4double result;
0131       G4int bin = GetBin(l, costh);
0132       if (bin != theNbin[l] - 1) {
0133         G4double y1, y2;
0134         y1 = value[l][bin];
0135         y2 = value[l][bin + 1];
0136         result = Interpolate(bin, l, y1, y2, costh);
0137       }
0138       else {
0139         result = value[l][bin];
0140       }
0141       return result;
0142     }
0143
0144   private:
0145     G4double regularEvaluate(int l, double x);
0146     G4double regularIntegrate(int l, double x);
0147
0148     inline G4int GetBin(G4int l, G4double costh)
0149     {
0150       G4int bin = 0;
0151       bin = G4int((theNbin[l] - 1) * (costh + 1) / 2.);
0152       if (bin == theNbin[l] - 1) bin--;
0153       return bin;
0154     }
0155
0156     inline G4double Interpolate(G4int bin, G4int l, G4double y1, G4double y2, G4double x)
0157     {
0158       G4double slope = 0, off = 0, x2 = 0, x1mx2;
0159       G4int half = (theNbin[l] - 1) / 2;
0160       //    x1 = (bin-half)/G4double(half);
0161       x2 = (bin + 1 - half) / G4double(half);
0162       x1mx2 = 1. / G4double((theNbin[l] - 1) / 2);
0163       //    slope = (y2-y1)/(x2-x1);
0164       slope = (y2 - y1) / x1mx2;
0165       off = y2 - x2 * slope;
0166       return x * slope + off;
0167     }
0168
0169     const G4double** value;
0170     const G4double** integral;
0171     G4int theNbin[31];
0172     static const G4double l0[201];
0173     static const G4double i0[201];
0174     static const G4double l1[401];
0175     static const G4double i1[401];
0176     static const G4double l2[601];
0177     static const G4double i2[601];
0178     static const G4double l3[801];
0179     static const G4double i3[801];
0180     static const G4double l4[1001];
0181     static const G4double i4[1001];
0182     static const G4double l5[1201];
0183     static const G4double i5[1201];
0184     static const G4double l6[1401];
0185     static const G4double i6[1401];
0186     static const G4double l7[1601];
0187     static const G4double i7[1601];
0188     static const G4double l8[1801];
0189     static const G4double i8[1801];
0190     static const G4double l9[2001];
0191     static const G4double i9[2001];
0192     static const G4double l10[2201];
0193     static const G4double i10[2201];
0194     static const G4double l11[2401];
0195     static const G4double i11[2401];
0196     static const G4double l12[2601];
0197     static const G4double i12[2601];
0198     static const G4double l13[2801];
0199     static const G4double i13[2801];
0200     static const G4double l14[3001];
0201     static const G4double i14[3001];
0202     static const G4double l15[3201];
0203     static const G4double i15[3201];
0204     static const G4double l16[3401];
0205     static const G4double i16[3401];
0206     static const G4double l17[3601];
0207     static const G4double i17[3601];
0208     static const G4double l18[3801];
0209     static const G4double i18[3801];
0210     static const G4double l19[4001];
0211     static const G4double i19[4001];
0212     static const G4double l20[4201];
0213     static const G4double i20[4201];
0214     static const G4double l21[4401];
0215     static const G4double i21[4401];
0216     static const G4double l22[4601];
0217     static const G4double i22[4601];
0218     static const G4double l23[4801];
0219     static const G4double i23[4801];
0220     static const G4double l24[5001];
0221     static const G4double i24[5001];
0222     static const G4double l25[5201];
0223     static const G4double i25[5201];
0224     static const G4double l26[5401];
0225     static const G4double i26[5401];
0226     static const G4double l27[5601];
0227     static const G4double i27[5601];
0228     static const G4double l28[5801];
0229     static const G4double i28[5801];
0230     static const G4double l29[6001];
0231     static const G4double i29[6001];
0232     static const G4double l30[6201];
0233     static const G4double i30[6201];
0234 };
0235 #endif