//////////////////////////////////////////////////////////////////// /// \class RAT::BetaFunction /// /// \brief Reads data/beta_decays.dat and sets the beta kinematic information /// /// \author Joe Formaggion and Jason Detwiler - not on SNO+ /// \author Logan Sibley -- contact person /// /// REVISION HISTORY:\n /// 04 Jan 2006 : WGS - Drag this code into the 21st century: use CLHEP, /// use strings, use STL.\n /// 28 May 2010 : Logan Sibley - changes to accommodate screening option\n /// 06 Aug 2010 : Logan Sibley - new method to sample from the Fermi /// function by storing histograms in TH1D's. \n /// 23 Sep 2010 : Logan Sibley - allow gamma decay to be handled using the /// same input format as for beta decay, but skipping the /// beta particle /// 30 Nov 2012 : Logan Sibley - enable the simulation of conversion /// electrons by specifying negative energies in the gamma /// ray array /// 03 Jan 2013 : Logan Sibley - added methods for setting parent mass and /// charge so they may be pulled into the MC data structure. /// Now requires GAMMA type decays that are electron captures /// to have Z < 0 in the data file. /// /// \details This class reads in the information stored in beta_decays.dat and /// uses it to generate the kinematic information for the electron or /// positron emitted in a beta decay. The kinetic energies of the /// electrons are sampled from TH1D's containing the spectrum shape for /// each decay branch. The spectra are calculated using /// FermiFunction.cc. Further comments are contained therein. This /// class also handles gamma decays (like electron capture), where no /// electron is emitted. In this case, the class only uses the gamma /// information, skipping the beta information. For gamma decays, the /// decay key in beta_decays.dat should be set to GAMMA: with a 0 beta /// endpoint and Z should be negative for electron captures or positive /// for isomeric transitions or internal conversions. /// //////////////////////////////////////////////////////////////////// #ifndef __RAT_BetaFunction__ #define __RAT_BetaFunction__ #include #include #include #include #include #include #include #include #include #include namespace RAT { enum EReactions { iBAD = -9999, NullParticle = -2, // for stable particles DecayBeta = +1, // for beta decay DecayEC = -1, // for positron emission DecayGamma = 0, // for gamma decay and electron capture DecayAlpha = +2 // for alpha decay }; class BetaFunction { enum EBetaFunctionConsts { PartMax = 1000, BrMax = 1000 }; public: BetaFunction(); BetaFunction(const std::string Name); BetaFunction(const std::string Name, double Z, double A = 0., int reaction = DecayBeta, double tau = 0.); ~BetaFunction(); void Reset(); void Show(); void SetTarget(const std::string Name, double Z, double A = 0., int reaction = DecayBeta, double tau = 0.); void SetBranches(double Branch, int Spin, double EndPoint, float Corr[]); void RemoveBranch(int iBranch); void SetGammas(int iBranch, int pid, double energy); void SetGammas(double energy); void SetParticleID(int iBranch, int n, int pid); void SetTargetMass(double A); void SetNorm(int iBranch, float Corr[]); int GetNParticles(int iBranch); int GetParticleID(int iBranch, int n); int GetSpin(int iBranch); double GetBranch(int iBranch); double GetEndPoint(int iBranch); double GetEnergy(int iBranch, int n); double GetValue(double energy, int iBranch, float Corr[]); void GenerateEvent(); void SetEventTime(); int GetEventID(int n); double GetEventEnergy(int n); double GetEventTotE(); Double_t GetRandomEnergy(TH1D *h); double GetRandomNumber(double rmin = 0., double rmax = 1.); bool ReadInputFile(const std::string dName); bool ReadInputFile(const std::string dName, int iType, bool screen); bool ReadInputFile(const std::string dName, int iZ, int iA, int iType = DecayBeta, bool screen = false); void ErrorLog(int iFlag = -1); inline void SetInputFileName(const std::string Name) { fInputFileName = Name; } inline void SetIsVerbose(bool iSet = true) {fIsVerbose = iSet;} inline void SetIsCulmulative(bool iSel = true) {fIsCulmulative = iSel;} inline void SetTargetName(const std::string Name) { fTargetName = Name; } inline void SetTargetCharge(double Z) {fTargetCharge = Z;} inline void SetTargetDecayType(int iType) {fTargetDecayType = iType;} inline void SetNeutrinoMass(double Mass) {fNeutrinoMass = Mass;} inline void SetTargetLifeTime(double Tau) {fTargetLifeTime = Tau;} inline void SetScreening(bool Screen) {fScreen = Screen;} inline void SetParentMass(double Mass) {fParentMass = Mass;} inline void SetParentCharge(double Z) {fParentCharge = Z;} inline int GetTargetDecayType() {return fTargetDecayType;} inline int GetNuberOfBranches() {return fNumberOfBranches;} inline int GetNGenerated() {return fNGenerated;} inline const std::string GetTargetName() {return fTargetName;} inline const std::string GetInputFileName() {return fInputFileName;} inline double GetTargetLifeTime() {return fTargetLifeTime;} inline double GetTargetMass() {return fTargetMass;} inline double GetTargetCharge() {return fTargetCharge;} inline double GetNeutrinoMass() {return fNeutrinoMass;} inline bool GetScreening() {return fScreen;} inline double GetEventTime() {return fEventTime;} inline double GetParentMass() {return fParentMass;} inline double GetParentCharge() {return fParentCharge;} //For the old sampling method /* // WGS: Initialize fDecayNorm if it hasn't been set before. inline double GetNorm(int iBranch) { std::map< size_t, double >::iterator i = fDecayNorm.find(iBranch); if ( i == fDecayNorm.end() ){ fDecayNorm[iBranch] = 1.; return 1.; // We could return fDecayNorm[iBranch], but that takes longer. } // Otherwise return fDecayNorm[iBranch]... which we already point // to with i. return (*i).second; }*/ private: bool fIsVerbose; bool fIsCulmulative; std::string fTargetName; std::string fInputFileName; int fNumberOfBranches; std::map< size_t, int > fParticleSpin; std::map< size_t, double > fParticleEndPoint; std::map< size_t, double > fParticleBranch; std::map< size_t, int > fNumberOfParticles; std::map< size_t, std::map< size_t, int > > fParticleID; std::map< size_t, std::map< size_t, double > > fParticleEnergy; double fTargetMass; double fTargetCharge; double fNeutrinoMass; double fTargetLifeTime; double fParentMass; double fParentCharge; bool fScreen; int fTargetDecayType; int fNGenerated; std::map< size_t, int > fIDn; std::map< size_t, double > fEn; double fEventTime; static TDirectory *fDirectory; // So we know where the beta histograms are }; } // namespace RAT #endif // __RAT_BetaFunction__