Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4LEPTSDiffXS.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4LEPTSDiffXS.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4LEPTSDiffXS.cc (Version 10.6.p2)


  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   (void) 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     (void) 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       (void) 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   (void) 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   (void) 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