diff --git a/DataFormats/EcalDigi/src/EcalTimeDigi.cc b/DataFormats/EcalDigi/src/EcalTimeDigi.cc index 3e0dfe72a82d7..7aa6eaedd5bb9 100644 --- a/DataFormats/EcalDigi/src/EcalTimeDigi.cc +++ b/DataFormats/EcalDigi/src/EcalTimeDigi.cc @@ -10,11 +10,9 @@ EcalTimeDigi::EcalTimeDigi(const DetId& id) : id_(id), size_(0), sampleOfInterest_(-1), waveform_(WAVEFORMSAMPLES), data_(MAXSAMPLES) {} void EcalTimeDigi::setSize(unsigned int size) { + size_ = size; if (size > MAXSAMPLES) - size_ = MAXSAMPLES; - else - size_ = size; - data_.resize(size_); + data_.resize(size_); } void EcalTimeDigi::setWaveform(float* waveform) { diff --git a/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc b/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc index 62800f8ee5461..958470774d10a 100644 --- a/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc +++ b/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc @@ -221,14 +221,16 @@ void EcalTimeMapDigitizer::run(EcalTimeDigiCollection& output) { output.back().setWaveform(vSamAll(m_index[i])->waveform); unsigned int nTimeHits = 0; - float timeHits[vSamAll(m_index[i])->time_average_capacity]; - unsigned int timeBX[vSamAll(m_index[i])->time_average_capacity]; + unsigned int nTimeHitsMax = vSamAll(m_index[i])->time_average_capacity; + float timeHits[nTimeHitsMax]; + unsigned int timeBX[nTimeHitsMax]; - for (unsigned int j(0); j != vSamAll(m_index[i])->time_average_capacity; ++j) //here sampling on the OOTPU + for (unsigned int j(0); j != nTimeHitsMax; ++j) //here sampling on the OOTPU { if (vSamAll(m_index[i])->nhits[j] > 0) { timeHits[nTimeHits] = vSamAll(m_index[i])->average_time[j]; - timeBX[nTimeHits++] = m_minBunch + j; + timeBX[nTimeHits] = m_minBunch + j; + nTimeHits++; } } diff --git a/SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h b/SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h index b94512fcdc228..c6ab0667a853b 100644 --- a/SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h +++ b/SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h @@ -24,11 +24,7 @@ class FullModelReactionDynamics { public: FullModelReactionDynamics() {} - virtual ~FullModelReactionDynamics() {} - - virtual G4double FindInelasticity() { return 0.0; } - - virtual G4double FindTimeDelay() { return 0.0; } + ~FullModelReactionDynamics() = default; G4bool GenerateXandPt( // derived from GENXPT G4FastVector &vec, @@ -100,6 +96,9 @@ class FullModelReactionDynamics { const G4double theAtomicMass, const G4double *massVec); + FullModelReactionDynamics(const FullModelReactionDynamics &) = delete; + FullModelReactionDynamics &operator=(const FullModelReactionDynamics &) = delete; + private: void Rotate(const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, @@ -137,8 +136,6 @@ class FullModelReactionDynamics { G4int &vecLen); G4double normal(); - - G4int Poisson(G4double x); }; #endif diff --git a/SimG4Core/CustomPhysics/src/FullModelReactionDynamics.cc b/SimG4Core/CustomPhysics/src/FullModelReactionDynamics.cc index 3707086877e9f..ad011f0a1e888 100644 --- a/SimG4Core/CustomPhysics/src/FullModelReactionDynamics.cc +++ b/SimG4Core/CustomPhysics/src/FullModelReactionDynamics.cc @@ -54,6 +54,7 @@ #include "Randomize.hh" #include #include "G4ParticleTable.hh" +#include "G4Poisson.hh" using namespace CLHEP; @@ -86,8 +87,6 @@ G4bool FullModelReactionDynamics::GenerateXandPt( return false; G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus(); - // G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton(); - // G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron(); G4ParticleDefinition *aProton = G4Proton::Proton(); G4ParticleDefinition *aNeutron = G4Neutron::Neutron(); G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus(); @@ -98,7 +97,6 @@ G4bool FullModelReactionDynamics::GenerateXandPt( G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong(); G4int i, l; - // G4double forVeryForward = 0.; G4bool veryForward = false; const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV; @@ -162,8 +160,6 @@ G4bool FullModelReactionDynamics::GenerateXandPt( 0.312 + 0.200 * std::log(std::log(centerofmassEnergy * centerofmassEnergy)) + std::pow(centerofmassEnergy * centerofmassEnergy, 1.5) / 6000.0); - // PROBLEMET ER HER!!! - G4double freeEnergy = centerofmassEnergy - currentMass - targetMass; if (freeEnergy < 0) { @@ -215,11 +211,10 @@ G4bool FullModelReactionDynamics::GenerateXandPt( xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount); if (xtarg <= 0.0) xtarg = 0.01; - G4int nuclearExcitationCount = Poisson(xtarg); + G4int nuclearExcitationCount = G4Poisson(xtarg); if (atomicWeight < 1.0001) nuclearExcitationCount = 0; G4int extraNucleonCount = 0; - //G4double extraNucleonMass = 0.0; if (nuclearExcitationCount > 0) { const G4double nucsup[] = {1.00, 0.7, 0.5, 0.4, 0.35, 0.3}; const G4double psup[] = {3., 6., 20., 50., 100., 1000.}; @@ -241,7 +236,6 @@ G4bool FullModelReactionDynamics::GenerateXandPt( pVec->SetSide(-2); // -2 means backside nucleon ++extraNucleonCount; backwardEnergy += pVec->GetMass() / GeV; - //extraNucleonMass += pVec->GetMass() / GeV; } else { G4double ran = G4UniformRand(); if (ran < 0.3181) @@ -606,7 +600,6 @@ G4bool FullModelReactionDynamics::GenerateXandPt( } else { if (vec[i]->GetSide() == -2) { --extraNucleonCount; - //extraNucleonMass -= vecMass; backwardEnergy -= vecMass; } --backwardCount; @@ -855,8 +848,8 @@ G4bool FullModelReactionDynamics::GenerateXandPt( const G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20}; // Replaced the following min function to get correct behaviour on DEC. G4int tempCount = std::max(1, std::min(5, backwardNucleonCount)) - 1; - //cout << "backwardNucleonCount " << backwardNucleonCount << G4endl; - //cout << "tempCount " << tempCount << G4endl; + //G4cout << "backwardNucleonCount " << backwardNucleonCount << G4endl; + //G4cout << "tempCount " << tempCount << G4endl; G4double rmb0 = 0.0; if (targetParticle.GetSide() == -3) rmb0 += targetParticle.GetMass() / GeV; @@ -889,14 +882,14 @@ G4bool FullModelReactionDynamics::GenerateXandPt( tempV.SetElement(tempLen++, vec[i]); } if (tempLen != backwardNucleonCount) { - G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl; - G4cerr << "tempLen = " << tempLen; - G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl; - G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl; - G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl; + G4ExceptionDescription ed; + ed << "tempLen is not the same as backwardNucleonCount" << G4endl; + ed << "tempLen = " << tempLen << ", backwardNucleonCount = " << backwardNucleonCount << G4endl; + ed << "targetParticle side = " << targetParticle.GetSide() << G4endl; + ed << "currentParticle side = " << currentParticle.GetSide() << G4endl; for (i = 0; i < vecLen; ++i) - G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl; - exit(EXIT_FAILURE); + ed << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl; + G4Exception("FullModelReactionDynamics::GenerateXandPt", "had064", FatalException, ed); } constantCrossSection = true; // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); @@ -924,71 +917,13 @@ G4bool FullModelReactionDynamics::GenerateXandPt( return false; // all the secondaries have been eliminated // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); - G4int numberofFinalStateNucleons = 0; - if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron || - currentParticle.GetDefinition() == G4SigmaMinus::SigmaMinus() || - currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || - currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() || - currentParticle.GetDefinition() == G4XiZero::XiZero() || - currentParticle.GetDefinition() == G4XiMinus::XiMinus() || - currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() || - currentParticle.GetDefinition() == G4Lambda::Lambda()) - ++numberofFinalStateNucleons; + G4int numberofFinalStateNucleons = + currentParticle.GetDefinition()->GetBaryonNumber() + targetParticle.GetDefinition()->GetBaryonNumber(); currentParticle.Lorentz(currentParticle, pseudoParticle[1]); - - if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron || - targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() || - targetParticle.GetDefinition() == G4XiMinus::XiMinus() || - targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() || - targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() || - targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || - targetParticle.GetDefinition() == G4SigmaMinus::SigmaMinus()) - ++numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiProton::AntiProton()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda()) - --numberofFinalStateNucleons; targetParticle.Lorentz(targetParticle, pseudoParticle[1]); for (i = 0; i < vecLen; ++i) { - if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron || - vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() || - vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() || - vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() || - vec[i]->GetDefinition() == G4SigmaMinus::SigmaMinus()) - ++numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiProton::AntiProton()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus()) - --numberofFinalStateNucleons; + numberofFinalStateNucleons += vec[i]->GetDefinition()->GetBaryonNumber(); vec[i]->Lorentz(*vec[i], pseudoParticle[1]); } // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); @@ -1191,13 +1126,13 @@ G4bool FullModelReactionDynamics::GenerateXandPt( if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0)); if (epnb >= pnCutOff) { - npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta)); + npnb = G4Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta)); if (numberofFinalStateNucleons + npnb > atomicWeight) npnb = G4int(atomicWeight + 0.00001 - numberofFinalStateNucleons); npnb = std::min(npnb, 127 - vecLen); } if (edta >= dtaCutOff) { - ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta)); + ndta = G4Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta)); ndta = std::min(ndta, 127 - vecLen); } G4double spall = numberofFinalStateNucleons; @@ -1245,7 +1180,6 @@ void FullModelReactionDynamics::SuppressChargedPions(G4FastVector atomicNumber / atomicWeight) currentParticle.SetDefinitionAndUpdateE(aNeutron); @@ -1267,10 +1198,7 @@ void FullModelReactionDynamics::SuppressChargedPions(G4FastVectorGetDefinition() == aGamma || - vec[i]->GetDefinition() == aPiPlus || vec[i]->GetDefinition() == aPiMinus) && + if (antiTest && (vec[i]->GetDefinition() == aPiPlus || vec[i]->GetDefinition() == aPiMinus) && (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) { if (G4UniformRand() > atomicNumber / atomicWeight) vec[i]->SetDefinitionAndUpdateE(aNeutron); @@ -1313,7 +1241,6 @@ G4bool FullModelReactionDynamics::TwoCluster( G4ParticleDefinition *aNeutron = G4Neutron::Neutron(); G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus(); G4ParticleDefinition *aPiZero = G4PionZero::PionZero(); - G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus(); const G4double protonMass = aProton->GetPDGMass() / MeV; const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV; const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV; @@ -1394,12 +1321,10 @@ G4bool FullModelReactionDynamics::TwoCluster( xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount); if (xtarg <= 0.0) xtarg = 0.01; - G4int nuclearExcitationCount = Poisson(xtarg); + G4int nuclearExcitationCount = G4Poisson(xtarg); if (atomicWeight < 1.0001) nuclearExcitationCount = 0; G4int extraNucleonCount = 0; - //G4double extraMass = 0.0; - //G4double extraNucleonMass = 0.0; if (nuclearExcitationCount > 0) { G4int momentumBin = std::min(4, G4int(pOriginal / 3.0)); const G4double nucsup[] = {1.0, 0.8, 0.6, 0.5, 0.4}; @@ -1429,15 +1354,12 @@ G4bool FullModelReactionDynamics::TwoCluster( pVec->SetDefinition(aPiMinus); } pVec->SetSide(-2); // backside particle - //extraMass += pVec->GetMass() / GeV; pVec->SetNewlyAdded(true); vec.SetElement(vecLen++, pVec); // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); } } // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); - //G4double forwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + cMass - forwardMass; - //G4double backwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + bMass - backwardMass; G4double eAvailable = centerofmassEnergy - (forwardMass + backwardMass); G4bool secondaryDeleted; G4double pMass; @@ -1745,67 +1667,13 @@ G4bool FullModelReactionDynamics::TwoCluster( // // Lorentz transformation in lab system // - G4int numberofFinalStateNucleons = 0; - if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron || - currentParticle.GetDefinition() == aSigmaMinus || currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || - currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() || - currentParticle.GetDefinition() == G4XiZero::XiZero() || - currentParticle.GetDefinition() == G4XiMinus::XiMinus() || - currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() || - currentParticle.GetDefinition() == G4Lambda::Lambda()) - ++numberofFinalStateNucleons; + G4int numberofFinalStateNucleons = + currentParticle.GetDefinition()->GetBaryonNumber() + targetParticle.GetDefinition()->GetBaryonNumber(); currentParticle.Lorentz(currentParticle, pseudoParticle[2]); - if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron || - targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() || - targetParticle.GetDefinition() == G4XiMinus::XiMinus() || - targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() || - targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() || - targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || targetParticle.GetDefinition() == aSigmaMinus) - ++numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiProton::AntiProton()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus()) - --numberofFinalStateNucleons; - if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda()) - --numberofFinalStateNucleons; targetParticle.Lorentz(targetParticle, pseudoParticle[2]); + for (i = 0; i < vecLen; ++i) { - if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron || - vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() || - vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() || - vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() || - vec[i]->GetDefinition() == aSigmaMinus) - ++numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiProton::AntiProton()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus()) - --numberofFinalStateNucleons; - if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus()) - --numberofFinalStateNucleons; + numberofFinalStateNucleons += vec[i]->GetDefinition()->GetBaryonNumber(); vec[i]->Lorentz(*vec[i], pseudoParticle[2]); } // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); @@ -2023,13 +1891,13 @@ G4bool FullModelReactionDynamics::TwoCluster( sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0)); if (epnb >= pnCutOff) { - npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta)); + npnb = G4Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta)); if (numberofFinalStateNucleons + npnb > atomicWeight) npnb = G4int(atomicWeight - numberofFinalStateNucleons); npnb = std::min(npnb, 127 - vecLen); } if (edta >= dtaCutOff) { - ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta)); + ndta = G4Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta)); ndta = std::min(ndta, 127 - vecLen); } G4double spall = numberofFinalStateNucleons; @@ -2116,35 +1984,6 @@ void FullModelReactionDynamics::TwoBody(G4FastVectorGetParticleName()<GetParticleName()<= pnCutOff) { - npnb = Poisson(epnb / 0.02); + npnb = G4Poisson(epnb / 0.02); /* - G4cout<<"A couple of Poisson numbers:"< atomicWeight) @@ -2527,52 +2365,6 @@ G4double FullModelReactionDynamics::normal() { return ran; } -G4int FullModelReactionDynamics::Poisson(G4double x) // generation of poisson distribution -{ - G4int iran; - G4double ran; - - if (x > 9.9) // use normal distribution with sigma^2 = - iran = static_cast(std::max(0.0, x + normal() * std::sqrt(x))); - else { - G4int mm = G4int(5.0 * x); - if (mm <= 0) // for very small x try iran=1,2,3 - { - G4double p1 = x * std::exp(-x); - G4double p2 = x * p1 / 2.0; - G4double p3 = x * p2 / 3.0; - ran = G4UniformRand(); - if (ran < p3) - iran = 3; - else if (ran < p2) // this is original Geisha, it should be ran < p2+p3 - iran = 2; - else if (ran < p1) // should be ran < p1+p2+p3 - iran = 1; - else - iran = 0; - } else { - iran = 0; - G4double r = std::exp(-x); - ran = G4UniformRand(); - if (ran > r) { - G4double rrr; - G4double rr = r; - for (G4int i = 1; i <= mm; ++i) { - iran++; - if (i > 5) // Stirling's formula for large numbers - rrr = std::exp(i * std::log(x) - (i + 0.5) * std::log((G4double)i) + i - 0.9189385); - else - rrr = std::pow(x, i) / Factorial(i); - rr += r * rrr; - if (ran <= rr) - break; - } - } - } - } - return iran; -} - G4int FullModelReactionDynamics::Factorial(G4int n) { // calculates factorial( n ) = n*(n-1)*(n-2)*...*1 G4int m = std::min(n, 10); G4int result = 1; @@ -3261,11 +3053,6 @@ void FullModelReactionDynamics::ProduceStrangeParticlePairs(G4FastVector