272{
273
274 if ( actualMult.
Get() ==
nullptr ) {
275 actualMult.
Get() =
new std::vector<G4int>( nDiscrete );
276 }
278 G4int nSecondaries = 0;
280
281 if (repFlag==1) {
283 for (i = 0; i < nDiscrete; ++i) {
284 current = theYield[i].
GetY(anEnergy);
286 if (nDiscrete == 1 && current < 1.0001) {
287 actualMult.
Get()->at(i) =
static_cast<G4int>(current);
288 if(current<1)
289 {
290 actualMult.
Get()->at(i) = 0;
292 }
293 }
294 nSecondaries += actualMult.
Get()->at(i);
295 }
296 for (i = 0; i < nSecondaries; ++i) {
299 thePhotons->push_back(theOne);
300 }
301
303
304 if (nDiscrete == 1 && nPartials == 1) {
305 if (actualMult.
Get()->at(0) > 0) {
306 if (disType[0] == 1) {
307
309 temp = partials[ 0 ]->
GetY(anEnergy);
311
312 std::vector< G4double > photons_e_best( actualMult.
Get()->at(0) , 0.0 );
315 for (
G4int j = 0; j < maxTry; ++j) {
316 std::vector<G4double> photons_e(actualMult.
Get()->at(0), 0.0);
317 for (auto it = photons_e.begin(); it < photons_e.end(); ++it) {
319 }
320
321 if (std::accumulate(photons_e.cbegin(), photons_e.cend(), 0.0) > maximumE) {
322 if (std::accumulate(photons_e.cbegin(), photons_e.cend(), 0.0) < best)
323 photons_e_best = photons_e;
324 continue;
325
326 } else {
328 for (auto it = photons_e.cbegin(); it < photons_e.cend(); ++it) {
329 thePhotons->operator[](iphot)->SetKineticEnergy(*it);
330
331
332 ++iphot;
333 }
334
335 break;
336 }
337 }
338 delete temp;
339
340 } else {
341
342 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
343 }
344 ++count;
345 if (count > nSecondaries)
346 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPPhotonDist::GetPhotons inconsistency");
347 }
348
349 } else {
350 for (i=0; i<nDiscrete; ++i) {
351 for (ii=0; ii< actualMult.
Get()->at(i); ++ii) {
352 if (disType[i] == 1) {
353
355 for (iii = 0; iii < nPartials; ++iii)
356 sum+=probs[iii].GetY(anEnergy);
359 for (iii = 0; iii < nPartials; ++iii) {
360 run+=probs[iii].
GetY(anEnergy);
361 theP = iii;
362 if(random<run/sum) break;
363 }
364
365 if (theP == nPartials) theP=nPartials-1;
366 sum = 0;
368 temp = partials[theP]->
GetY(anEnergy);
370 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
371 delete temp;
372
373 } else {
374
375 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
376 }
377 ++count;
378 if (count > nSecondaries)
379 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPPhotonDist::GetPhotons inconsistency");
380 }
381 }
382 }
383
384
385 if (isoFlag == 1) {
386 for (i=0; i< nSecondaries; ++i)
387 {
389 G4double theta = std::acos(costheta);
392 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
393 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
394 thePhotons->operator[](i)->SetMomentum( temp ) ;
395 }
396 }
397 else
398 {
399 for(i=0; i<nSecondaries; ++i)
400 {
401 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
402 for(ii=0; ii<nDiscrete2; ++ii)
403 {
404 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
405 }
406 if(ii==nDiscrete2) --ii;
407 if(ii<nIso)
408 {
409
410
411
412
414 G4double theta = std::acos(costheta);
417 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
418
419 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costheta );
420 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
421 }
422 else if(tabulationType==1)
423 {
424
426 for (iii=0; iii<nNeu[ii-nIso]; ++iii)
427 {
428 it = iii;
429 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
430 break;
431 }
433 aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));
434 if ( it > 0 )
435 {
436 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
437 }
438 else
439 {
440 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it]));
441 }
442 G4double cosTh = aStore.SampleMax(anEnergy);
446 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
447 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
448 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
449 }
450 else
451 {
452
454 for (iii=0; iii<nNeu[ii-nIso]; ++iii)
455 {
456 it = iii;
457 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
458 break;
459 }
464 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
465 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
466 thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
467 }
468 }
469 }
470
471 } else if (repFlag == 2) {
473 running[0]=theTransitionProbabilities[0];
474 for(i=1; i<nGammaEnergies; ++i)
475 {
476 running[i]=running[i-1]+theTransitionProbabilities[i];
477 }
480 for(i=0; i<nGammaEnergies; ++i)
481 {
482 it = i;
483 if(random < running[i]/running[nGammaEnergies-1]) break;
484 }
485 delete [] running;
486 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
490 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
491 {
493
494
496 }
498 if( isoFlag == 1 )
499 {
501 G4double theta = std::acos(costheta);
504
505
507
508 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
510 }
511 else
512 {
514 for(ii=0; ii<nDiscrete2; ++ii)
515 {
516 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
517 }
518 if(ii==nDiscrete2) --ii;
519 if(ii<nIso)
520 {
521
522
523
525
528
529
531
532 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
534 }
535 else if(tabulationType==1)
536 {
537
539 for (iii=0; iii<nNeu[ii-nIso]; ++iii)
540 {
541 itt = iii;
542 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
543 break;
544 }
546 aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));
547
548
549 if ( itt > 0 )
550 {
551 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
552 }
553 else
554 {
555 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt]));
556 }
557 G4double cosTh = aStore.SampleMax(anEnergy);
561
562
564
565 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
567 }
568 else
569 {
570
572 for (iii=0; iii<nNeu[ii-nIso]; ++iii)
573 {
574 itt = iii;
575 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
576 break;
577 }
582
583
585
586 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
588 }
589 }
590 thePhotons->push_back(theOne);
591 }
592 else if( repFlag==0 )
593 {
594 if ( thePartialXsec == 0 )
595 {
596 return thePhotons;
597 }
598
599
600
603 thePhotons->push_back( theOne );
604
605
606
608 std::vector < G4double > dif( nDiscrete , 0.0 );
609 for (
G4int j = 0 ; j < nDiscrete ; ++j )
610 {
612 if ( x > 0 )
613 {
614 sum += x;
615 }
616 dif [ j ] = sum;
617 }
618
620
622 for (
G4int j = 0 ; j < nDiscrete ; ++j )
623 {
625 if ( dif [ j ] > y )
626 {
627 iphoton = j;
628 break;
629 }
630 }
631
632
633
634 if (theReactionXsec) {
635 if (thePartialXsec[iphoton].GetXsec(anEnergy)/theReactionXsec->
GetXsec(anEnergy) <
G4UniformRand() ) {
636 delete thePhotons;
637 thePhotons = nullptr;
638 return thePhotons;
639 }
640 }
641
642
644
645 if ( isoFlag == 1 )
646 {
647
648
650 }
651 else
652 {
653 if ( iphoton < nIso )
654 {
655
656
658 }
659 else
660 {
661 if ( tabulationType == 1 )
662 {
663
664
666 for (
G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; ++j )
667 {
668 iangle = j;
669 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
670 }
671
673 aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
674 aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
675
676 cosTheta = aStore.SampleMax( anEnergy );
677 }
678 else if ( tabulationType == 2 )
679 {
680
681
683 for (
G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; ++j )
684 {
685 iangle = j;
686 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
687 }
688 cosTheta = theAngular[iphoton-nIso][ iangle ].
GetCosTh();
689
690 }
691 }
692 }
693
694
696 G4double theta = std::acos( cosTheta );
697 G4double sinTheta = std::sin( theta );
698
699 G4double photonE = theGammas[ iphoton ];
700 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
702 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
703 }
704 else
705 {
706 delete thePhotons;
707 thePhotons = nullptr;
708 }
709 return thePhotons;
710}
G4long G4Poisson(G4double mean)
std::vector< G4ReactionProduct * > G4ReactionProductVector
static G4Electron * Electron()
G4double GetPDGMass() const
G4double GetCosTh(G4int l)
G4double GetY(G4int i, G4int j)
G4double GetXsec(G4int i)
G4double GetY(G4double x)
G4double GetX(G4int i) const
G4int GetVectorLength() const
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
G4double GetTotalMomentum() const
G4double GetTotalEnergy() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)