| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // AMR Simplification /4 26 // AMR Simplification /4
27 // read Diff XSection & Interpolate 27 // read Diff XSection & Interpolate
>> 28 #include <string.h>
>> 29 #include <stdio.h>
>> 30 #include <string>
28 31
>> 32 #include <cmath>
29 #include "globals.hh" 33 #include "globals.hh"
30 #include "Randomize.hh" << 34 #include <iostream>
31 << 35 using namespace std;
32 #include "CLHEP/Units/PhysicalConstants.h" 36 #include "CLHEP/Units/PhysicalConstants.h"
33 37
34 #include "G4LEPTSDiffXS.hh" 38 #include "G4LEPTSDiffXS.hh"
35 #include "G4Exp.hh" <<
36 39
37 G4LEPTSDiffXS::G4LEPTSDiffXS(const G4String& f << 40
>> 41 G4LEPTSDiffXS::G4LEPTSDiffXS(string file) {
38 fileName = file; 42 fileName = file;
39 43
40 readDXS(); 44 readDXS();
41 BuildCDXS(); 45 BuildCDXS();
42 //BuildCDXS(1.0, 0.5); 46 //BuildCDXS(1.0, 0.5);
43 NormalizeCDXS(); 47 NormalizeCDXS();
44 InterpolateCDXS(); 48 InterpolateCDXS();
45 } 49 }
46 50
47 51
48 52
49 //DXS y KT 53 //DXS y KT
50 void G4LEPTSDiffXS::readDXS( ) { 54 void G4LEPTSDiffXS::readDXS( ) {
51 55
52 FILE *fp; 56 FILE *fp;
53 G4float data, data2; << 57 float data, data2;
54 58
55 if ((fp=fopen(fileName.c_str(), "r"))==nullp << 59 if ((fp=fopen(fileName.c_str(), "r"))==NULL){
56 //G4cout << "Error reading " << fileName < 60 //G4cout << "Error reading " << fileName << G4endl;
57 NumEn = 0; 61 NumEn = 0;
58 bFileFound = false; 62 bFileFound = false;
59 return; 63 return;
60 } 64 }
61 65
62 bFileFound = true; 66 bFileFound = true;
63 67
64 //G4cout << "Reading2 " << fileName << G4end 68 //G4cout << "Reading2 " << fileName << G4endl;
65 69
66 //NumAng = 181; 70 //NumAng = 181;
67 (void) fscanf(fp, "%d %d %s", &NumAng, &NumE << 71 fscanf(fp, "%d %d %s", &NumAng, &NumEn, DXSTypeName);
68 if( strcmp(DXSTypeName, "KTC") == 0 ) DX << 72 if( !strcmp(DXSTypeName, "KTC") ) DXSType = 2; // read DXS & calculate KT
69 else if( strcmp(DXSTypeName, "KT") == 0 ) DX << 73 else if( !strcmp(DXSTypeName, "KT") ) DXSType = 1; // read DXS & KT
70 else DXSTyp 74 else DXSType = 0;
71 75
72 // if( verboseLevel >= 1 ) G4cout << "Read 76 // if( verboseLevel >= 1 ) G4cout << "Read DXS (" << fileName << ")\t NEg " << NumEn << " NAng " << NumAng
73 // << "DXSType " << DXSTypeName << " " 77 // << "DXSType " << DXSTypeName << " " << DXSType << G4endl;
74 78
75 for (G4int eBin=1; eBin<=NumEn; eBin++){ 79 for (G4int eBin=1; eBin<=NumEn; eBin++){
76 (void) fscanf(fp,"%f ",&data); << 80 fscanf(fp,"%f ",&data);
77 Eb[eBin] = (G4double)data; 81 Eb[eBin] = (G4double)data;
78 } 82 }
79 83
80 84
81 //for (aBin=1;aBin<NumAng;aBin++){ 85 //for (aBin=1;aBin<NumAng;aBin++){
82 86
83 if(DXSType==1) { 87 if(DXSType==1) {
84 G4cout << "DXSTYpe 1" << G4endl; 88 G4cout << "DXSTYpe 1" << G4endl;
85 for (G4int aBin=0;aBin<NumAng;aBin++){ 89 for (G4int aBin=0;aBin<NumAng;aBin++){
86 (void) fscanf(fp,"%f ",&data); << 90 fscanf(fp,"%f ",&data);
87 DXS[0][aBin]=(G4double)data; 91 DXS[0][aBin]=(G4double)data;
88 for (G4int eBin=1;eBin<=NumEn;eBin++){ 92 for (G4int eBin=1;eBin<=NumEn;eBin++){
89 (void) fscanf(fp,"%f %f ",&data2, &data); << 93 fscanf(fp,"%f %f ",&data2, &data);
90 DXS[eBin][aBin]=(G4double)data; 94 DXS[eBin][aBin]=(G4double)data;
91 KT[eBin][aBin]=(G4double)data2; 95 KT[eBin][aBin]=(G4double)data2;
92 } 96 }
93 } 97 }
94 } 98 }
95 else { 99 else {
96 for(G4int aBin=0; aBin<NumAng; aBin++){ 100 for(G4int aBin=0; aBin<NumAng; aBin++){
97 for(G4int eBin=0; eBin<=NumEn; eBin++){ 101 for(G4int eBin=0; eBin<=NumEn; eBin++){
98 (void) fscanf(fp,"%f ",&data); << 102 fscanf(fp,"%f ",&data);
99 DXS[eBin][aBin] = (G4double)data; 103 DXS[eBin][aBin] = (G4double)data;
100 } 104 }
101 } 105 }
102 for(G4int aBin=0; aBin<NumAng; aBin++){ 106 for(G4int aBin=0; aBin<NumAng; aBin++){
103 for(G4int eBin=1; eBin<=NumEn; eBin++){ 107 for(G4int eBin=1; eBin<=NumEn; eBin++){
104 G4double A = DXS[0][aBin]; 108 G4double A = DXS[0][aBin]; // Angle
105 G4double E = Eb[eBin]; 109 G4double E = Eb[eBin]; // Energy
106 G4double p = std::sqrt(std::pow( (E/27.2/137 << 110 G4double p = sqrt(pow( (E/27.2/137),2) +2*E/27.2); // Momentum
107 KT[eBin][aBin] = p *std::sqrt(2.-2.*std::cos << 111 KT[eBin][aBin] = p *sqrt(2.-2.*cos(A*CLHEP::twopi/360.)); // Mom. Transfer
108 //G4cout << "aEpKt " << aBin << " " << A << 112 //G4cout << "aEpKt " << aBin << " " << A << " E " << E << " p " << p << " KT "
109 // << KT[eBin][aBin] << " DXS " << DXS[eBi 113 // << KT[eBin][aBin] << " DXS " << DXS[eBin][aBin] << G4endl;
110 } 114 }
111 } 115 }
112 } 116 }
113 117
114 fclose(fp); 118 fclose(fp);
115 } 119 }
116 120
117 121
118 122
119 // CDXS from DXS 123 // CDXS from DXS
120 void G4LEPTSDiffXS::BuildCDXS(G4double E, G4do 124 void G4LEPTSDiffXS::BuildCDXS(G4double E, G4double El) {
121 125
122 for(G4int aBin=0;aBin<NumAng;aBin++) { 126 for(G4int aBin=0;aBin<NumAng;aBin++) {
123 for(G4int eBin=0;eBin<=NumEn;eBin++){ 127 for(G4int eBin=0;eBin<=NumEn;eBin++){
124 CDXS[eBin][aBin]=0.0; 128 CDXS[eBin][aBin]=0.0;
125 } 129 }
126 } 130 }
127 131
128 for(G4int aBin=0;aBin<NumAng;aBin++) 132 for(G4int aBin=0;aBin<NumAng;aBin++)
129 CDXS[0][aBin] = DXS[0][aBin]; 133 CDXS[0][aBin] = DXS[0][aBin];
130 134
131 for (G4int eBin=1;eBin<=NumEn;eBin++){ 135 for (G4int eBin=1;eBin<=NumEn;eBin++){
132 G4double sum=0.0; 136 G4double sum=0.0;
133 for (G4int aBin=0;aBin<NumAng;aBin++){ 137 for (G4int aBin=0;aBin<NumAng;aBin++){
134 sum += std::pow(DXS[eBin][aBin], (1.0-El << 138 sum += pow(DXS[eBin][aBin], (1.0-El/E) );
135 CDXS[eBin][aBin]=sum; 139 CDXS[eBin][aBin]=sum;
136 } 140 }
137 } 141 }
138 } 142 }
139 143
140 144
141 145
142 // CDXS from DXS 146 // CDXS from DXS
143 void G4LEPTSDiffXS::BuildCDXS() { 147 void G4LEPTSDiffXS::BuildCDXS() {
144 148
145 BuildCDXS(1.0, 0.0); // El = 0 149 BuildCDXS(1.0, 0.0); // El = 0
146 } 150 }
147 151
148 152
149 153
150 // CDXS & DXS 154 // CDXS & DXS
151 void G4LEPTSDiffXS::NormalizeCDXS() { 155 void G4LEPTSDiffXS::NormalizeCDXS() {
152 156
153 // Normalize: 1/area 157 // Normalize: 1/area
154 for (G4int eBin=1; eBin<=NumEn; eBin++){ 158 for (G4int eBin=1; eBin<=NumEn; eBin++){
155 G4double area = CDXS[eBin][NumAng-1]; 159 G4double area = CDXS[eBin][NumAng-1];
156 //G4cout << eBin << " area = " << area << 160 //G4cout << eBin << " area = " << area << G4endl;
157 161
158 for (G4int aBin=0; aBin<NumAng; aBin++) { 162 for (G4int aBin=0; aBin<NumAng; aBin++) {
159 CDXS[eBin][aBin] /= area; 163 CDXS[eBin][aBin] /= area;
160 //DXS[eBin][aBin] /= area; 164 //DXS[eBin][aBin] /= area;
161 } 165 }
162 } 166 }
163 } 167 }
164 168
165 169
166 170
167 //ICDXS from CDXS & IKT from KT 171 //ICDXS from CDXS & IKT from KT
168 void G4LEPTSDiffXS::InterpolateCDXS() { // *1 172 void G4LEPTSDiffXS::InterpolateCDXS() { // *10 angles, linear
169 173
170 G4double eps = 1e-5; 174 G4double eps = 1e-5;
171 G4int ia = 0; 175 G4int ia = 0;
172 176
173 for( G4int aBin=0; aBin<NumAng-1; aBin++) { 177 for( G4int aBin=0; aBin<NumAng-1; aBin++) {
174 G4double x1 = CDXS[0][aBin] + eps; 178 G4double x1 = CDXS[0][aBin] + eps;
175 G4double x2 = CDXS[0][aBin+1] + eps; 179 G4double x2 = CDXS[0][aBin+1] + eps;
176 G4double dx = (x2-x1)/100; 180 G4double dx = (x2-x1)/100;
177 181
178 //if( x1<10 || x1) dx = (x2-x1)/100; 182 //if( x1<10 || x1) dx = (x2-x1)/100;
179 183
180 for( G4double x=x1; x < (x2-dx/10); x += d 184 for( G4double x=x1; x < (x2-dx/10); x += dx) {
181 for( G4int eBin=0; eBin<=NumEn; eBin++) 185 for( G4int eBin=0; eBin<=NumEn; eBin++) {
182 G4double y1 = CDXS[eBin][aBin]; 186 G4double y1 = CDXS[eBin][aBin];
183 G4double y2 = CDXS[eBin][aBin+1]; 187 G4double y2 = CDXS[eBin][aBin+1];
184 G4double z1 = KT[eBin][aBin]; 188 G4double z1 = KT[eBin][aBin];
185 G4double z2 = KT[eBin][aBin+1]; 189 G4double z2 = KT[eBin][aBin+1];
186 190
187 if( aBin==0 ) { 191 if( aBin==0 ) {
188 y1 /=100; 192 y1 /=100;
189 z1 /=100; 193 z1 /=100;
190 } 194 }
191 195
192 if( eBin==0 ) { //linear abscisa 196 if( eBin==0 ) { //linear abscisa
193 ICDXS[eBin][ia] = (y1*(x2-x) + y2*(x-x1))/ 197 ICDXS[eBin][ia] = (y1*(x2-x) + y2*(x-x1))/(x2-x1);
194 } 198 }
195 else { //log-log ordenada 199 else { //log-log ordenada
196 ICDXS[eBin][ia] = G4Exp( (std::log(y1)*std << 200 ICDXS[eBin][ia] = exp( (log(y1)*log(x2/x)+log(y2)*log(x/x1))/log(x2/x1) );
197 } 201 }
198 202
199 IKT[eBin][ia] = (z1*(x2-x) + z2*(x-x1))/(x2- 203 IKT[eBin][ia] = (z1*(x2-x) + z2*(x-x1))/(x2-x1);
200 //IKT[eBin][ia] = exp( (log(z1)*log(x2/x)+lo 204 //IKT[eBin][ia] = exp( (log(z1)*log(x2/x)+log(z2)*log(x/x1))/log(x2/x1) );
201 } 205 }
202 206
203 ia++; 207 ia++;
204 } 208 }
205 209
206 } 210 }
207 211
208 INumAng = ia; 212 INumAng = ia;
209 } 213 }
210 214
211 215
212 216
213 // from ICDXS 217 // from ICDXS
>> 218 #include "Randomize.hh"
214 G4double G4LEPTSDiffXS::SampleAngle(G4double E 219 G4double G4LEPTSDiffXS::SampleAngle(G4double Energy) {
215 G4int ii,jj,kk=0, Ebin; 220 G4int ii,jj,kk=0, Ebin;
216 221
217 Ebin=1; 222 Ebin=1;
218 for(ii=2; ii<=NumEn; ii++) 223 for(ii=2; ii<=NumEn; ii++)
219 if(Energy >= Eb[ii]) 224 if(Energy >= Eb[ii])
220 Ebin=ii; 225 Ebin=ii;
221 if(Energy > Eb[NumEn]) Ebin=NumEn; 226 if(Energy > Eb[NumEn]) Ebin=NumEn;
222 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) 227 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
223 228
224 //G4cout << "SampleAngle E " << Energy << " 229 //G4cout << "SampleAngle E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
225 230
226 ii=0; 231 ii=0;
227 jj=INumAng-1; 232 jj=INumAng-1;
228 G4double rnd=G4UniformRand(); 233 G4double rnd=G4UniformRand();
229 234
230 while ((jj-ii)>1){ 235 while ((jj-ii)>1){
231 kk=(ii+jj)/2; 236 kk=(ii+jj)/2;
232 G4double dxs = ICDXS[Ebin][kk]; 237 G4double dxs = ICDXS[Ebin][kk];
233 if (dxs < rnd) ii=kk; 238 if (dxs < rnd) ii=kk;
234 else jj=kk; 239 else jj=kk;
235 } 240 }
236 241
237 242
238 //G4double x = ICDXS[0][jj]; 243 //G4double x = ICDXS[0][jj];
239 G4double x = ICDXS[0][kk] *CLHEP::twopi/360. 244 G4double x = ICDXS[0][kk] *CLHEP::twopi/360.;
240 245
241 return(x); 246 return(x);
242 } 247 }
243 248
244 249
245 250
246 G4double G4LEPTSDiffXS::SampleAngleEthylene(G4 251 G4double G4LEPTSDiffXS::SampleAngleEthylene(G4double E, G4double El) {
247 252
248 BuildCDXS(E, El); 253 BuildCDXS(E, El);
249 NormalizeCDXS(); 254 NormalizeCDXS();
250 InterpolateCDXS(); 255 InterpolateCDXS();
251 256
252 return( SampleAngle(E) ); 257 return( SampleAngle(E) );
253 } 258 }
254 259
255 260
256 261
257 //Momentum Transfer formula 262 //Momentum Transfer formula
258 G4double G4LEPTSDiffXS::SampleAngleMT(G4double 263 G4double G4LEPTSDiffXS::SampleAngleMT(G4double Energy, G4double Elost) {
259 G4int ii, jj, kk=0, Ebin, iMin, iMax; 264 G4int ii, jj, kk=0, Ebin, iMin, iMax;
260 265
261 G4double Ei = Energy; 266 G4double Ei = Energy;
262 G4double Ed = Energy - Elost; 267 G4double Ed = Energy - Elost;
263 G4double Pi = std::sqrt( std::pow( (Ei/27.2/ << 268 G4double Pi = sqrt( pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
264 G4double Pd = std::sqrt( std::pow( (Ed/27.2/ << 269 G4double Pd = sqrt( pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
265 G4double Kmin = Pi - Pd; 270 G4double Kmin = Pi - Pd;
266 G4double Kmax = Pi + Pd; 271 G4double Kmax = Pi + Pd;
267 272
268 if(Pd <= 1e-9 ) return (0.0); 273 if(Pd <= 1e-9 ) return (0.0);
269 274
270 275
271 // locate Energy bin 276 // locate Energy bin
272 Ebin=1; 277 Ebin=1;
273 for(ii=2; ii<=NumEn; ii++) 278 for(ii=2; ii<=NumEn; ii++)
274 if(Energy > Eb[ii]) Ebin=ii; 279 if(Energy > Eb[ii]) Ebin=ii;
275 if(Energy > Eb[NumEn]) Ebin=NumEn; 280 if(Energy > Eb[NumEn]) Ebin=NumEn;
276 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) 281 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
277 282
278 //G4cout << "SampleAngle2 E " << Energy << " 283 //G4cout << "SampleAngle2 E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
279 284
280 ii=0; jj=INumAng-1; 285 ii=0; jj=INumAng-1;
281 while ((jj-ii)>1) { 286 while ((jj-ii)>1) {
282 kk=(ii+jj)/2; 287 kk=(ii+jj)/2;
283 if( IKT[Ebin][kk] < Kmin ) ii=kk; 288 if( IKT[Ebin][kk] < Kmin ) ii=kk;
284 else jj=kk; 289 else jj=kk;
285 } 290 }
286 iMin = ii; 291 iMin = ii;
287 292
288 ii=0; jj=INumAng-1; 293 ii=0; jj=INumAng-1;
289 while ((jj-ii)>1) { 294 while ((jj-ii)>1) {
290 kk=(ii+jj)/2; 295 kk=(ii+jj)/2;
291 if( IKT[Ebin][kk] < Kmax ) ii=kk; 296 if( IKT[Ebin][kk] < Kmax ) ii=kk;
292 else jj=kk; 297 else jj=kk;
293 } 298 }
294 iMax = ii; 299 iMax = ii;
295 300
296 301
297 // r -> a + (b-a)*r = a*(1-r) + b*r 302 // r -> a + (b-a)*r = a*(1-r) + b*r
298 G4double rnd = G4UniformRand(); 303 G4double rnd = G4UniformRand();
299 rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[ 304 rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[Ebin][iMax];
300 //G4double rnd = (ICDXS[Ebin][iMax] - ICDXS[ 305 //G4double rnd = (ICDXS[Ebin][iMax] - ICDXS[Ebin][iMin]) * G4UniformRand()
301 // + ICDXS[Ebin][iMin]; 306 // + ICDXS[Ebin][iMin];
302 307
303 ii=0; jj=INumAng-1; 308 ii=0; jj=INumAng-1;
304 while ((jj-ii)>1){ 309 while ((jj-ii)>1){
305 kk=(ii+jj)/2; 310 kk=(ii+jj)/2;
306 if( ICDXS[Ebin][kk] < rnd) ii=kk; 311 if( ICDXS[Ebin][kk] < rnd) ii=kk;
307 else jj=kk; 312 else jj=kk;
308 } 313 }
309 314
310 //Sampled 315 //Sampled
311 G4double KR = IKT[Ebin][kk]; 316 G4double KR = IKT[Ebin][kk];
312 317
313 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*P 318 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
314 if(co > 1) co =1; 319 if(co > 1) co =1;
315 G4double x = std::acos(co); //*360/twopi; << 320 G4double x = acos(co); //*360/twopi; //ang. dispers.
316 321
317 // Elastic aprox. 322 // Elastic aprox.
318 //x = 2*asin(KR/Pi/2)*360/twopi; 323 //x = 2*asin(KR/Pi/2)*360/twopi;
319 324
320 return(x); 325 return(x);
321 } 326 }
322 327
323 328
324 329
325 void G4LEPTSDiffXS::PrintDXS(G4int NE) { 330 void G4LEPTSDiffXS::PrintDXS(G4int NE) {
>> 331 // Debug
>> 332 //#include <string>
>> 333 //using namespace std;
>> 334
326 335
327 G4double dxs; 336 G4double dxs;
328 337
329 G4cout << G4endl<< "DXS & CDXS: " << fileNam 338 G4cout << G4endl<< "DXS & CDXS: " << fileName << G4endl<< G4endl;
330 339
331 for (G4int aBin=0; aBin<NumAng; aBin++) { 340 for (G4int aBin=0; aBin<NumAng; aBin++) {
332 if( aBin>0) 341 if( aBin>0)
333 dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1] 342 dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
334 else 343 else
335 dxs = CDXS[NE][aBin]; 344 dxs = CDXS[NE][aBin];
336 345
337 G4cout << CDXS[0][aBin] << " " << dxs << " 346 G4cout << CDXS[0][aBin] << " " << dxs << " " << CDXS[NE][aBin] << G4endl;
338 } 347 }
339 348
340 G4cout << G4endl<< "IDXS & ICDXS: " << fileN 349 G4cout << G4endl<< "IDXS & ICDXS: " << fileName << G4endl<< G4endl;
341 350
342 for (G4int aBin=0; aBin<INumAng; aBin++) { 351 for (G4int aBin=0; aBin<INumAng; aBin++) {
343 if( aBin>0) 352 if( aBin>0)
344 dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin- 353 dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
345 else 354 else
346 dxs = ICDXS[NE][aBin]; 355 dxs = ICDXS[NE][aBin];
347 356
348 G4cout << ICDXS[0][aBin] << " " << dxs << 357 G4cout << ICDXS[0][aBin] << " " << dxs << " " << ICDXS[NE][aBin] << G4endl;
349 } 358 }
350 359
351 360
352 // if(jmGlobals->VerboseHeaders) { 361 // if(jmGlobals->VerboseHeaders) {
353 // G4cout << G4endl << "dxskt1" << G4endl; 362 // G4cout << G4endl << "dxskt1" << G4endl;
354 // for (G4int aBin=0;aBin<NumAng;aBin++){ 363 // for (G4int aBin=0;aBin<NumAng;aBin++){
355 // G4cout << DXS[0][aBin] << "\t" << DXS 364 // G4cout << DXS[0][aBin] << "\t" << DXS[1][aBin] << "\t" << DXS[2][aBin] << "\t"
356 // << CDXS[1][aBin] << "\t" << KT[12][ 365 // << CDXS[1][aBin] << "\t" << KT[12][aBin] << G4endl;
357 // } 366 // }
358 // } 367 // }
359 368
360 } 369 }
361 370