00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include <numeric>
00046
00047 #include "G4NeutronHPPhotonDist.hh"
00048 #include "G4PhysicalConstants.hh"
00049 #include "G4SystemOfUnits.hh"
00050 #include "G4NeutronHPLegendreStore.hh"
00051 #include "G4Electron.hh"
00052 #include "G4Poisson.hh"
00053
00054 G4bool G4NeutronHPPhotonDist::InitMean(std::ifstream & aDataFile)
00055 {
00056
00057 G4bool result = true;
00058 if(aDataFile >> repFlag)
00059 {
00060
00061 aDataFile >> targetMass;
00062 if(repFlag==1)
00063 {
00064
00065 aDataFile >> nDiscrete;
00066 disType = new G4int[nDiscrete];
00067 energy = new G4double[nDiscrete];
00068 actualMult = new G4int[nDiscrete];
00069 theYield = new G4NeutronHPVector[nDiscrete];
00070 for (G4int i=0; i<nDiscrete; i++)
00071 {
00072 aDataFile >> disType[i]>>energy[i];
00073 energy[i]*=eV;
00074 theYield[i].Init(aDataFile, eV);
00075 }
00076 }
00077 else if(repFlag == 2)
00078 {
00079 aDataFile >> theInternalConversionFlag;
00080 aDataFile >> theBaseEnergy;
00081 theBaseEnergy*=eV;
00082 aDataFile >> theInternalConversionFlag;
00083
00084 aDataFile >> nGammaEnergies;
00085 theLevelEnergies = new G4double[nGammaEnergies];
00086 theTransitionProbabilities = new G4double[nGammaEnergies];
00087 if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies];
00088 for(G4int ii=0; ii<nGammaEnergies; ii++)
00089 {
00090 if(theInternalConversionFlag == 1)
00091 {
00092 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
00093 theLevelEnergies[ii]*=eV;
00094 }
00095 else if(theInternalConversionFlag == 2)
00096 {
00097 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
00098 theLevelEnergies[ii]*=eV;
00099 }
00100 else
00101 {
00102 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: Unknown conversion flag");
00103 }
00104 }
00105
00106
00107 }
00108 else
00109 {
00110 G4cout << "Data representation in G4NeutronHPPhotonDist: "<<repFlag<<G4endl;
00111 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist: This data representation is not implemented.");
00112 }
00113 }
00114 else
00115 {
00116 result = false;
00117 }
00118 return result;
00119 }
00120
00121 void G4NeutronHPPhotonDist::InitAngular(std::ifstream & aDataFile)
00122 {
00123
00124 G4int i, ii;
00125
00126 aDataFile >> isoFlag;
00127 if (isoFlag != 1)
00128 {
00129 if ( repFlag == 2 ) G4cout << "G4NeutronHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 Hyper News. Thanks." << G4endl;
00130 aDataFile >> tabulationType >> nDiscrete2 >> nIso;
00131
00132 if ( theGammas != NULL && nDiscrete2 != nDiscrete ) G4cout << "080731c G4NeutronHPPhotonDist nDiscrete2 != nDiscrete, It looks like something wrong in your NDL files. Please update the latest. If you still have this messages after the update, then please report to Geant4 Hyper News." << G4endl;
00133
00134
00135 std::vector < G4double > vct_gammas_par;
00136 std::vector < G4double > vct_shells_par;
00137 std::vector < G4int > vct_primary_par;
00138 std::vector < G4int > vct_distype_par;
00139 std::vector < G4NeutronHPVector* > vct_pXS_par;
00140 if ( theGammas != NULL )
00141 {
00142
00143 for ( i = 0 ; i < nDiscrete ; i++ )
00144 {
00145 vct_gammas_par.push_back( theGammas[ i ] );
00146 vct_shells_par.push_back( theShells[ i ] );
00147 vct_primary_par.push_back( isPrimary[ i ] );
00148 vct_distype_par.push_back( disType[ i ] );
00149 G4NeutronHPVector* hpv = new G4NeutronHPVector;
00150 *hpv = thePartialXsec[ i ];
00151 vct_pXS_par.push_back( hpv );
00152 }
00153 }
00154 if ( theGammas == NULL ) theGammas = new G4double[nDiscrete2];
00155 if ( theShells == NULL ) theShells = new G4double[nDiscrete2];
00156
00157
00158 for (i=0; i< nIso; i++)
00159 {
00160 aDataFile >> theGammas[i] >> theShells[i];
00161 theGammas[i]*=eV;
00162 theShells[i]*=eV;
00163 }
00164 nNeu = new G4int [nDiscrete2-nIso];
00165 if(tabulationType==1)theLegendre=new G4NeutronHPLegendreTable *[nDiscrete2-nIso];
00166 if(tabulationType==2)theAngular =new G4NeutronHPAngularP *[nDiscrete2-nIso];
00167 for(i=nIso; i< nDiscrete2; i++)
00168 {
00169 if(tabulationType==1)
00170 {
00171 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
00172 theGammas[i]*=eV;
00173 theShells[i]*=eV;
00174 theLegendre[i-nIso]=new G4NeutronHPLegendreTable[nNeu[i-nIso]];
00175 theLegendreManager.Init(aDataFile);
00176 for (ii=0; ii<nNeu[i-nIso]; ii++)
00177 {
00178 theLegendre[i-nIso][ii].Init(aDataFile);
00179 }
00180 }
00181 else if(tabulationType==2)
00182 {
00183 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
00184 theGammas[i]*=eV;
00185 theShells[i]*=eV;
00186 theAngular[i-nIso]=new G4NeutronHPAngularP[nNeu[i-nIso]];
00187 for (ii=0; ii<nNeu[i-nIso]; ii++)
00188 {
00189 theAngular[i-nIso][ii].Init(aDataFile);
00190 }
00191 }
00192 else
00193 {
00194 G4cout << "tabulation type: tabulationType"<<G4endl;
00195 throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
00196 }
00197 }
00198
00199 if ( vct_gammas_par.size() > 0 )
00200 {
00201
00202 for ( i = 0 ; i < nDiscrete ; i++ )
00203 {
00204 for ( G4int j = 0 ; j < nDiscrete ; j++ )
00205 {
00206
00207 if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
00208 {
00209 isPrimary [ i ] = vct_primary_par [ j ];
00210 disType [ i ] = vct_distype_par [ j ];
00211 thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
00212 }
00213 }
00214 }
00215
00216 for ( std::vector < G4NeutronHPVector* >::iterator
00217 it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ )
00218 {
00219 delete *it;
00220 }
00221 }
00222
00223 }
00224 }
00225
00226
00227 void G4NeutronHPPhotonDist::InitEnergies(std::ifstream & aDataFile)
00228 {
00229 G4int i, energyDistributionsNeeded = 0;
00230 for (i=0; i<nDiscrete; i++)
00231 {
00232 if( disType[i]==1) energyDistributionsNeeded =1;
00233 }
00234 if(!energyDistributionsNeeded) return;
00235 aDataFile >> nPartials;
00236 distribution = new G4int[nPartials];
00237 probs = new G4NeutronHPVector[nPartials];
00238 partials = new G4NeutronHPPartial * [nPartials];
00239 G4int nen;
00240 G4int dummy;
00241 for (i=0; i<nPartials; i++)
00242 {
00243 aDataFile >> dummy;
00244 probs[i].Init(aDataFile, eV);
00245 aDataFile >> nen;
00246 partials[i] = new G4NeutronHPPartial(nen);
00247 partials[i]->InitInterpolation(aDataFile);
00248 partials[i]->Init(aDataFile);
00249 }
00250 }
00251
00252 void G4NeutronHPPhotonDist::InitPartials(std::ifstream & aDataFile)
00253 {
00254
00255
00256 aDataFile >> nDiscrete >> targetMass;
00257 if(nDiscrete != 1)
00258 {
00259 theTotalXsec.Init(aDataFile, eV);
00260 }
00261 G4int i;
00262 theGammas = new G4double[nDiscrete];
00263 theShells = new G4double[nDiscrete];
00264 isPrimary = new G4int[nDiscrete];
00265 disType = new G4int[nDiscrete];
00266 thePartialXsec = new G4NeutronHPVector[nDiscrete];
00267 for(i=0; i<nDiscrete; i++)
00268 {
00269 aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
00270 theGammas[i]*=eV;
00271 theShells[i]*=eV;
00272 thePartialXsec[i].Init(aDataFile, eV);
00273 }
00274
00275
00276
00277
00278
00279 }
00280
00281 G4ReactionProductVector * G4NeutronHPPhotonDist::GetPhotons(G4double anEnergy)
00282 {
00283
00284
00285
00286 G4int i, ii, iii;
00287 G4int nSecondaries = 0;
00288 G4ReactionProductVector * thePhotons = new G4ReactionProductVector;
00289 if(repFlag==1)
00290 {
00291 G4double current=0;
00292 for(i=0; i<nDiscrete; i++)
00293 {
00294 current = theYield[i].GetY(anEnergy);
00295 actualMult[i] = G4Poisson(current);
00296 if(nDiscrete==1&¤t<1.0001)
00297 {
00298 actualMult[i] = static_cast<G4int>(current);
00299 if(current<1)
00300 {
00301 actualMult[i] = 0;
00302 if(G4UniformRand()<current) actualMult[i] = 1;
00303 }
00304 }
00305 nSecondaries += actualMult[i];
00306 }
00307
00308 for(i=0;i<nSecondaries;i++)
00309 {
00310 G4ReactionProduct * theOne = new G4ReactionProduct;
00311 theOne->SetDefinition(G4Gamma::Gamma());
00312 thePhotons->push_back(theOne);
00313 }
00314 G4int count=0;
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324 if ( nDiscrete == 1 && nPartials == 1 )
00325 {
00326 if ( actualMult[ 0 ] > 0 )
00327 {
00328 if ( disType[0] == 1 )
00329 {
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 G4NeutronHPVector * temp;
00365 temp = partials[ 0 ]->GetY(anEnergy);
00366 G4double maximumE = temp->GetX( temp->GetVectorLength()-1 );
00367
00368
00369
00370 std::vector< G4double > photons_e_best( actualMult[ 0 ] , 0.0 );
00371 G4double best = DBL_MAX;
00372 G4int maxTry = 1000;
00373 for ( G4int j = 0 ; j < maxTry ; j++ )
00374 {
00375 std::vector< G4double > photons_e( actualMult[ 0 ] , 0.0 );
00376 for ( std::vector< G4double >::iterator
00377 it = photons_e.begin() ; it < photons_e.end() ; it++ )
00378 {
00379 *it = temp->Sample();
00380 }
00381 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) > maximumE )
00382 {
00383 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) < best )
00384 photons_e_best = photons_e;
00385 continue;
00386 }
00387 else
00388 {
00389 for ( std::vector< G4double >::iterator
00390 it = photons_e.begin() ; it < photons_e.end() ; it++ )
00391 {
00392 thePhotons->operator[](count)->SetKineticEnergy( *it );
00393 }
00394
00395
00396
00397
00398 break;
00399 }
00400 G4cout << "NeutronHPPhotonDist could not find fitted energy set for multiplicity of " << actualMult[0] << "." << G4endl;
00401 G4cout << "NeutronHPPhotonDist will use the best set." << G4endl;
00402 for ( std::vector< G4double >::iterator
00403 it = photons_e_best.begin() ; it < photons_e_best.end() ; it++ )
00404 {
00405 thePhotons->operator[](count)->SetKineticEnergy( *it );
00406 }
00407
00408
00409
00410 }
00411
00412 delete temp;
00413 }
00414 else
00415 {
00416 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
00417 }
00418 count++;
00419 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
00420 }
00421
00422 }
00423 else
00424 {
00425 for(i=0; i<nDiscrete; i++)
00426 {
00427 for(ii=0; ii< actualMult[i]; ii++)
00428 {
00429 if(disType[i]==1)
00430 {
00431 G4double sum=0, run=0;
00432 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
00433 G4double random = G4UniformRand();
00434 G4int theP = 0;
00435 for(iii=0; iii<nPartials; iii++)
00436 {
00437 run+=probs[iii].GetY(anEnergy);
00438 theP = iii;
00439 if(random<run/sum) break;
00440 }
00441 if(theP==nPartials) theP=nPartials-1;
00442 sum=0;
00443 G4NeutronHPVector * temp;
00444 temp = partials[theP]->GetY(anEnergy);
00445 G4double eGamm = temp->Sample();
00446 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
00447 delete temp;
00448 }
00449 else
00450 {
00451 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
00452 }
00453 count++;
00454 if(count > nSecondaries) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPPhotonDist::GetPhotons inconsistancy");
00455 }
00456 }
00457 }
00458
00459 if( isoFlag == 1)
00460 {
00461 for (i=0; i< nSecondaries; i++)
00462 {
00463 G4double costheta = 2.*G4UniformRand()-1;
00464 G4double theta = std::acos(costheta);
00465 G4double phi = twopi*G4UniformRand();
00466 G4double sinth = std::sin(theta);
00467 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00468 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00469 thePhotons->operator[](i)->SetMomentum( temp ) ;
00470
00471 }
00472 }
00473 else
00474 {
00475 for(i=0; i<nSecondaries; i++)
00476 {
00477 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
00478 for(ii=0; ii<nDiscrete2; ii++)
00479 {
00480 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
00481 }
00482 if(ii==nDiscrete2) ii--;
00483 if(ii<nIso)
00484 {
00485
00486 G4double theta = pi*G4UniformRand();
00487 G4double phi = twopi*G4UniformRand();
00488 G4double sinth = std::sin(theta);
00489 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00490 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00491 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
00492 }
00493 else if(tabulationType==1)
00494 {
00495
00496 G4int it(0);
00497 for (iii=0; iii<nNeu[ii-nIso]; iii++)
00498 {
00499 it = iii;
00500 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
00501 break;
00502 }
00503 G4NeutronHPLegendreStore aStore(2);
00504 aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));
00505
00506
00507 if ( it > 0 )
00508 {
00509 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
00510 }
00511 else
00512 {
00513 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it]));
00514 }
00515 G4double cosTh = aStore.SampleMax(anEnergy);
00516 G4double theta = std::acos(cosTh);
00517 G4double phi = twopi*G4UniformRand();
00518 G4double sinth = std::sin(theta);
00519 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00520 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00521 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
00522 }
00523 else
00524 {
00525
00526 G4int it(0);
00527 for (iii=0; iii<nNeu[ii-nIso]; iii++)
00528 {
00529 it = iii;
00530 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
00531 break;
00532 }
00533 G4double costh = theAngular[ii-nIso][it].GetCosTh();
00534 G4double theta = std::acos(costh);
00535 G4double phi = twopi*G4UniformRand();
00536 G4double sinth = std::sin(theta);
00537 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
00538 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
00539 thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
00540 }
00541 }
00542 }
00543 }
00544 else if(repFlag == 2)
00545 {
00546 G4double * running = new G4double[nGammaEnergies];
00547 running[0]=theTransitionProbabilities[0];
00548
00549 for(i=1; i<nGammaEnergies; i++)
00550 {
00551 running[i]=running[i-1]+theTransitionProbabilities[i];
00552 }
00553 G4double random = G4UniformRand();
00554 G4int it=0;
00555 for(i=0; i<nGammaEnergies; i++)
00556 {
00557 it = i;
00558 if(random < running[i]/running[nGammaEnergies-1]) break;
00559 }
00560 delete [] running;
00561 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
00562 G4ReactionProduct * theOne = new G4ReactionProduct;
00563 theOne->SetDefinition(G4Gamma::Gamma());
00564 random = G4UniformRand();
00565 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
00566 {
00567 theOne->SetDefinition(G4Electron::Electron());
00568
00569
00570 totalEnergy += G4Electron::Electron()->GetPDGMass();
00571 }
00572 theOne->SetTotalEnergy(totalEnergy);
00573 if( isoFlag == 1 )
00574 {
00575 G4double costheta = 2.*G4UniformRand()-1;
00576 G4double theta = std::acos(costheta);
00577 G4double phi = twopi*G4UniformRand();
00578 G4double sinth = std::sin(theta);
00579
00580
00581 G4double en = theOne->GetTotalMomentum();
00582
00583 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00584 theOne->SetMomentum( temp ) ;
00585 }
00586 else
00587 {
00588 G4double currentEnergy = theOne->GetTotalEnergy();
00589 for(ii=0; ii<nDiscrete2; ii++)
00590 {
00591 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
00592 }
00593 if(ii==nDiscrete2) ii--;
00594 if(ii<nIso)
00595 {
00596
00597
00598
00599 G4double theta = std::acos(2.*G4UniformRand()-1.);
00600
00601 G4double phi = twopi*G4UniformRand();
00602 G4double sinth = std::sin(theta);
00603
00604
00605 G4double en = theOne->GetTotalMomentum();
00606
00607 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00608 theOne->SetMomentum( tempVector ) ;
00609 }
00610 else if(tabulationType==1)
00611 {
00612
00613 G4int itt(0);
00614 for (iii=0; iii<nNeu[ii-nIso]; iii++)
00615 {
00616 itt = iii;
00617 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
00618 break;
00619 }
00620 G4NeutronHPLegendreStore aStore(2);
00621 aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));
00622
00623
00624 if ( itt > 0 )
00625 {
00626 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
00627 }
00628 else
00629 {
00630 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt]));
00631 }
00632 G4double cosTh = aStore.SampleMax(anEnergy);
00633 G4double theta = std::acos(cosTh);
00634 G4double phi = twopi*G4UniformRand();
00635 G4double sinth = std::sin(theta);
00636
00637
00638 G4double en = theOne->GetTotalMomentum();
00639
00640 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00641 theOne->SetMomentum( tempVector ) ;
00642 }
00643 else
00644 {
00645
00646 G4int itt(0);
00647 for (iii=0; iii<nNeu[ii-nIso]; iii++)
00648 {
00649 itt = iii;
00650 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
00651 break;
00652 }
00653 G4double costh = theAngular[ii-nIso][itt].GetCosTh();
00654 G4double theta = std::acos(costh);
00655 G4double phi = twopi*G4UniformRand();
00656 G4double sinth = std::sin(theta);
00657
00658
00659 G4double en = theOne->GetTotalMomentum();
00660
00661 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
00662 theOne->SetMomentum( tmpVector ) ;
00663 }
00664 }
00665 thePhotons->push_back(theOne);
00666 }
00667 else if( repFlag==0 )
00668 {
00669
00670
00671 if ( thePartialXsec == 0 )
00672 {
00673
00674
00675 return thePhotons;
00676 }
00677
00678
00679
00680 G4ReactionProduct * theOne = new G4ReactionProduct;
00681 theOne->SetDefinition( G4Gamma::Gamma() );
00682 thePhotons->push_back( theOne );
00683
00684
00685
00686
00687 G4double sum = 0.0;
00688 std::vector < G4double > dif( nDiscrete , 0.0 );
00689 for ( G4int j = 0 ; j < nDiscrete ; j++ )
00690 {
00691 G4double x = thePartialXsec[ j ].GetXsec( anEnergy );
00692 if ( x > 0 )
00693 {
00694 sum += x;
00695 }
00696 dif [ j ] = sum;
00697
00698 }
00699
00700 G4double rand = G4UniformRand();
00701
00702 G4int iphoton = 0;
00703 for ( G4int j = 0 ; j < nDiscrete ; j++ )
00704 {
00705 G4double y = rand*sum;
00706 if ( dif [ j ] > y )
00707 {
00708 iphoton = j;
00709 break;
00710 }
00711 }
00712
00713
00714
00715
00716 G4double cosTheta = 0.0;
00717
00718 if ( isoFlag == 1 )
00719 {
00720
00721
00722
00723 cosTheta = 2.*G4UniformRand()-1;
00724
00725 }
00726 else
00727 {
00728
00729 if ( iphoton < nIso )
00730 {
00731
00732
00733
00734 cosTheta = 2.*G4UniformRand()-1;
00735
00736 }
00737 else
00738 {
00739
00740
00741
00742
00743
00744
00745
00746
00747 if ( tabulationType == 1 )
00748 {
00749
00750
00751 G4int iangle = 0;
00752 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
00753 {
00754 iangle = j;
00755 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
00756 }
00757
00758 G4NeutronHPLegendreStore aStore( 2 );
00759 aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
00760 aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
00761
00762 cosTheta = aStore.SampleMax( anEnergy );
00763
00764 }
00765 else if ( tabulationType == 2 )
00766 {
00767
00768
00769
00770 G4int iangle = 0;
00771 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
00772 {
00773 iangle = j;
00774 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
00775 }
00776
00777 cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh();
00778
00779 }
00780 }
00781 }
00782
00783
00784 G4double phi = twopi*G4UniformRand();
00785 G4double theta = std::acos( cosTheta );
00786 G4double sinTheta = std::sin( theta );
00787
00788 G4double photonE = theGammas[ iphoton ];
00789 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
00790 G4ThreeVector photonP = photonE * direction;
00791 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
00792
00793 }
00794 else
00795 {
00796 delete thePhotons;
00797 thePhotons = 0;
00798 }
00799 return thePhotons;
00800 }
00801