| 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 // AMR Simplification /4
27 // read Diff XSection & Interpolate
28
29 #include "globals.hh"
30 #include "Randomize.hh"
31
32 #include "CLHEP/Units/PhysicalConstants.h"
33
34 #include "G4LEPTSDiffXS.hh"
35 #include "G4Exp.hh"
36
37 G4LEPTSDiffXS::G4LEPTSDiffXS(const G4String& file) {
38 fileName = file;
39
40 readDXS();
41 BuildCDXS();
42 //BuildCDXS(1.0, 0.5);
43 NormalizeCDXS();
44 InterpolateCDXS();
45 }
46
47
48
49 //DXS y KT
50 void G4LEPTSDiffXS::readDXS( ) {
51
52 FILE *fp;
53 G4float data, data2;
54
55 if ((fp=fopen(fileName.c_str(), "r"))==nullptr){
56 //G4cout << "Error reading " << fileName << G4endl;
57 NumEn = 0;
58 bFileFound = false;
59 return;
60 }
61
62 bFileFound = true;
63
64 //G4cout << "Reading2 " << fileName << G4endl;
65
66 //NumAng = 181;
67 (void) fscanf(fp, "%d %d %s", &NumAng, &NumEn, DXSTypeName);
68 if( strcmp(DXSTypeName, "KTC") == 0 ) DXSType = 2; // read DXS & calculate KT
69 else if( strcmp(DXSTypeName, "KT") == 0 ) DXSType = 1; // read DXS & KT
70 else DXSType = 0;
71
72 // if( verboseLevel >= 1 ) G4cout << "Read DXS (" << fileName << ")\t NEg " << NumEn << " NAng " << NumAng
73 // << "DXSType " << DXSTypeName << " " << DXSType << G4endl;
74
75 for (G4int eBin=1; eBin<=NumEn; eBin++){
76 (void) fscanf(fp,"%f ",&data);
77 Eb[eBin] = (G4double)data;
78 }
79
80
81 //for (aBin=1;aBin<NumAng;aBin++){
82
83 if(DXSType==1) {
84 G4cout << "DXSTYpe 1" << G4endl;
85 for (G4int aBin=0;aBin<NumAng;aBin++){
86 (void) fscanf(fp,"%f ",&data);
87 DXS[0][aBin]=(G4double)data;
88 for (G4int eBin=1;eBin<=NumEn;eBin++){
89 (void) fscanf(fp,"%f %f ",&data2, &data);
90 DXS[eBin][aBin]=(G4double)data;
91 KT[eBin][aBin]=(G4double)data2;
92 }
93 }
94 }
95 else {
96 for(G4int aBin=0; aBin<NumAng; aBin++){
97 for(G4int eBin=0; eBin<=NumEn; eBin++){
98 (void) fscanf(fp,"%f ",&data);
99 DXS[eBin][aBin] = (G4double)data;
100 }
101 }
102 for(G4int aBin=0; aBin<NumAng; aBin++){
103 for(G4int eBin=1; eBin<=NumEn; eBin++){
104 G4double A = DXS[0][aBin]; // Angle
105 G4double E = Eb[eBin]; // Energy
106 G4double p = std::sqrt(std::pow( (E/27.2/137),2) +2*E/27.2); // Momentum
107 KT[eBin][aBin] = p *std::sqrt(2.-2.*std::cos(A*CLHEP::twopi/360.)); // Mom. Transfer
108 //G4cout << "aEpKt " << aBin << " " << A << " E " << E << " p " << p << " KT "
109 // << KT[eBin][aBin] << " DXS " << DXS[eBin][aBin] << G4endl;
110 }
111 }
112 }
113
114 fclose(fp);
115 }
116
117
118
119 // CDXS from DXS
120 void G4LEPTSDiffXS::BuildCDXS(G4double E, G4double El) {
121
122 for(G4int aBin=0;aBin<NumAng;aBin++) {
123 for(G4int eBin=0;eBin<=NumEn;eBin++){
124 CDXS[eBin][aBin]=0.0;
125 }
126 }
127
128 for(G4int aBin=0;aBin<NumAng;aBin++)
129 CDXS[0][aBin] = DXS[0][aBin];
130
131 for (G4int eBin=1;eBin<=NumEn;eBin++){
132 G4double sum=0.0;
133 for (G4int aBin=0;aBin<NumAng;aBin++){
134 sum += std::pow(DXS[eBin][aBin], (1.0-El/E) );
135 CDXS[eBin][aBin]=sum;
136 }
137 }
138 }
139
140
141
142 // CDXS from DXS
143 void G4LEPTSDiffXS::BuildCDXS() {
144
145 BuildCDXS(1.0, 0.0); // El = 0
146 }
147
148
149
150 // CDXS & DXS
151 void G4LEPTSDiffXS::NormalizeCDXS() {
152
153 // Normalize: 1/area
154 for (G4int eBin=1; eBin<=NumEn; eBin++){
155 G4double area = CDXS[eBin][NumAng-1];
156 //G4cout << eBin << " area = " << area << G4endl;
157
158 for (G4int aBin=0; aBin<NumAng; aBin++) {
159 CDXS[eBin][aBin] /= area;
160 //DXS[eBin][aBin] /= area;
161 }
162 }
163 }
164
165
166
167 //ICDXS from CDXS & IKT from KT
168 void G4LEPTSDiffXS::InterpolateCDXS() { // *10 angles, linear
169
170 G4double eps = 1e-5;
171 G4int ia = 0;
172
173 for( G4int aBin=0; aBin<NumAng-1; aBin++) {
174 G4double x1 = CDXS[0][aBin] + eps;
175 G4double x2 = CDXS[0][aBin+1] + eps;
176 G4double dx = (x2-x1)/100;
177
178 //if( x1<10 || x1) dx = (x2-x1)/100;
179
180 for( G4double x=x1; x < (x2-dx/10); x += dx) {
181 for( G4int eBin=0; eBin<=NumEn; eBin++) {
182 G4double y1 = CDXS[eBin][aBin];
183 G4double y2 = CDXS[eBin][aBin+1];
184 G4double z1 = KT[eBin][aBin];
185 G4double z2 = KT[eBin][aBin+1];
186
187 if( aBin==0 ) {
188 y1 /=100;
189 z1 /=100;
190 }
191
192 if( eBin==0 ) { //linear abscisa
193 ICDXS[eBin][ia] = (y1*(x2-x) + y2*(x-x1))/(x2-x1);
194 }
195 else { //log-log ordenada
196 ICDXS[eBin][ia] = G4Exp( (std::log(y1)*std::log(x2/x)+std::log(y2)*std::log(x/x1))/std::log(x2/x1) );
197 }
198
199 IKT[eBin][ia] = (z1*(x2-x) + z2*(x-x1))/(x2-x1);
200 //IKT[eBin][ia] = exp( (log(z1)*log(x2/x)+log(z2)*log(x/x1))/log(x2/x1) );
201 }
202
203 ia++;
204 }
205
206 }
207
208 INumAng = ia;
209 }
210
211
212
213 // from ICDXS
214 G4double G4LEPTSDiffXS::SampleAngle(G4double Energy) {
215 G4int ii,jj,kk=0, Ebin;
216
217 Ebin=1;
218 for(ii=2; ii<=NumEn; ii++)
219 if(Energy >= Eb[ii])
220 Ebin=ii;
221 if(Energy > Eb[NumEn]) Ebin=NumEn;
222 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
223
224 //G4cout << "SampleAngle E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
225
226 ii=0;
227 jj=INumAng-1;
228 G4double rnd=G4UniformRand();
229
230 while ((jj-ii)>1){
231 kk=(ii+jj)/2;
232 G4double dxs = ICDXS[Ebin][kk];
233 if (dxs < rnd) ii=kk;
234 else jj=kk;
235 }
236
237
238 //G4double x = ICDXS[0][jj];
239 G4double x = ICDXS[0][kk] *CLHEP::twopi/360.;
240
241 return(x);
242 }
243
244
245
246 G4double G4LEPTSDiffXS::SampleAngleEthylene(G4double E, G4double El) {
247
248 BuildCDXS(E, El);
249 NormalizeCDXS();
250 InterpolateCDXS();
251
252 return( SampleAngle(E) );
253 }
254
255
256
257 //Momentum Transfer formula
258 G4double G4LEPTSDiffXS::SampleAngleMT(G4double Energy, G4double Elost) {
259 G4int ii, jj, kk=0, Ebin, iMin, iMax;
260
261 G4double Ei = Energy;
262 G4double Ed = Energy - Elost;
263 G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
264 G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
265 G4double Kmin = Pi - Pd;
266 G4double Kmax = Pi + Pd;
267
268 if(Pd <= 1e-9 ) return (0.0);
269
270
271 // locate Energy bin
272 Ebin=1;
273 for(ii=2; ii<=NumEn; ii++)
274 if(Energy > Eb[ii]) Ebin=ii;
275 if(Energy > Eb[NumEn]) Ebin=NumEn;
276 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
277
278 //G4cout << "SampleAngle2 E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
279
280 ii=0; jj=INumAng-1;
281 while ((jj-ii)>1) {
282 kk=(ii+jj)/2;
283 if( IKT[Ebin][kk] < Kmin ) ii=kk;
284 else jj=kk;
285 }
286 iMin = ii;
287
288 ii=0; jj=INumAng-1;
289 while ((jj-ii)>1) {
290 kk=(ii+jj)/2;
291 if( IKT[Ebin][kk] < Kmax ) ii=kk;
292 else jj=kk;
293 }
294 iMax = ii;
295
296
297 // r -> a + (b-a)*r = a*(1-r) + b*r
298 G4double rnd = G4UniformRand();
299 rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[Ebin][iMax];
300 //G4double rnd = (ICDXS[Ebin][iMax] - ICDXS[Ebin][iMin]) * G4UniformRand()
301 // + ICDXS[Ebin][iMin];
302
303 ii=0; jj=INumAng-1;
304 while ((jj-ii)>1){
305 kk=(ii+jj)/2;
306 if( ICDXS[Ebin][kk] < rnd) ii=kk;
307 else jj=kk;
308 }
309
310 //Sampled
311 G4double KR = IKT[Ebin][kk];
312
313 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
314 if(co > 1) co =1;
315 G4double x = std::acos(co); //*360/twopi; //ang. dispers.
316
317 // Elastic aprox.
318 //x = 2*asin(KR/Pi/2)*360/twopi;
319
320 return(x);
321 }
322
323
324
325 void G4LEPTSDiffXS::PrintDXS(G4int NE) {
326
327 G4double dxs;
328
329 G4cout << G4endl<< "DXS & CDXS: " << fileName << G4endl<< G4endl;
330
331 for (G4int aBin=0; aBin<NumAng; aBin++) {
332 if( aBin>0)
333 dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
334 else
335 dxs = CDXS[NE][aBin];
336
337 G4cout << CDXS[0][aBin] << " " << dxs << " " << CDXS[NE][aBin] << G4endl;
338 }
339
340 G4cout << G4endl<< "IDXS & ICDXS: " << fileName << G4endl<< G4endl;
341
342 for (G4int aBin=0; aBin<INumAng; aBin++) {
343 if( aBin>0)
344 dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
345 else
346 dxs = ICDXS[NE][aBin];
347
348 G4cout << ICDXS[0][aBin] << " " << dxs << " " << ICDXS[NE][aBin] << G4endl;
349 }
350
351
352 // if(jmGlobals->VerboseHeaders) {
353 // G4cout << G4endl << "dxskt1" << G4endl;
354 // for (G4int aBin=0;aBin<NumAng;aBin++){
355 // G4cout << DXS[0][aBin] << "\t" << DXS[1][aBin] << "\t" << DXS[2][aBin] << "\t"
356 // << CDXS[1][aBin] << "\t" << KT[12][aBin] << G4endl;
357 // }
358 // }
359
360 }
361