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" 39 #include "G4Exp.hh" 36 40 37 G4LEPTSDiffXS::G4LEPTSDiffXS(const G4String& f << 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] = G4Exp( (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