Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-14 10:32:28 

0001 // @(#)root/geom:$Id$
0002 // Author: Andrei Gheata   17/06/04
0003 // Added support for radionuclides: Mihaela Gheata 24/08/2006
0004 /*************************************************************************
0005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011 
0012 #ifndef ROOT_TGeoElement
0013 #define ROOT_TGeoElement
0014 
0015 #include "TNamed.h"
0016 
0017 #include "TAttLine.h"
0018 
0019 #include "TAttFill.h"
0020 
0021 #include "TAttMarker.h"
0022 
0023 #include "TObjArray.h"
0024 
0025 #include <map>
0026 
0027 class TGeoElementTable;
0028 class TGeoIsotope;
0029 
0030 // a chemical element
0031 class TGeoElement : public TNamed {
0032 protected:
0033    enum EGeoElement { kElemUsed = BIT(17), kElemDefined = BIT(18), kElementChecked = BIT(19) };
0034 
0035    Int_t fZ;              // Z of element
0036    Int_t fN;              // Number of nucleons
0037    Int_t fNisotopes;      // Number of isotopes for the element
0038    Double_t fA;           // A of element
0039    TObjArray *fIsotopes;  // List of isotopes
0040    Double_t *fAbundances; //[fNisotopes] Array of relative isotope abundances
0041    Double_t fCoulomb;     // Coulomb correction factor
0042    Double_t fRadTsai;     // Tsai formula for the radiation length
0043 
0044 private:
0045    TGeoElement(const TGeoElement &other) = delete;
0046    TGeoElement &operator=(const TGeoElement &other) = delete;
0047 
0048    // Compute the Coulomb correction factor
0049    void ComputeCoulombFactor();
0050    // Compute the Tsai formula for the radiation length
0051    void ComputeLradTsaiFactor();
0052 
0053 public:
0054    // constructors
0055    TGeoElement();
0056    TGeoElement(const char *name, const char *title, Int_t z, Double_t a);
0057    TGeoElement(const char *name, const char *title, Int_t nisotopes);
0058    TGeoElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a);
0059    // destructor
0060    ~TGeoElement() override;
0061    // methods
0062    virtual Int_t ENDFCode() const { return 0; }
0063    Int_t Z() const { return fZ; }
0064    Int_t N() const { return fN; }
0065    Double_t Neff() const;
0066    Double_t A() const { return fA; }
0067    void AddIsotope(TGeoIsotope *isotope, Double_t relativeAbundance);
0068    Int_t GetNisotopes() const { return fNisotopes; }
0069    TGeoIsotope *GetIsotope(Int_t i) const;
0070    Double_t GetRelativeAbundance(Int_t i) const;
0071    // Calculate properties for an atomic number
0072    void ComputeDerivedQuantities();
0073    // Specific activity (in Bq/gram)
0074    virtual Double_t GetSpecificActivity() const { return 0.; }
0075    Bool_t HasIsotopes() const { return (fNisotopes == 0) ? kFALSE : kTRUE; }
0076    Bool_t IsDefined() const { return TObject::TestBit(kElemDefined); }
0077    virtual Bool_t IsRadioNuclide() const { return kFALSE; }
0078    Bool_t IsUsed() const { return TObject::TestBit(kElemUsed); }
0079    void Print(Option_t *option = "") const override;
0080    void SetDefined(Bool_t flag = kTRUE) { TObject::SetBit(kElemDefined, flag); }
0081    void SetUsed(Bool_t flag = kTRUE) { TObject::SetBit(kElemUsed, flag); }
0082    static TGeoElementTable *GetElementTable();
0083    // Coulomb correction factor
0084    inline Double_t GetfCoulomb() const { return fCoulomb; }
0085    // Tsai formula for the radiation length
0086    inline Double_t GetfRadTsai() const { return fRadTsai; }
0087 
0088    ClassDefOverride(TGeoElement, 3) // base element class
0089 };
0090 
0091 /// an isotope defined by the atomic number, number of nucleons and atomic weight (g/mole)
0092 class TGeoIsotope : public TNamed {
0093 protected:
0094    Int_t fZ;    // atomic number
0095    Int_t fN;    // number of nucleons
0096    Double_t fA; // atomic mass (g/mole)
0097 
0098 public:
0099    TGeoIsotope();
0100    TGeoIsotope(const char *name, Int_t z, Int_t n, Double_t a);
0101    ~TGeoIsotope() override {}
0102 
0103    Int_t GetZ() const { return fZ; }
0104    Int_t GetN() const { return fN; }
0105    Double_t GetA() const { return fA; }
0106    static TGeoIsotope *FindIsotope(const char *name);
0107    void Print(Option_t *option = "") const override;
0108 
0109    ClassDefOverride(TGeoIsotope, 1) // Isotope class defined by Z,N,A
0110 };
0111 
0112 class TGeoDecayChannel;
0113 class TGeoBatemanSol;
0114 
0115 /// a radionuclide
0116 class TGeoElementRN : public TGeoElement {
0117 protected:
0118    Int_t fENDFcode;        // ENDF element code
0119    Int_t fIso;             // Isomer number
0120    Double_t fLevel;        // Isomeric level
0121    Double_t fDeltaM;       // Mass excess
0122    Double_t fHalfLife;     // Half life
0123    Double_t fNatAbun;      // Natural Abundance
0124                            //   char                     fJP[11];   // Spin-parity
0125    Double_t fTH_F;         // Hynalation toxicity
0126    Double_t fTG_F;         // Ingestion toxicity
0127    Double_t fTH_S;         // Hynalation toxicity
0128    Double_t fTG_S;         // Ingestion toxicity
0129    Int_t fStatus;          // Status code
0130    TGeoBatemanSol *fRatio; // Time evolution of proportion by number
0131 
0132    TObjArray *fDecays; // List of decay modes
0133 
0134    void MakeName(Int_t a, Int_t z, Int_t iso);
0135 
0136 private:
0137    TGeoElementRN(const TGeoElementRN &elem) = delete;
0138    TGeoElementRN &operator=(const TGeoElementRN &elem) = delete;
0139 
0140 public:
0141    TGeoElementRN();
0142    TGeoElementRN(Int_t A, Int_t Z, Int_t iso, Double_t level, Double_t deltaM, Double_t halfLife, const char *JP,
0143                  Double_t natAbun, Double_t th_f, Double_t tg_f, Double_t th_s, Double_t tg_s, Int_t status);
0144    ~TGeoElementRN() override;
0145 
0146    void AddDecay(Int_t decay, Int_t diso, Double_t branchingRatio, Double_t qValue);
0147    void AddDecay(TGeoDecayChannel *dc);
0148    void AddRatio(TGeoBatemanSol &ratio);
0149    void ResetRatio();
0150    static Int_t ENDF(Int_t a, Int_t z, Int_t iso) { return 10000 * z + 10 * a + iso; }
0151 
0152    // Getters
0153    Int_t ENDFCode() const override { return fENDFcode; }
0154    Double_t GetSpecificActivity() const override;
0155    Bool_t IsRadioNuclide() const override { return kTRUE; }
0156    Int_t MassNo() const { return (Int_t)fA; }
0157    Int_t AtomicNo() const { return fZ; }
0158    Int_t IsoNo() const { return fIso; }
0159    Double_t Level() const { return fLevel; }
0160    Double_t MassEx() const { return fDeltaM; }
0161    Double_t HalfLife() const { return fHalfLife; }
0162    Double_t NatAbun() const { return fNatAbun; }
0163    const char *PJ() const { return fTitle.Data(); }
0164    Double_t TH_F() const { return fTH_F; }
0165    Double_t TG_F() const { return fTG_F; }
0166    Double_t TH_S() const { return fTH_S; }
0167    Double_t TG_S() const { return fTG_S; }
0168    Double_t Status() const { return fStatus; }
0169    Bool_t Stable() const { return !fDecays; }
0170    TObjArray *Decays() const { return fDecays; }
0171    Int_t GetNdecays() const;
0172    TGeoBatemanSol *Ratio() const { return fRatio; }
0173 
0174    // Utilities
0175    Bool_t CheckDecays() const;
0176    Int_t DecayResult(TGeoDecayChannel *dc) const;
0177    void FillPopulation(TObjArray *population, Double_t precision = 0.001, Double_t factor = 1.);
0178    void Print(Option_t *option = "") const override;
0179    static TGeoElementRN *ReadElementRN(const char *record, Int_t &ndecays);
0180    void SavePrimitive(std::ostream &out, Option_t *option = "") override;
0181 
0182    ClassDefOverride(TGeoElementRN, 2) // radionuclides class
0183 };
0184 
0185 /// decay channel utility class.
0186 class TGeoDecayChannel : public TObject {
0187 private:
0188    UInt_t fDecay;            // Decay mode
0189    Int_t fDiso;              // Delta isomeric number
0190    Double_t fBranchingRatio; // Branching Ratio
0191    Double_t fQvalue;         // Qvalue in GeV
0192    TGeoElementRN *fParent;   // Parent element
0193    TGeoElementRN *fDaughter; // Daughter element
0194 public:
0195    enum ENuclearDecayMode {
0196       kBitMask32 = 0xffffffff,
0197       k2BetaMinus = BIT(0),
0198       kBetaMinus = BIT(1),
0199       kNeutronEm = BIT(2),
0200       kProtonEm = BIT(3),
0201       kAlpha = BIT(4),
0202       kECF = BIT(5),
0203       kElecCapt = BIT(6),
0204       kIsoTrans = BIT(7),
0205       kI = BIT(8),
0206       kSpontFiss = BIT(9),
0207       k2P = BIT(10),
0208       k2N = BIT(11),
0209       k2A = BIT(12),
0210       kCarbon12 = BIT(13),
0211       kCarbon14 = BIT(14)
0212    };
0213    TGeoDecayChannel() : fDecay(0), fDiso(0), fBranchingRatio(0), fQvalue(0), fParent(nullptr), fDaughter(nullptr) {}
0214    TGeoDecayChannel(Int_t decay, Int_t diso, Double_t branchingRatio, Double_t qValue)
0215       : fDecay(decay),
0216         fDiso(diso),
0217         fBranchingRatio(branchingRatio),
0218         fQvalue(qValue),
0219         fParent(nullptr),
0220         fDaughter(nullptr)
0221    {
0222    }
0223    TGeoDecayChannel(const TGeoDecayChannel &dc)
0224       : TObject(dc),
0225         fDecay(dc.fDecay),
0226         fDiso(dc.fDiso),
0227         fBranchingRatio(dc.fBranchingRatio),
0228         fQvalue(dc.fQvalue),
0229         fParent(dc.fParent),
0230         fDaughter(dc.fDaughter)
0231    {
0232    }
0233    ~TGeoDecayChannel() override {}
0234 
0235    TGeoDecayChannel &operator=(const TGeoDecayChannel &dc);
0236 
0237    // Getters
0238    Int_t GetIndex() const;
0239    const char *GetName() const override;
0240    UInt_t Decay() const { return fDecay; }
0241    Double_t BranchingRatio() const { return fBranchingRatio; }
0242    Double_t Qvalue() const { return fQvalue; }
0243    Int_t DeltaIso() const { return fDiso; }
0244    TGeoElementRN *Daughter() const { return fDaughter; }
0245    TGeoElementRN *Parent() const { return fParent; }
0246    static void DecayName(UInt_t decay, TString &name);
0247    // Setters
0248    void SetParent(TGeoElementRN *parent) { fParent = parent; }
0249    void SetDaughter(TGeoElementRN *daughter) { fDaughter = daughter; }
0250    // Services
0251    void Print(Option_t *opt = " ") const override;
0252    static TGeoDecayChannel *ReadDecay(const char *record);
0253    void SavePrimitive(std::ostream &out, Option_t *option = "") override;
0254    virtual void DecayShift(Int_t &dA, Int_t &dZ, Int_t &dI) const;
0255 
0256    ClassDefOverride(TGeoDecayChannel, 1) // Decay channel for Elements
0257 };
0258 
0259 /// Class representing the Bateman solution for a decay branch
0260 class TGeoBatemanSol : public TObject, public TAttLine, public TAttFill, public TAttMarker {
0261 private:
0262    typedef struct {
0263       Double_t cn;     // Concentration for element 'i': Ni/Ntop
0264       Double_t lambda; // Decay coef. for element 'i'
0265    } BtCoef_t;
0266    TGeoElementRN *fElem;    // Referred RN element
0267    TGeoElementRN *fElemTop; // Top RN element
0268    Int_t fCsize;            // Size of the array of coefficients
0269    Int_t fNcoeff;           // Number of coefficients
0270    Double_t fFactor;        // Constant factor that applies to all coefficients
0271    Double_t fTmin;          // Minimum value of the time interval
0272    Double_t fTmax;          // Maximum value of the time interval
0273    BtCoef_t *fCoeff;        //[fNcoeff] Array of coefficients
0274 public:
0275    TGeoBatemanSol()
0276       : TObject(),
0277         TAttLine(),
0278         TAttFill(),
0279         TAttMarker(),
0280         fElem(nullptr),
0281         fElemTop(nullptr),
0282         fCsize(0),
0283         fNcoeff(0),
0284         fFactor(1.),
0285         fTmin(0.),
0286         fTmax(0),
0287         fCoeff(nullptr)
0288    {
0289    }
0290    TGeoBatemanSol(TGeoElementRN *elem);
0291    TGeoBatemanSol(const TObjArray *chain);
0292    TGeoBatemanSol(const TGeoBatemanSol &other);
0293    ~TGeoBatemanSol() override;
0294 
0295    TGeoBatemanSol &operator=(const TGeoBatemanSol &other);
0296    TGeoBatemanSol &operator+=(const TGeoBatemanSol &other);
0297 
0298    Double_t Concentration(Double_t time) const;
0299    void Draw(Option_t *option = "") override;
0300    void GetCoeff(Int_t i, Double_t &cn, Double_t &lambda) const
0301    {
0302       cn = fCoeff[i].cn;
0303       lambda = fCoeff[i].lambda;
0304    }
0305    void GetRange(Double_t &tmin, Double_t &tmax) const
0306    {
0307       tmin = fTmin;
0308       tmax = fTmax;
0309    }
0310    TGeoElementRN *GetElement() const { return fElem; }
0311    TGeoElementRN *GetTopElement() const { return fElemTop; }
0312    Int_t GetNcoeff() const { return fNcoeff; }
0313    void Print(Option_t *option = "") const override;
0314    void SetRange(Double_t tmin = 0., Double_t tmax = 0.)
0315    {
0316       fTmin = tmin;
0317       fTmax = tmax;
0318    }
0319    void SetFactor(Double_t factor) { fFactor = factor; }
0320    void FindSolution(const TObjArray *array);
0321    void Normalize(Double_t factor);
0322 
0323    ClassDefOverride(TGeoBatemanSol, 1) // Solution for the Bateman equation
0324 };
0325 
0326 /// iterator for decay chains.
0327 class TGeoElemIter {
0328 private:
0329    const TGeoElementRN *fTop;  // Top element of the iteration
0330    const TGeoElementRN *fElem; // Current element
0331    TObjArray *fBranch;         // Current branch
0332    Int_t fLevel;               // Current level
0333    Double_t fLimitRatio;       // Minimum cumulative branching ratio
0334    Double_t fRatio;            // Current ratio
0335 
0336 protected:
0337    TGeoElemIter() : fTop(nullptr), fElem(nullptr), fBranch(nullptr), fLevel(0), fLimitRatio(0), fRatio(0) {}
0338    TGeoElementRN *Down(Int_t ibranch);
0339    TGeoElementRN *Up();
0340 
0341 public:
0342    TGeoElemIter(TGeoElementRN *top, Double_t limit = 1.e-4);
0343    TGeoElemIter(const TGeoElemIter &iter);
0344    virtual ~TGeoElemIter();
0345 
0346    TGeoElemIter &operator=(const TGeoElemIter &iter);
0347    TGeoElementRN *operator()();
0348    TGeoElementRN *Next();
0349 
0350    TObjArray *GetBranch() const { return fBranch; }
0351    const TGeoElementRN *GetTop() const { return fTop; }
0352    const TGeoElementRN *GetElement() const { return fElem; }
0353    Int_t GetLevel() const { return fLevel; }
0354    Double_t GetRatio() const { return fRatio; }
0355    virtual void Print(Option_t *option = "") const;
0356    void SetLimitRatio(Double_t limit) { fLimitRatio = limit; }
0357 
0358    ClassDef(TGeoElemIter, 0) // Iterator for radionuclide chains.
0359 };
0360 
0361 /// table of elements
0362 class TGeoElementTable : public TObject {
0363 private:
0364    // data members
0365    Int_t fNelements;     // number of elements
0366    Int_t fNelementsRN;   // number of RN elements
0367    Int_t fNisotopes;     // number of isotopes
0368    TObjArray *fList;     // list of elements
0369    TObjArray *fListRN;   // list of RN elements
0370    TObjArray *fIsotopes; // list of user-defined isotopes
0371    // Map of radionuclides
0372    typedef std::map<Int_t, TGeoElementRN *> ElementRNMap_t;
0373    typedef ElementRNMap_t::iterator ElementRNMapIt_t;
0374    ElementRNMap_t fElementsRN; //! map of RN elements with ENDF key
0375 
0376 protected:
0377    TGeoElementTable(const TGeoElementTable &);
0378    TGeoElementTable &operator=(const TGeoElementTable &);
0379 
0380 public:
0381    // constructors
0382    TGeoElementTable();
0383    TGeoElementTable(Int_t nelements);
0384    // destructor
0385    ~TGeoElementTable() override;
0386    // methods
0387 
0388    enum EGeoETStatus { kETDefaultElements = BIT(14), kETRNElements = BIT(15) };
0389    void AddElement(const char *name, const char *title, Int_t z, Double_t a);
0390    void AddElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a);
0391    void AddElement(TGeoElement *elem);
0392    void AddElementRN(TGeoElementRN *elem);
0393    void AddIsotope(TGeoIsotope *isotope);
0394    void BuildDefaultElements();
0395    void ImportElementsRN();
0396    Bool_t CheckTable() const;
0397    TGeoElement *FindElement(const char *name) const;
0398    TGeoIsotope *FindIsotope(const char *name) const;
0399    TGeoElement *GetElement(Int_t z) { return (TGeoElement *)fList->At(z); }
0400    TGeoElementRN *GetElementRN(Int_t ENDFcode) const;
0401    TGeoElementRN *GetElementRN(Int_t a, Int_t z, Int_t iso = 0) const;
0402    TObjArray *GetElementsRN() const { return fListRN; }
0403    Bool_t HasDefaultElements() const { return TObject::TestBit(kETDefaultElements); }
0404    Bool_t HasRNElements() const { return TObject::TestBit(kETRNElements); }
0405 
0406    Int_t GetNelements() const { return fNelements; }
0407    Int_t GetNelementsRN() const { return fNelementsRN; }
0408    void ExportElementsRN(const char *filename = "");
0409    void Print(Option_t *option = "") const override;
0410 
0411    ClassDefOverride(TGeoElementTable, 4) // table of elements
0412 };
0413 
0414 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!