| Geant4 Cross Reference |
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 // Author: Alexei Sytov
27
28 /// \file G4ChannelingFastSimInterpolation.cc
29 /// \brief Implementation of the G4ChannelingFastSimInterpolation class
30
31 #include "G4ChannelingFastSimInterpolation.hh"
32 #include "G4SystemOfUnits.hh"
33
34 G4ChannelingFastSimInterpolation::G4ChannelingFastSimInterpolation(G4double dx0,
35 G4double dy0,
36 G4int nPointsx0,
37 G4int nPointsy0,
38 G4int iModel0)
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 coefficients
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<G4double>(nPointsy+1));
60 fBI3D.resize(nPointsx+1, std::vector<G4double>(nPointsy+1));
61 fCI3D.resize(nPointsx+1, std::vector<G4double>(nPointsy+1));
62 fAI3D3.resize(nPointsx+1, std::vector<G4double>(nPointsy+1));
63 fBI3D3.resize(nPointsx+1, std::vector<G4double>(nPointsy+1));
64 fCI3D3.resize(nPointsx+1, std::vector<G4double>(nPointsy+1));
65
66 for(G4int i=0; i<=nPointsx; i++)
67 {
68 fCI3D[i][nPointsy]=0.;
69 }
70 }
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75 G4double G4ChannelingFastSimInterpolation::GetIF(G4double xx, G4double yy){
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........oooOO0OOooo........oooOO0OOooo....
89
90 void G4ChannelingFastSimInterpolation::SetCoefficients1D(G4double AI0,
91 G4double BI0,
92 G4double CI0,
93 G4double DI0,
94 G4int i){
95 fAI[i] = AI0;
96 fBI[i] = BI0/cm;
97 fCI[i] = CI0/cm2;
98 fDI[i] = DI0/cm3;
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103 void G4ChannelingFastSimInterpolation::SetCoefficients2D(G4double AI3D0,
104 G4double BI3D0,
105 G4double CI3D0,
106 G4int i,
107 G4int j,
108 G4int k){
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./cm2;
114 }
115 else if (k==1)
116 {
117 fAI3D3[i][j] = AI3D0/fStepj/fStepi/6./cm2;
118 fBI3D3[i][j] = BI3D0/fStepj/fStepi/6./cm2;
119 fCI3D3[i][j] = CI3D0/fStepj/fStepi/6./cm2/cm2;
120 }
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125 G4double G4ChannelingFastSimInterpolation::Spline1D(G4double xx) //cubic spline of
126 //1-variable function
127 {
128 G4double x1 = xx;
129
130 //if a particle escapes the interpolation area
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]+x1*(fCI[mmx]+fDI[mmx]*x1));
144 return Spline3;
145 }
146
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
149 G4double G4ChannelingFastSimInterpolation::Spline2D(G4double xx, G4double yy)//cubic
150 //spline of 2-variable function
151 {
152 G4double x1 = xx;
153 G4double y1 = yy;
154
155 //if a particle escapes the interpolation area
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]*tt2y3 + fCI3D3[mmx ][mmy+1]*tt1y3 +
188 fAI3D3[mmx ][mmy]*tt2y + fBI3D3[mmx ][mmy]*tt1y;
189 G4double spl3dxx2 = fCI3D3[mmx+1][mmy]*tt2y3 + fCI3D3[mmx+1][mmy+1]*tt1y3 +
190 fAI3D3[mmx+1][mmy]*tt2y + fBI3D3[mmx+1][mmy]*tt1y;
191 G4double spl3d1 = fCI3D[ mmx ][mmy]*tt2y3 + fCI3D[mmx ][mmy+1]*tt1y3 +
192 fAI3D[ mmx ][mmy]*tt2y + fBI3D[mmx ][mmy]*tt1y;
193 G4double spl3d2 = fCI3D[ mmx+1][mmy]*tt2y3 + fCI3D[mmx+1][mmy+1]*tt1y3 +
194 fAI3D[ mmx+1][mmy]*tt2y + fBI3D[mmx+1][mmy]*tt1y;
195
196 G4double spl3d = spl3dxx1*tt23 + spl3dxx2*tt13 +
197 (spl3d1*6.-spl3dxx1*fStepi2)*tt2 + (spl3d2*6.-spl3dxx2*fStepi2)*tt1;
198
199 return spl3d;
200 }
201
202
203