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