58 if(aDataFile >> repFlag)
61 aDataFile >> targetMass;
65 aDataFile >> nDiscrete;
66 disType =
new G4int[nDiscrete];
68 actualMult =
new G4int[nDiscrete];
70 for (
G4int i=0; i<nDiscrete; i++)
72 aDataFile >> disType[i]>>energy[i];
74 theYield[i].
Init(aDataFile,
eV);
79 aDataFile >> theInternalConversionFlag;
80 aDataFile >> theBaseEnergy;
82 aDataFile >> theInternalConversionFlag;
84 aDataFile >> nGammaEnergies;
85 theLevelEnergies =
new G4double[nGammaEnergies];
86 theTransitionProbabilities =
new G4double[nGammaEnergies];
87 if(theInternalConversionFlag == 2) thePhotonTransitionFraction =
new G4double[nGammaEnergies];
88 for(
G4int ii=0; ii<nGammaEnergies; ii++)
90 if(theInternalConversionFlag == 1)
92 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
93 theLevelEnergies[ii]*=
eV;
95 else if(theInternalConversionFlag == 2)
97 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
98 theLevelEnergies[ii]*=
eV;
102 throw G4HadronicException(__FILE__, __LINE__,
"G4NeutronHPPhotonDist: Unknown conversion flag");
110 G4cout <<
"Data representation in G4NeutronHPPhotonDist: "<<repFlag<<
G4endl;
111 throw G4HadronicException(__FILE__, __LINE__,
"G4NeutronHPPhotonDist: This data representation is not implemented.");
126 aDataFile >> isoFlag;
129 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;
130 aDataFile >> tabulationType >> nDiscrete2 >> nIso;
132 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;
135 std::vector < G4double > vct_gammas_par;
136 std::vector < G4double > vct_shells_par;
137 std::vector < G4int > vct_primary_par;
138 std::vector < G4int > vct_distype_par;
139 std::vector < G4NeutronHPVector* > vct_pXS_par;
140 if ( theGammas != NULL )
143 for ( i = 0 ; i < nDiscrete ; i++ )
145 vct_gammas_par.push_back( theGammas[ i ] );
146 vct_shells_par.push_back( theShells[ i ] );
147 vct_primary_par.push_back( isPrimary[ i ] );
148 vct_distype_par.push_back( disType[ i ] );
150 *hpv = thePartialXsec[ i ];
151 vct_pXS_par.push_back( hpv );
154 if ( theGammas == NULL ) theGammas =
new G4double[nDiscrete2];
155 if ( theShells == NULL ) theShells =
new G4double[nDiscrete2];
158 for (i=0; i< nIso; i++)
160 aDataFile >> theGammas[i] >> theShells[i];
164 nNeu =
new G4int [nDiscrete2-nIso];
167 for(i=nIso; i< nDiscrete2; i++)
169 if(tabulationType==1)
171 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
175 theLegendreManager.
Init(aDataFile);
176 for (ii=0; ii<nNeu[i-nIso]; ii++)
178 theLegendre[i-nIso][ii].
Init(aDataFile);
181 else if(tabulationType==2)
183 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
187 for (ii=0; ii<nNeu[i-nIso]; ii++)
189 theAngular[i-nIso][ii].
Init(aDataFile);
195 throw G4HadronicException(__FILE__, __LINE__,
"cannot deal with this tabulation type for angular distributions.");
199 if ( vct_gammas_par.size() > 0 )
202 for ( i = 0 ; i < nDiscrete ; i++ )
204 for (
G4int j = 0 ; j < nDiscrete ; j++ )
207 if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
209 isPrimary [ i ] = vct_primary_par [ j ];
210 disType [ i ] = vct_distype_par [ j ];
211 thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
216 for ( std::vector < G4NeutronHPVector* >::iterator
217 it = vct_pXS_par.begin() ; it != vct_pXS_par.end() ; it++ )
229 G4int i, energyDistributionsNeeded = 0;
230 for (i=0; i<nDiscrete; i++)
232 if( disType[i]==1) energyDistributionsNeeded =1;
234 if(!energyDistributionsNeeded)
return;
235 aDataFile >> nPartials;
236 distribution =
new G4int[nPartials];
241 for (i=0; i<nPartials; i++)
244 probs[i].
Init(aDataFile,
eV);
248 partials[i]->
Init(aDataFile);
256 aDataFile >> nDiscrete >> targetMass;
259 theTotalXsec.
Init(aDataFile,
eV);
262 theGammas =
new G4double[nDiscrete];
263 theShells =
new G4double[nDiscrete];
264 isPrimary =
new G4int[nDiscrete];
265 disType =
new G4int[nDiscrete];
267 for(i=0; i<nDiscrete; i++)
269 aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
272 thePartialXsec[i].
Init(aDataFile,
eV);
287 G4int nSecondaries = 0;
292 for(i=0; i<nDiscrete; i++)
294 current = theYield[i].
GetY(anEnergy);
296 if(nDiscrete==1&¤t<1.0001)
298 actualMult[i] =
static_cast<G4int>(current);
305 nSecondaries += actualMult[i];
308 for(i=0;i<nSecondaries;i++)
312 thePhotons->push_back(theOne);
324 if ( nDiscrete == 1 && nPartials == 1 )
326 if ( actualMult[ 0 ] > 0 )
328 if ( disType[0] == 1 )
365 temp = partials[ 0 ]->
GetY(anEnergy);
370 std::vector< G4double > photons_e_best( actualMult[ 0 ] , 0.0 );
373 for (
G4int j = 0 ; j < maxTry ; j++ )
375 std::vector< G4double > photons_e( actualMult[ 0 ] , 0.0 );
376 for ( std::vector< G4double >::iterator
377 it = photons_e.begin() ; it < photons_e.end() ; it++ )
381 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) > maximumE )
383 if ( std::accumulate( photons_e.begin() , photons_e.end() , 0.0 ) < best )
384 photons_e_best = photons_e;
389 for ( std::vector< G4double >::iterator
390 it = photons_e.begin() ; it < photons_e.end() ; it++ )
392 thePhotons->operator[](count)->SetKineticEnergy( *it );
400 G4cout <<
"NeutronHPPhotonDist could not find fitted energy set for multiplicity of " << actualMult[0] <<
"." <<
G4endl;
401 G4cout <<
"NeutronHPPhotonDist will use the best set." <<
G4endl;
402 for ( std::vector< G4double >::iterator
403 it = photons_e_best.begin() ; it < photons_e_best.end() ; it++ )
405 thePhotons->operator[](count)->SetKineticEnergy( *it );
416 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
419 if(count > nSecondaries)
throw G4HadronicException(__FILE__, __LINE__,
"G4NeutronHPPhotonDist::GetPhotons inconsistancy");
425 for(i=0; i<nDiscrete; i++)
427 for(ii=0; ii< actualMult[i]; ii++)
432 for(iii=0; iii<nPartials; iii++) sum+=probs[iii].GetY(anEnergy);
435 for(iii=0; iii<nPartials; iii++)
437 run+=probs[iii].
GetY(anEnergy);
439 if(random<run/sum)
break;
441 if(theP==nPartials) theP=nPartials-1;
444 temp = partials[theP]->
GetY(anEnergy);
446 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
451 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
454 if(count > nSecondaries)
throw G4HadronicException(__FILE__, __LINE__,
"G4NeutronHPPhotonDist::GetPhotons inconsistancy");
461 for (i=0; i< nSecondaries; i++)
464 G4double theta = std::acos(costheta);
467 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
468 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
469 thePhotons->operator[](i)->SetMomentum( temp ) ;
475 for(i=0; i<nSecondaries; i++)
477 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
478 for(ii=0; ii<nDiscrete2; ii++)
480 if (std::abs(currentEnergy-theGammas[ii])<0.1*
keV)
break;
482 if(ii==nDiscrete2) ii--;
489 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
490 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
491 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
493 else if(tabulationType==1)
497 for (iii=0; iii<nNeu[ii-nIso]; iii++)
500 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
504 aStore.
SetCoeff(1, &(theLegendre[ii-nIso][it]));
509 aStore.
SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
513 aStore.
SetCoeff(0, &(theLegendre[ii-nIso][it]));
519 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
520 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
521 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
527 for (iii=0; iii<nNeu[ii-nIso]; iii++)
530 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
537 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
538 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
539 thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
544 else if(repFlag == 2)
547 running[0]=theTransitionProbabilities[0];
549 for(i=1; i<nGammaEnergies; i++)
551 running[i]=running[i-1]+theTransitionProbabilities[i];
555 for(i=0; i<nGammaEnergies; i++)
558 if(random < running[i]/running[nGammaEnergies-1])
break;
561 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
565 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
576 G4double theta = std::acos(costheta);
583 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
589 for(ii=0; ii<nDiscrete2; ii++)
591 if (std::abs(currentEnergy-theGammas[ii])<0.1*
keV)
break;
593 if(ii==nDiscrete2) ii--;
607 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
610 else if(tabulationType==1)
614 for (iii=0; iii<nNeu[ii-nIso]; iii++)
617 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
621 aStore.
SetCoeff(1, &(theLegendre[ii-nIso][itt]));
626 aStore.
SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
630 aStore.
SetCoeff(0, &(theLegendre[ii-nIso][itt]));
640 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
647 for (iii=0; iii<nNeu[ii-nIso]; iii++)
650 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
661 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
665 thePhotons->push_back(theOne);
667 else if( repFlag==0 )
671 if ( thePartialXsec == 0 )
682 thePhotons->push_back( theOne );
688 std::vector < G4double > dif( nDiscrete , 0.0 );
689 for (
G4int j = 0 ; j < nDiscrete ; j++ )
703 for (
G4int j = 0 ; j < nDiscrete ; j++ )
729 if ( iphoton < nIso )
747 if ( tabulationType == 1 )
752 for (
G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
755 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy )
break;
759 aStore.
SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
760 aStore.
SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
765 else if ( tabulationType == 2 )
771 for (
G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; j++ )
774 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy )
break;
777 cosTheta = theAngular[iphoton-nIso][ iangle ].
GetCosTh();
785 G4double theta = std::acos( cosTheta );
786 G4double sinTheta = std::sin( theta );
788 G4double photonE = theGammas[ iphoton ];
789 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
791 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
G4ReactionProductVector * GetPhotons(G4double anEnergy)
G4double GetY(G4double x)
G4int GetVectorLength() const
G4long G4Poisson(G4double mean)
G4double SampleMax(G4double energy)
G4double GetTotalMomentum() const
void Init(G4int aScheme, G4int aRange)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void Init(std::istream &aDataFile)
void SetCoeff(G4int i, G4int l, G4double coeff)
G4double G4NeutronHPJENDLHEData::G4double result
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4double GetX(G4int i) const
G4double GetCosTh(G4int l)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void InitAngular(std::istream &aDataFile)
G4GLOB_DLL std::ostream G4cout
void Init(std::istream &aDataFile)
void SetTotalEnergy(const G4double en)
void InitEnergies(std::istream &aDataFile)
G4double GetY(G4int i, G4int j)
G4double GetTotalEnergy() const
G4bool InitMean(std::istream &aDataFile)
G4double GetPDGMass() const
G4double GetXsec(G4int i)
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
static G4Electron * Electron()
void InitInterpolation(G4int i, std::istream &aDataFile)
void Init(std::istream &aDataFile)
void InitPartials(std::istream &aDataFile)