Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Author: Alexei Sytov 27 28 /// \file G4ChannelingFastSimInterpolation.cc 29 /// \brief Implementation of the G4ChannelingF 30 31 #include "G4ChannelingFastSimInterpolation.hh" 32 #include "G4SystemOfUnits.hh" 33 34 G4ChannelingFastSimInterpolation::G4Channeling 35 36 37 38 39 { 40 fDx = dx0; 41 fDy = dy0; 42 nPointsx = nPointsx0; 43 nPointsy = nPointsy0; 44 fStepi = fDx/nPointsx; 45 fStepi2= fStepi*fStepi; 46 iModel = iModel0; 47 48 //resize vectors for interpolation coeffic 49 if(iModel==1) 50 { 51 fAI.resize(nPointsx+1); 52 fBI.resize(nPointsx+1); 53 fCI.resize(nPointsx+1); 54 fDI.resize(nPointsx+1); 55 } 56 else if(iModel==2) 57 { 58 fStepj = fDy/nPointsy; 59 fAI3D.resize(nPointsx+1, std::vector<G 60 fBI3D.resize(nPointsx+1, std::vector<G 61 fCI3D.resize(nPointsx+1, std::vector<G 62 fAI3D3.resize(nPointsx+1, std::vector< 63 fBI3D3.resize(nPointsx+1, std::vector< 64 fCI3D3.resize(nPointsx+1, std::vector< 65 66 for(G4int i=0; i<=nPointsx; i++) 67 { 68 fCI3D[i][nPointsy]=0.; 69 } 70 } 71 } 72 73 //....oooOO0OOooo........oooOO0OOooo........oo 74 75 G4double G4ChannelingFastSimInterpolation::Get 76 G4double SplineF=0.; 77 if(iModel==1) 78 { 79 SplineF = Spline1D(xx); 80 } 81 else if(iModel==2) 82 { 83 SplineF = Spline2D(xx, yy); 84 } 85 return SplineF; 86 } 87 88 //....oooOO0OOooo........oooOO0OOooo........oo 89 90 void G4ChannelingFastSimInterpolation::SetCoef 91 92 93 94 95 fAI[i] = AI0; 96 fBI[i] = BI0/cm; 97 fCI[i] = CI0/cm2; 98 fDI[i] = DI0/cm3; 99 } 100 101 //....oooOO0OOooo........oooOO0OOooo........oo 102 103 void G4ChannelingFastSimInterpolation::SetCoef 104 105 106 107 108 109 if (k==0) 110 { 111 fAI3D[i][j] = AI3D0/fStepj/fStepi/6.; 112 fBI3D[i][j] = BI3D0/fStepj/fStepi/6.; 113 fCI3D[i][j] = CI3D0/fStepj/fStepi/6./c 114 } 115 else if (k==1) 116 { 117 fAI3D3[i][j] = AI3D0/fStepj/fStepi/6./ 118 fBI3D3[i][j] = BI3D0/fStepj/fStepi/6./ 119 fCI3D3[i][j] = CI3D0/fStepj/fStepi/6./ 120 } 121 } 122 123 //....oooOO0OOooo........oooOO0OOooo........oo 124 125 G4double G4ChannelingFastSimInterpolation::Spl 126 127 { 128 G4double x1 = xx; 129 130 //if a particle escapes the interpolation 131 if (x1<0.) 132 { 133 x1 += fDx; 134 } 135 else if (x1>=fDx) 136 { 137 x1 -= fDx; 138 } 139 140 // calculation of interpolation function 141 G4int mmx = std::floor(x1/fStepi); 142 x1 -= (mmx+1)*fStepi; 143 G4double Spline3 = fAI[mmx]+x1*(fBI[mmx]+ 144 return Spline3; 145 } 146 147 //....oooOO0OOooo........oooOO0OOooo........oo 148 149 G4double G4ChannelingFastSimInterpolation:: 150 151 { 152 G4double x1 = xx; 153 G4double y1 = yy; 154 155 //if a particle escapes the interpolation 156 if (x1<0.) 157 { 158 x1 += fDx; 159 } 160 else if (x1>=fDx) 161 { 162 x1 -= fDx; 163 } 164 if (y1<0.) 165 { 166 y1 += fDy; 167 } 168 else if (y1>=fDy) 169 { 170 y1 -= fDy; 171 } 172 173 // calculation of interpolation function 174 G4int mmx = std::floor(x1/fStepi); 175 G4int mmy = std::floor(y1/fStepj); 176 177 G4double tt1 = x1-mmx*fStepi; 178 G4double tt2 = fStepi-tt1; 179 G4double tt13 = tt1*tt1*tt1; 180 G4double tt23 = tt2*tt2*tt2; 181 182 G4double tt1y = y1-mmy*fStepj; 183 G4double tt2y = fStepj-tt1y; 184 G4double tt1y3 = tt1y*tt1y*tt1y; 185 G4double tt2y3 = tt2y*tt2y*tt2y; 186 187 G4double spl3dxx1 = fCI3D3[mmx ][mmy]*tt 188 fAI3D3[mmx ][mmy]*tt 189 G4double spl3dxx2 = fCI3D3[mmx+1][mmy]*tt 190 fAI3D3[mmx+1][mmy]*tt 191 G4double spl3d1 = fCI3D[ mmx ][mmy]*tt 192 fAI3D[ mmx ][mmy]*tt 193 G4double spl3d2 = fCI3D[ mmx+1][mmy]*tt 194 fAI3D[ mmx+1][mmy]*tt 195 196 G4double spl3d = spl3dxx1*tt23 + spl3dxx2 197 (spl3d1*6.-spl3dxx1*fStepi2)*tt2 198 199 return spl3d; 200 } 201 202 203