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 // neutron_hp -- source file 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model. 29 // 29 // 30 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi) 31 // 080808 bug fix in Sample() and GetXsec() by 31 // 080808 bug fix in Sample() and GetXsec() by T. Koi 32 // 32 // 33 // P. Arce, June-2014 Conversion neutron_hp to 33 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 34 // 34 // 35 #include "G4ParticleHPVector.hh" 35 #include "G4ParticleHPVector.hh" 36 << 37 #include "G4SystemOfUnits.hh" 36 #include "G4SystemOfUnits.hh" 38 #include "G4Threading.hh" 37 #include "G4Threading.hh" 39 38 40 // if the ranges do not match, constant extrap << 39 // if the ranges do not match, constant extrapolation is used. 41 G4ParticleHPVector& operator+(G4ParticleHPVect << 40 G4ParticleHPVector & operator + (G4ParticleHPVector & left, G4ParticleHPVector & right) 42 { << 41 { 43 auto result = new G4ParticleHPVector; << 42 G4ParticleHPVector * result = new G4ParticleHPVector; 44 G4int j = 0; << 43 G4int j=0; 45 G4double x; << 44 G4double x; 46 G4double y; << 45 G4double y; 47 G4int running = 0; << 46 G4int running = 0; 48 for (G4int i = 0; i < left.GetVectorLength() << 47 for(G4int i=0; i<left.GetVectorLength(); i++) 49 while (j < right.GetVectorLength()) // Lo << 48 { 50 { << 49 while(j<right.GetVectorLength()) // Loop checking, 11.05.2015, T. Koi 51 if (right.GetX(j) < left.GetX(i) * 1.001 << 52 x = right.GetX(j); << 53 y = right.GetY(j) + left.GetY(x); << 54 result->SetData(running++, x, y); << 55 j++; << 56 } << 57 // else if(std::abs((right.GetX(j)-left. << 58 else if (left.GetX(i) + right.GetX(j) == << 59 || std::abs((right.GetX(j) - le << 60 { 50 { 61 x = left.GetX(i); << 51 if(right.GetX(j)<left.GetX(i)*1.001) 62 y = left.GetY(i) + right.GetY(x); << 52 { 63 result->SetData(running++, x, y); << 53 x = right.GetX(j); 64 break; << 54 y = right.GetY(j)+left.GetY(x); 65 } << 55 result->SetData(running++, x, y); 66 else { << 56 j++; 67 break; << 57 } 68 } << 58 //else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001) 69 } << 59 else if( left.GetX(i)+right.GetX(j) == 0 70 if (j == right.GetVectorLength()) { << 60 || std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j))) > 0.001 ) 71 x = left.GetX(i); << 61 { 72 y = left.GetY(i) + right.GetY(x); << 62 x = left.GetX(i); 73 result->SetData(running++, x, y); << 63 y = left.GetY(i)+right.GetY(x); 74 } << 64 result->SetData(running++, x, y); 75 } << 76 result->ThinOut(0.02); << 77 return *result; << 78 } << 79 << 80 G4ParticleHPVector::G4ParticleHPVector() << 81 { << 82 theData = new G4ParticleHPDataPoint[20]; << 83 nPoints = 20; << 84 nEntries = 0; << 85 Verbose = 0; << 86 theIntegral = nullptr; << 87 totalIntegral = -1; << 88 isFreed = 0; << 89 maxValue = -DBL_MAX; << 90 the15percentBorderCash = -DBL_MAX; << 91 the50percentBorderCash = -DBL_MAX; << 92 label = -DBL_MAX; << 93 } << 94 << 95 G4ParticleHPVector::G4ParticleHPVector(G4int n << 96 { << 97 nPoints = std::max(n, 20); << 98 theData = new G4ParticleHPDataPoint[nPoints] << 99 nEntries = 0; << 100 Verbose = 0; << 101 theIntegral = nullptr; << 102 totalIntegral = -1; << 103 isFreed = 0; << 104 maxValue = -DBL_MAX; << 105 the15percentBorderCash = -DBL_MAX; << 106 the50percentBorderCash = -DBL_MAX; << 107 label = -DBL_MAX; << 108 } << 109 << 110 G4ParticleHPVector::~G4ParticleHPVector() << 111 { << 112 // if(Verbose==1)G4cout <<"G4ParticleHPVe << 113 delete[] theData; << 114 // if(Verbose==1)G4cout <<"Vector: delete << 115 delete[] theIntegral; << 116 // if(Verbose==1)G4cout <<"Vector: delete << 117 theHash.Clear(); << 118 isFreed = 1; << 119 } << 120 << 121 G4ParticleHPVector& G4ParticleHPVector::operat << 122 { << 123 if (&right == this) return *this; << 124 << 125 G4int i; << 126 << 127 totalIntegral = right.totalIntegral; << 128 if (right.theIntegral != nullptr) theIntegra << 129 for (i = 0; i < right.nEntries; i++) { << 130 SetPoint(i, right.GetPoint(i)); // copy t << 131 if (right.theIntegral != nullptr) theInteg << 132 } << 133 theManager = right.theManager; << 134 label = right.label; << 135 << 136 Verbose = right.Verbose; << 137 the15percentBorderCash = right.the15percentB << 138 the50percentBorderCash = right.the50percentB << 139 theHash = right.theHash; << 140 return *this; << 141 } << 142 << 143 G4double G4ParticleHPVector::GetXsec(G4double << 144 { << 145 if (nEntries == 0) return 0; << 146 // if(!theHash.Prepared()) Hash(); << 147 if (!theHash.Prepared()) { << 148 if (G4Threading::IsWorkerThread()) { << 149 ; << 150 } << 151 else { << 152 Hash(); << 153 } << 154 } << 155 G4int min = theHash.GetMinIndex(e); << 156 G4int i; << 157 for (i = min; i < nEntries; i++) { << 158 // if(theData[i].GetX()>e) break; << 159 if (theData[i].GetX() >= e) break; << 160 } << 161 G4int low = i - 1; << 162 G4int high = i; << 163 if (i == 0) { << 164 low = 0; << 165 high = 1; << 166 } << 167 else if (i == nEntries) { << 168 low = nEntries - 2; << 169 high = nEntries - 1; << 170 } << 171 G4double y; << 172 if (e < theData[nEntries - 1].GetX()) { << 173 // Protect against doubled-up x values << 174 // if( (theData[high].GetX()-theData[low]. << 175 if (theData[high].GetX() != 0 << 176 // 080808 TKDB << 177 //&&( theData[high].GetX()-theData[low << 178 && (std::abs((theData[high].GetX() - t << 179 < 0.000001)) << 180 { << 181 y = theData[low].GetY(); << 182 } << 183 else { << 184 y = theInt.Interpolate(theManager.GetSch << 185 theData[high].Get << 186 } << 187 } << 188 else { << 189 y = theData[nEntries - 1].GetY(); << 190 } << 191 return y; << 192 } << 193 << 194 G4int G4ParticleHPVector::GetEnergyIndex(G4dou << 195 { << 196 // returns energy index of the bin in which << 197 // returns energy index of emax for NJOY << 198 if (nEntries == 0) return 0; << 199 if ( ( ! theHash.Prepared() ) && ( ! G4Threa << 200 G4int min = theHash.GetMinIndex(e); << 201 G4int i = 0; << 202 for (i=min ; i < nEntries; i++) if (theData[ << 203 return i; << 204 } << 205 << 206 void G4ParticleHPVector::Dump() << 207 { << 208 G4cout << nEntries << G4endl; << 209 for (G4int i = 0; i < nEntries; i++) { << 210 G4cout << theData[i].GetX() << " "; << 211 G4cout << theData[i].GetY() << " "; << 212 // if (i!=1&&i==5*(i/5)) G4cout << G4 << 213 G4cout << G4endl; << 214 } << 215 G4cout << G4endl; << 216 } << 217 << 218 void G4ParticleHPVector::Check(G4int i) << 219 { << 220 if (i > nEntries) << 221 throw G4HadronicException(__FILE__, __LINE << 222 "Skipped some in << 223 if (i == nPoints) { << 224 nPoints = static_cast<G4int>(1.2 * nPoints << 225 auto buff = new G4ParticleHPDataPoint[nPoi << 226 for (G4int j = 0; j < nEntries; j++) << 227 buff[j] = theData[j]; << 228 delete[] theData; << 229 theData = buff; << 230 } << 231 if (i == nEntries) nEntries = i + 1; << 232 } << 233 << 234 void G4ParticleHPVector::Merge(G4Interpolation << 235 G4ParticleHPVec << 236 { << 237 // interpolate between labels according to a << 238 // continue in unknown areas by substraction << 239 << 240 CleanUp(); << 241 G4int s_tmp = 0, n = 0, m_tmp = 0; << 242 G4ParticleHPVector* tmp; << 243 G4int a = s_tmp, p = n, t; << 244 while (a < active->GetVectorLength()) // Lo << 245 { << 246 if (active->GetEnergy(a) <= passive->GetEn << 247 G4double xa = active->GetEnergy(a); << 248 G4double yy = theInt.Interpolate(aScheme << 249 active- << 250 SetData(m_tmp, xa, yy); << 251 theManager.AppendScheme(m_tmp, active->G << 252 m_tmp++; << 253 a++; << 254 G4double xp = passive->GetEnergy(p); << 255 // if( std::abs(std::abs(xp-xa)/xa)<0.00 << 256 if (xa != 0 && std::abs(std::abs(xp - xa << 257 { << 258 p++; << 259 tmp = active; << 260 t = a; << 261 active = passive; << 262 a = p; << 263 passive = tmp; << 264 p = t; << 265 } << 266 } << 267 else { << 268 tmp = active; << 269 t = a; << 270 active = passive; << 271 a = p; << 272 passive = tmp; << 273 p = t; << 274 } << 275 } << 276 << 277 G4double deltaX = passive->GetXsec(GetEnergy << 278 while (p != passive->GetVectorLength() << 279 && passive->GetEnergy(p) <= aValue) << 280 { << 281 G4double anX; << 282 anX = passive->GetXsec(p) - deltaX; << 283 if (anX > 0) { << 284 // if(std::abs(GetEnergy(m-1)-passive->G << 285 if (passive->GetEnergy(p) == 0 << 286 || std::abs(GetEnergy(m_tmp - 1) - p << 287 > 0.0000001) << 288 { << 289 SetData(m_tmp, passive->GetEnergy(p), << 290 theManager.AppendScheme(m_tmp++, passi << 291 } << 292 } << 293 p++; << 294 } << 295 // Rebuild the Hash; << 296 if (theHash.Prepared()) { << 297 ReHash(); << 298 } << 299 } << 300 << 301 void G4ParticleHPVector::ThinOut(G4double prec << 302 { << 303 // anything in there? << 304 if (GetVectorLength() == 0) return; << 305 // make the new vector << 306 auto aBuff = new G4ParticleHPDataPoint[nPoin << 307 G4double x, x1, x2, y, y1, y2; << 308 G4int count = 0, current = 2, start = 1; << 309 << 310 // First element always goes and is never te << 311 aBuff[0] = theData[0]; << 312 << 313 // Find the rest << 314 while (current < GetVectorLength()) // Loop << 315 { << 316 x1 = aBuff[count].GetX(); << 317 y1 = aBuff[count].GetY(); << 318 x2 = theData[current].GetX(); << 319 y2 = theData[current].GetY(); << 320 << 321 if (x1 - x2 == 0) { << 322 // Following block added for avoiding di << 323 for (G4int j = start; j < current; j++) << 324 y = (y2 + y1) / 2.; << 325 if (std::abs(y - theData[j].GetY()) > << 326 aBuff[++count] = theData[current - 1 << 327 start = current; // the next candid << 328 break; 65 break; 329 } 66 } 330 } << 331 } << 332 else { << 333 for (G4int j = start; j < current; j++) << 334 x = theData[j].GetX(); << 335 if (x1 - x2 == 0) << 336 y = (y2 + y1) / 2.; << 337 else 67 else 338 y = theInt.Lin(x, x1, x2, y1, y2); << 68 { 339 if (std::abs(y - theData[j].GetY()) > << 340 aBuff[++count] = theData[current - 1 << 341 start = current; // the next candid << 342 break; 69 break; 343 } 70 } 344 } 71 } >> 72 if(j==right.GetVectorLength()) >> 73 { >> 74 x = left.GetX(i); >> 75 y = left.GetY(i)+right.GetY(x); >> 76 result->SetData(running++, x, y); >> 77 } 345 } 78 } 346 current++; << 79 result->ThinOut(0.02); >> 80 return *result; 347 } 81 } 348 // The last one also always goes, and is nev << 82 349 aBuff[++count] = theData[GetVectorLength() - << 83 G4ParticleHPVector::G4ParticleHPVector() 350 delete[] theData; << 84 { 351 theData = aBuff; << 85 theData = new G4ParticleHPDataPoint[20]; 352 nEntries = count + 1; << 86 nPoints=20; 353 << 87 nEntries=0; 354 // Rebuild the Hash; << 88 Verbose=0; 355 if (theHash.Prepared()) { << 89 theIntegral=0; 356 ReHash(); << 90 totalIntegral=-1; 357 } << 91 isFreed = 0; 358 } << 92 maxValue = -DBL_MAX; 359 << 93 the15percentBorderCash = -DBL_MAX; 360 G4bool G4ParticleHPVector::IsBlocked(G4double << 94 the50percentBorderCash = -DBL_MAX; 361 { << 95 label = -DBL_MAX; 362 G4bool result = false; << 96 } 363 std::vector<G4double>::iterator i; << 97 364 for (i = theBlocked.begin(); i != theBlocked << 98 G4ParticleHPVector::G4ParticleHPVector(G4int n) 365 G4double aBlock = *i; << 99 { 366 if (std::abs(aX - aBlock) < 0.1 * MeV) { << 100 nPoints=std::max(n, 20); 367 result = true; << 101 theData = new G4ParticleHPDataPoint[nPoints]; 368 theBlocked.erase(i); << 102 nEntries=0; 369 break; << 103 Verbose=0; 370 } << 104 theIntegral=0; 371 } << 105 totalIntegral=-1; 372 return result; << 106 isFreed = 0; 373 } << 107 maxValue = -DBL_MAX; 374 << 108 the15percentBorderCash = -DBL_MAX; 375 G4double G4ParticleHPVector::Sample() // Samp << 109 the50percentBorderCash = -DBL_MAX; 376 { << 110 label = -DBL_MAX; 377 G4double result = 0.; << 111 } 378 G4int j; << 112 379 for (j = 0; j < GetVectorLength(); j++) { << 113 G4ParticleHPVector::~G4ParticleHPVector() 380 if (GetY(j) < 0) SetY(j, 0); << 114 { 381 } << 115 // if(Verbose==1)G4cout <<"G4ParticleHPVector::~G4ParticleHPVector"<<G4endl; 382 << 116 delete [] theData; 383 if (!theBuffered.empty() && G4UniformRand() << 117 // if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl; 384 result = theBuffered[0]; << 118 delete [] theIntegral; 385 theBuffered.erase(theBuffered.begin()); << 119 // if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl; 386 if (result < GetX(GetVectorLength() - 1)) << 120 theHash.Clear(); 387 } << 121 isFreed = 1; 388 if (GetVectorLength() == 1) { << 122 } 389 result = theData[0].GetX(); << 123 390 } << 124 G4ParticleHPVector & G4ParticleHPVector:: 391 else { << 125 operator = (const G4ParticleHPVector & right) 392 if (theIntegral == nullptr) { << 126 { 393 IntegrateAndNormalise(); << 127 if(&right == this) return *this; 394 } << 128 395 G4int icounter = 0; << 129 G4int i; 396 G4int icounter_max = 1024; << 130 397 do { << 131 totalIntegral = right.totalIntegral; 398 icounter++; << 132 if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries]; 399 if (icounter > icounter_max) { << 133 for(i=0; i<right.nEntries; i++) 400 G4cout << "Loop-counter exceeded the t << 134 { 401 << __FILE__ << "." << G4endl; << 135 SetPoint(i, right.GetPoint(i)); // copy theData 402 break; << 136 if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i]; 403 } << 137 } 404 // 080808 << 138 theManager = right.theManager; 405 /* << 139 label = right.label; 406 G4double rand; << 140 407 G4double value, test, baseline; << 141 Verbose = right.Verbose; 408 baseline = theData[GetVectorLeng << 142 the15percentBorderCash = right.the15percentBorderCash; 409 do << 143 the50percentBorderCash = right.the50percentBorderCash; 410 { << 144 theHash = right.theHash; 411 value = baseline*G4UniformRand << 145 return *this; 412 value += theData[0].GetX(); << 146 } 413 test = GetY(value)/maxValue; << 147 414 rand = G4UniformRand(); << 148 415 } << 149 G4double G4ParticleHPVector::GetXsec(G4double e) 416 //while(test<rand); << 150 { 417 while( test < rand && test > 0 ) << 151 if(nEntries == 0) return 0; 418 result = value; << 152 //if(!theHash.Prepared()) Hash(); 419 */ << 153 if ( !theHash.Prepared() ) { 420 G4double rand; << 154 if ( G4Threading::IsWorkerThread() ) { 421 G4double value = 0., test; << 155 ; 422 G4int jcounter = 0; << 156 } else { 423 G4int jcounter_max = 1024; << 157 Hash(); 424 do { << 158 } 425 jcounter++; << 159 } 426 if (jcounter > jcounter_max) { << 160 G4int min = theHash.GetMinIndex(e); 427 G4cout << "Loop-counter exceeded the << 161 G4int i; 428 << __FILE__ << "." << G4endl; << 162 for(i=min ; i<nEntries; i++) 429 break; << 163 { 430 } << 164 //if(theData[i].GetX()>e) break; 431 rand = G4UniformRand(); << 165 if(theData[i].GetX() >= e) break; 432 G4int ibin = -1; << 166 } 433 for (G4int i = 0; i < GetVectorLength( << 167 G4int low = i-1; 434 if (rand < theIntegral[i]) { << 168 G4int high = i; 435 ibin = i; << 169 if(i==0) 436 break; << 170 { 437 } << 171 low = 0; >> 172 high = 1; >> 173 } >> 174 else if(i==nEntries) >> 175 { >> 176 low = nEntries-2; >> 177 high = nEntries-1; >> 178 } >> 179 G4double y; >> 180 if(e<theData[nEntries-1].GetX()) >> 181 { >> 182 // Protect against doubled-up x values >> 183 //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001) >> 184 if ( theData[high].GetX() !=0 >> 185 //080808 TKDB >> 186 //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001) >> 187 &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) ) >> 188 { >> 189 y = theData[low].GetY(); >> 190 } >> 191 else >> 192 { >> 193 y = theInt.Interpolate(theManager.GetScheme(high), e, >> 194 theData[low].GetX(), theData[high].GetX(), >> 195 theData[low].GetY(), theData[high].GetY()); >> 196 } >> 197 } >> 198 else >> 199 { >> 200 y=theData[nEntries-1].GetY(); >> 201 } >> 202 return y; >> 203 } >> 204 >> 205 void G4ParticleHPVector::Dump() >> 206 { >> 207 G4cout << nEntries<<G4endl; >> 208 for(G4int i=0; i<nEntries; i++) >> 209 { >> 210 G4cout << theData[i].GetX()<<" "; >> 211 G4cout << theData[i].GetY()<<" "; >> 212 // if (i!=1&&i==5*(i/5)) G4cout << G4endl; >> 213 G4cout << G4endl; >> 214 } >> 215 G4cout << G4endl; >> 216 } >> 217 >> 218 void G4ParticleHPVector::Check(G4int i) >> 219 { >> 220 if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4ParticleHPVector"); >> 221 if(i==nPoints) >> 222 { >> 223 nPoints = static_cast<G4int>(1.2*nPoints); >> 224 G4ParticleHPDataPoint * buff = new G4ParticleHPDataPoint[nPoints]; >> 225 for (G4int j=0; j<nEntries; j++) buff[j] = theData[j]; >> 226 delete [] theData; >> 227 theData = buff; >> 228 } >> 229 if(i==nEntries) nEntries=i+1; >> 230 } >> 231 >> 232 void G4ParticleHPVector:: >> 233 Merge(G4InterpolationScheme aScheme, G4double aValue, >> 234 G4ParticleHPVector * active, G4ParticleHPVector * passive) >> 235 { >> 236 // interpolate between labels according to aScheme, cut at aValue, >> 237 // continue in unknown areas by substraction of the last difference. >> 238 >> 239 CleanUp(); >> 240 G4int s_tmp = 0, n=0, m_tmp=0; >> 241 G4ParticleHPVector * tmp; >> 242 G4int a = s_tmp, p = n, t; >> 243 while ( a<active->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi >> 244 { >> 245 if(active->GetEnergy(a) <= passive->GetEnergy(p)) >> 246 { >> 247 G4double xa = active->GetEnergy(a); >> 248 G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(), >> 249 active->GetXsec(a), passive->GetXsec(xa)); >> 250 SetData(m_tmp, xa, yy); >> 251 theManager.AppendScheme(m_tmp, active->GetScheme(a)); >> 252 m_tmp++; >> 253 a++; >> 254 G4double xp = passive->GetEnergy(p); >> 255 //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() ) >> 256 if ( xa != 0 >> 257 && std::abs(std::abs(xp-xa)/xa) < 0.0000001 >> 258 && a < active->GetVectorLength() ) >> 259 { >> 260 p++; >> 261 tmp = active; t=a; >> 262 active = passive; a=p; >> 263 passive = tmp; p=t; 438 } 264 } 439 if (ibin < 0) G4cout << "TKDB 080807 " << 265 } else { 440 // result << 266 tmp = active; t=a; 441 rand = G4UniformRand(); << 267 active = passive; a=p; 442 G4double x1, x2; << 268 passive = tmp; p=t; 443 if (ibin == 0) { << 269 } 444 x1 = theData[ibin].GetX(); << 270 } 445 value = x1; << 271 446 break; << 272 G4double deltaX = passive->GetXsec(GetEnergy(m_tmp-1)) - GetXsec(m_tmp-1); >> 273 while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue) // Loop checking, 11.05.2015, T. Koi >> 274 { >> 275 G4double anX; >> 276 anX = passive->GetXsec(p)-deltaX; >> 277 if(anX>0) >> 278 { >> 279 //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001) >> 280 if ( passive->GetEnergy(p) == 0 >> 281 || std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 ) >> 282 { >> 283 SetData(m_tmp, passive->GetEnergy(p), anX); >> 284 theManager.AppendScheme(m_tmp++, passive->GetScheme(p)); 447 } 285 } 448 x1 = theData[ibin - 1].GetX(); << 286 } >> 287 p++; >> 288 } >> 289 // Rebuild the Hash; >> 290 if(theHash.Prepared()) >> 291 { >> 292 ReHash(); >> 293 } >> 294 } >> 295 >> 296 void G4ParticleHPVector::ThinOut(G4double precision) >> 297 { >> 298 // anything in there? >> 299 if(GetVectorLength()==0) return; >> 300 // make the new vector >> 301 G4ParticleHPDataPoint * aBuff = new G4ParticleHPDataPoint[nPoints]; >> 302 G4double x, x1, x2, y, y1, y2; >> 303 G4int count = 0, current = 2, start = 1; >> 304 >> 305 // First element always goes and is never tested. >> 306 aBuff[0] = theData[0]; >> 307 >> 308 // Find the rest >> 309 while(current < GetVectorLength()) // Loop checking, 11.05.2015, T. Koi >> 310 { >> 311 x1=aBuff[count].GetX(); >> 312 y1=aBuff[count].GetY(); >> 313 x2=theData[current].GetX(); >> 314 y2=theData[current].GetY(); >> 315 >> 316 if ( x1-x2 == 0 ) { >> 317 //Following block added for avoiding div 0 error on Release + G4FPE_DEBUG >> 318 for ( G4int j=start; j<current; j++ ) { >> 319 y = (y2+y1)/2.; >> 320 if ( std::abs( y-theData[j].GetY() ) > precision*y ) { >> 321 aBuff[++count] = theData[current-1]; // for this one, everything was fine >> 322 start = current; // the next candidate >> 323 break; >> 324 } >> 325 } >> 326 } else { >> 327 for(G4int j=start; j<current; j++) >> 328 { >> 329 x = theData[j].GetX(); >> 330 if(x1-x2 == 0) y = (y2+y1)/2.; >> 331 else y = theInt.Lin(x, x1, x2, y1, y2); >> 332 if (std::abs(y-theData[j].GetY())>precision*y) >> 333 { >> 334 aBuff[++count] = theData[current-1]; // for this one, everything was fine >> 335 start = current; // the next candidate >> 336 break; >> 337 } >> 338 } >> 339 } >> 340 current++ ; >> 341 } >> 342 // The last one also always goes, and is never tested. >> 343 aBuff[++count] = theData[GetVectorLength()-1]; >> 344 delete [] theData; >> 345 theData = aBuff; >> 346 nEntries = count+1; >> 347 >> 348 // Rebuild the Hash; >> 349 if(theHash.Prepared()) >> 350 { >> 351 ReHash(); >> 352 } >> 353 } 449 354 450 x2 = theData[ibin].GetX(); << 355 G4bool G4ParticleHPVector::IsBlocked(G4double aX) 451 value = rand * (x2 - x1) + x1; << 356 { 452 //************************************ << 357 G4bool result = false; 453 /* << 358 std::vector<G4double>::iterator i; 454 test = GetY ( value ) / std::max << 359 for(i=theBlocked.begin(); i!=theBlocked.end(); i++) 455 */ << 360 { 456 //************************************ << 361 G4double aBlock = *i; 457 // EMendoza - Always linear interpolat << 362 if(std::abs(aX-aBlock) < 0.1*MeV) 458 G4double y1 = theData[ibin - 1].GetY() << 363 { 459 G4double y2 = theData[ibin].GetY(); << 364 result = true; 460 G4double mval = (y2 - y1) / (x2 - x1); << 365 theBlocked.erase(i); 461 G4double bval = y1 - mval * x1; << 366 break; 462 test = (mval * value + bval) / std::ma << 367 } 463 //************************************ << 464 } while (G4UniformRand() > test); // Lo << 465 result = value; << 466 // 080807 << 467 } while (IsBlocked(result)); // Loop chec << 468 } << 469 return result; << 470 } << 471 << 472 G4double G4ParticleHPVector::Get15percentBorde << 473 { << 474 if (the15percentBorderCash > -DBL_MAX / 2.) << 475 G4double result; << 476 if (GetVectorLength() == 1) { << 477 result = theData[0].GetX(); << 478 the15percentBorderCash = result; << 479 } << 480 else { << 481 if (theIntegral == nullptr) { << 482 IntegrateAndNormalise(); << 483 } 368 } 484 G4int i; << 369 return result; 485 result = theData[GetVectorLength() - 1].Ge << 370 } 486 for (i = 0; i < GetVectorLength(); i++) { << 371 487 if (theIntegral[i] / theIntegral[GetVect << 372 G4double G4ParticleHPVector::Sample() // Samples X according to distribution Y 488 result = theData[std::min(i + 1, GetVe << 373 { 489 the15percentBorderCash = result; << 374 G4double result; 490 break; << 375 G4int j; 491 } << 376 for(j=0; j<GetVectorLength(); j++) 492 } << 377 { 493 the15percentBorderCash = result; << 378 if(GetY(j)<0) SetY(j, 0); 494 } << 495 return result; << 496 } << 497 << 498 G4double G4ParticleHPVector::Get50percentBorde << 499 { << 500 if (the50percentBorderCash > -DBL_MAX / 2.) << 501 G4double result; << 502 if (GetVectorLength() == 1) { << 503 result = theData[0].GetX(); << 504 the50percentBorderCash = result; << 505 } << 506 else { << 507 if (theIntegral == nullptr) { << 508 IntegrateAndNormalise(); << 509 } 379 } 510 G4int i; << 380 511 G4double x = 0.5; << 381 if(theBuffered.size() !=0 && G4UniformRand()<0.5) 512 result = theData[GetVectorLength() - 1].Ge << 382 { 513 for (i = 0; i < GetVectorLength(); i++) { << 383 result = theBuffered[0]; 514 if (theIntegral[i] / theIntegral[GetVect << 384 theBuffered.erase(theBuffered.begin()); 515 G4int it; << 385 if(result < GetX(GetVectorLength()-1) ) return result; 516 it = i; << 386 } 517 if (it == GetVectorLength() - 1) { << 387 if(GetVectorLength()==1) 518 result = theData[GetVectorLength() - << 388 { >> 389 result = theData[0].GetX(); >> 390 } >> 391 else >> 392 { >> 393 if(theIntegral==0) { IntegrateAndNormalise(); } >> 394 G4int icounter=0; >> 395 G4int icounter_max=1024; >> 396 do >> 397 { >> 398 icounter++; >> 399 if ( icounter > icounter_max ) { >> 400 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; >> 401 break; 519 } 402 } 520 else { << 403 //080808 521 G4double x1, x2, y1, y2; << 404 /* 522 x1 = theIntegral[i - 1] / theIntegra << 405 G4double rand; 523 x2 = theIntegral[i] / theIntegral[Ge << 406 G4double value, test, baseline; 524 y1 = theData[i - 1].GetX(); << 407 baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX(); 525 y2 = theData[i].GetX(); << 408 do 526 result = theLin.Lin(x, x1, x2, y1, y << 409 { >> 410 value = baseline*G4UniformRand(); >> 411 value += theData[0].GetX(); >> 412 test = GetY(value)/maxValue; >> 413 rand = G4UniformRand(); >> 414 } >> 415 //while(test<rand); >> 416 while( test < rand && test > 0 ); >> 417 result = value; >> 418 */ >> 419 G4double rand; >> 420 G4double value, test; >> 421 G4int jcounter=0; >> 422 G4int jcounter_max=1024; >> 423 do >> 424 { >> 425 jcounter++; >> 426 if ( jcounter > jcounter_max ) { >> 427 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; >> 428 break; >> 429 } >> 430 rand = G4UniformRand(); >> 431 G4int ibin = -1; >> 432 for ( G4int i = 0 ; i < GetVectorLength() ; i++ ) >> 433 { >> 434 if ( rand < theIntegral[i] ) >> 435 { >> 436 ibin = i; >> 437 break; >> 438 } >> 439 } >> 440 if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl; >> 441 // result >> 442 rand = G4UniformRand(); >> 443 G4double x1, x2; >> 444 if ( ibin == 0 ) >> 445 { >> 446 x1 = theData[ ibin ].GetX(); >> 447 value = x1; >> 448 break; >> 449 } >> 450 else >> 451 { >> 452 x1 = theData[ ibin-1 ].GetX(); >> 453 } >> 454 >> 455 x2 = theData[ ibin ].GetX(); >> 456 value = rand * ( x2 - x1 ) + x1; >> 457 //*********************************************************************** >> 458 /* >> 459 test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) ); >> 460 */ >> 461 //*********************************************************************** >> 462 //EMendoza - Always linear interpolation: >> 463 G4double y1=theData[ ibin-1 ].GetY(); >> 464 G4double y2=theData[ ibin ].GetY(); >> 465 G4double mval=(y2-y1)/(x2-x1); >> 466 G4double bval=y1-mval*x1; >> 467 test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) ); >> 468 //*********************************************************************** 527 } 469 } 528 the50percentBorderCash = result; << 470 while ( G4UniformRand() > test ); // Loop checking, 11.05.2015, T. Koi 529 break; << 471 result = value; >> 472 //080807 530 } 473 } >> 474 while(IsBlocked(result)); // Loop checking, 11.05.2015, T. Koi 531 } 475 } 532 the50percentBorderCash = result; << 476 return result; 533 } 477 } 534 return result; << 535 } << 536 478 537 G4double G4ParticleHPVector::GetMaxY(G4double << 479 G4double G4ParticleHPVector::Get15percentBorder() 538 { << 480 { 539 G4double xsmax = 0.; << 481 if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash; 540 if (emin > emax || nEntries == 0) return xsm << 482 G4double result; 541 if (emin >= theData[nEntries - 1].GetX()) { << 483 if(GetVectorLength()==1) 542 xsmax = theData[nEntries - 1].GetY(); << 484 { 543 return xsmax; << 485 result = theData[0].GetX(); 544 } << 486 the15percentBorderCash = result; 545 if (emax <= theData[0].GetX()) { << 487 } 546 xsmax = theData[0].GetY(); << 488 else 547 return xsmax; << 489 { 548 } << 490 if(theIntegral==0) { IntegrateAndNormalise(); } 549 if (!theHash.Prepared() && !G4Threading::IsW << 491 G4int i; 550 // Find the lowest index, low, where x(energ << 492 result = theData[GetVectorLength()-1].GetX(); 551 G4int i = theHash.GetMinIndex(emin); << 493 for(i=0;i<GetVectorLength();i++) 552 for (; i < nEntries; ++i) { << 494 { 553 if (theData[i].GetX() >= emin) break; << 495 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15) 554 } << 496 { 555 G4int low = i; << 497 result = theData[std::min(i+1, GetVectorLength()-1)].GetX(); 556 // Find the lowest index, high, where x(ener << 498 the15percentBorderCash = result; 557 i = theHash.GetMinIndex(emax); << 499 break; 558 for (; i < nEntries; ++i) { << 500 } 559 if (theData[i].GetX() >= emax) break; << 501 } 560 } << 502 the15percentBorderCash = result; 561 G4int high = i; << 503 } 562 xsmax = GetXsec(emin); // Set xsmax as low << 504 return result; 563 // Find the highest cross-section << 564 for (i = low; i < high; ++i) { << 565 if (xsmax < theData[i].GetY()) xsmax = the << 566 } << 567 // Check if it is smaller than the high bord << 568 G4double highborder = GetXsec(emax); << 569 if (xsmax < highborder) xsmax = highborder; << 570 << 571 if (xsmax == 0.) { << 572 throw G4HadronicException(__FILE__, __LINE << 573 "G4ParticleHPVec << 574 "G4Nucleus::GetB << 575 } 505 } 576 506 577 return xsmax; << 507 G4double G4ParticleHPVector::Get50percentBorder() 578 } << 508 { >> 509 if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash; >> 510 G4double result; >> 511 if(GetVectorLength()==1) >> 512 { >> 513 result = theData[0].GetX(); >> 514 the50percentBorderCash = result; >> 515 } >> 516 else >> 517 { >> 518 if(theIntegral==0) { IntegrateAndNormalise(); } >> 519 G4int i; >> 520 G4double x = 0.5; >> 521 result = theData[GetVectorLength()-1].GetX(); >> 522 for(i=0;i<GetVectorLength();i++) >> 523 { >> 524 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x) >> 525 { >> 526 G4int it; >> 527 it = i; >> 528 if(it == GetVectorLength()-1) >> 529 { >> 530 result = theData[GetVectorLength()-1].GetX(); >> 531 } >> 532 else >> 533 { >> 534 G4double x1, x2, y1, y2; >> 535 x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1]; >> 536 x2 = theIntegral[i]/theIntegral[GetVectorLength()-1]; >> 537 y1 = theData[i-1].GetX(); >> 538 y2 = theData[i].GetX(); >> 539 result = theLin.Lin(x, x1, x2, y1, y2); >> 540 } >> 541 the50percentBorderCash = result; >> 542 break; >> 543 } >> 544 } >> 545 the50percentBorderCash = result; >> 546 } >> 547 return result; >> 548 } 579 549