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