ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEPTSDiffXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 }