gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
decay_chain.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17 
18 #ifndef __decay_chain_hpp__
19 #define __decay_chain_hpp__
20 
21 #include <vector>
22 #include <unordered_map>
23 #include <string>
24 #include <set>
25 #include <boost/shared_ptr.hpp>
27 
28 namespace Gambit
29 {
30  namespace DarkBit
31  {
32  struct TH_Channel;
33  struct TH_ProcessCatalog;
34  class SimYieldTable;
35  namespace DecayChain
36  {
37  using std::vector;
38  using std::ofstream;
39  using std::ostream;
40  using std::string;
41  using boost::shared_ptr;
42  using std::unordered_map;
43  using std::set;
44 
45  // *********************************************
46  // Generic 3-vector class
47  // *********************************************
48  class vec3
49  {
50  public:
51  double vals[3];
52  // Constructors
53  vec3(){}
54  vec3(double v0, double v1, double v2) {vals[0]=v0;vals[1]=v1;vals[2]=v2;}
55  vec3(double v0){vals[0]=v0;vals[1]=v0;vals[2]=v0;}
56  // Get the length of the vector
57  double length() const;
58  // Normalize vector to unit length
59  void normalize();
60  // Normalize vector to length len
61  void normalize(const double len);
62  // Operators
63  double& operator[](int i){ return vals[i];}
64  double operator[](int i)const{ return vals[i];}
65  vec3 operator-()const{return vec3(-vals[0],-vals[1],-vals[2]); }
66  };
67  vec3 operator* (double x, const vec3 &y);
68  vec3 operator* (const vec3 &y, double x);
69  vec3 operator/ (const vec3 &y, double x);
70  ostream& operator<<(ostream& os, const vec3& v);
71  double dot(const vec3 &a, const vec3 &b);
72 
73 
74  // *********************************************
75  // Generic 4-vector class
76  // *********************************************
77  class vec4
78  {
79  public:
80  double vals[4];
81  // Constructors
82  vec4(){}
83  vec4(double v0, double v1, double v2, double v3){vals[0]=v0;vals[1]=v1;vals[2]=v2;vals[3]=v3;}
84  vec4(double v0){vals[0]=v0;vals[1]=v0;vals[2]=v0;vals[3]=v0;}
85  vec4(double v0, vec3 v){vals[0]=v0;vals[1]=v[0];vals[2]=v[1];vals[3]=v[2];}
86  // Returns vec3 containting elements 1,2,3
87  vec3 xyz() const{return vec3(vals[1],vals[2],vals[3]);}
88  // Operators
89  double& operator[](int i){ return vals[i];}
90  double operator[](int i)const{ return vals[i];}
91  vec4 operator-(){return vec4(-vals[0],-vals[1],-vals[2],-vals[3]);}
92  };
93  vec4 operator* (double x, const vec4 &y);
94  vec4 operator* (const vec4 &y, double x);
95  vec4 operator+ (const vec4 &x, const vec4 &y);
96  vec4 operator- (const vec4 &x, const vec4 &y);
97  ostream& operator<<(ostream& os, const vec4& v);
98  double dot(const vec4 &a, const vec4 &b);
99  // Construct an energy-momentum 4-vector from the given 3-momentum and mass
100  vec4 Ep4vec(const vec3 p, double m);
101 
102 
103  // *********************************************
104  // Generic 4x4 class
105  // *********************************************
106  class mat4
107  {
108  public:
109  double vals[4][4];
110  mat4(){};
111  mat4( double v00, double v01, double v02, double v03,
112  double v10, double v11, double v12, double v13,
113  double v20, double v21, double v22, double v23,
114  double v30, double v31, double v32, double v33);
115  mat4(double v);
116  // Identity matrix
117  static mat4 identity();
118  };
119  vec4 operator* (const mat4 &m, const vec4 &v);
120  mat4 operator* (const mat4 &m1, const mat4 &m2);
121  ostream& operator<<(ostream& os, const mat4& m);
122 
123 
124  // *********************************************
125  // Utility functions
126  // *********************************************
127  // Generate a random number between -1 and 1
128  double rand_m1_1();
129  // Generate a random number between 0 and 1
130  inline double rand_0_1(){return Random::draw();}
131  // Generate a 3-vector to a random point on the unit sphere
132  vec3 randOnSphere();
133  // Calculate Lorentz boost matrix corresponding to beta_xyz.
134  // Always calculate gamma=E/m and use the second version when possible, as this is far more numerically stable.
135  void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat);
136  void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat, double gamma);
137  // Boost inVec according to beta_xyz.
138  vec4 lorentzBoost(const vec4 &inVec, const vec3 &beta_xyz);
139  // Boost inVec to the frame where a particle at rest (in this frame) would have 4-momentum p_parent
140  vec4 p_parentFrame(const vec4 &inVec, const vec4 &p_parent);
141  // Get Lorentz matrix for boosting to parent frame
142  void boostMatrixParentFrame(mat4 &mat, vec4 &p_parent);
143  // Calculate the invariant mass of the given 4-vector pair
144  double invariantMass(const vec4 &a, const vec4 &b);
145 
146 
147  // *********************************************
148  // Class containing all the allowed decay channels of a single particle,
149  // as well as the particle mass.
150  // Contains functonality for picking random decays according to their decay widths
151  // *********************************************
153  {
154  friend class DecayTable;
155  public:
156  // Particle mass
157  const double m;
158  // Is the particle stable in the context of the decay chain?
159  bool stable;
160  // Flags indicating whether or not the various decay states are endpoints
161  unordered_map<const TH_Channel*, bool> endpointFlags;
162  // Constructor
163  DecayTableEntry(string pID, double m, bool stable) :
164  m(m), stable(stable), enabledWidth(0),
165  totalWidth(0), invisibleWidth(0), useForcedTotalWidth(false), forcedTotalWidth(0), pID(pID), randInit(false){}
166  // Dummy constructor (for the DecayTable map)
168  m(0), stable(false), enabledWidth(0),
169  totalWidth(0), invisibleWidth(0), useForcedTotalWidth(false), forcedTotalWidth(0), pID(""), randInit(false){}
170  // Pick a random decay channel
171  bool randomDecay(const TH_Channel* &decay) const;
172  // Update decay widths and if necessary Monte Carlo table
173  void update();
174  // Functions for checking if a decay exists in the decay tables
175  bool isEnabled(const TH_Channel *in) const;
176  bool isDisabled(const TH_Channel *in) const;
177  bool isRegistered(const TH_Channel *in) const;
178  // Add decay channel as enabled or disabled
179  void addChannel(const TH_Channel *in);
180  void addDisabled(const TH_Channel *in);
181  void setInvisibleWidth(double width);
182  // Functions for enabling/disabling (registered) decays.
183  // Returns false if decay is in the list of disabled/enabled decays
184  bool enableDecay(const TH_Channel *in);
185  bool disableDecay(const TH_Channel *in);
186  // Get branching ratio of enabled decays
187  double getEnabledBranching() const;
188  // Set total decay width to a fixed value
189  void forceTotalWidth(bool enabled, double width);
190  // Get total decay width
191  double getTotalWidth() const;
192  bool hasEnabledDecays() const;
193  private:
194  // Lists of decays
195  vector<const TH_Channel*> enabledDecays;
196  vector<const TH_Channel*> disabledDecays;
197  // Table used for picking random decays
198  mutable vector<double> randLims;
199  // Decay widths of enabled channels and all channels
200  double enabledWidth;
201  double totalWidth;
205  // String particle identifier
206  string pID;
207  // Initialization status of Monte Carlo table
208  mutable bool randInit;
209  // Function used in picking a random decay
210  int findChannelIdx(double pick) const;
211  // Generate table for picking a random decay in Monte Carlo decay chains
212  void generateRandTable() const;
213  };
214 
215  // *********************************************
216  // Table of all particles and their decay channels.
217  // Uses particle PID as array index
218  // *********************************************
220  {
221  public:
222  DecayTable(const TH_ProcessCatalog &cat, const SimYieldTable &tab, set<string> disabledList);
224  bool hasEntry(string) const;
225  // Add particle to decay table, specifying particle ID, mass and whether or not it should be decayed in decay chains
226  void addEntry(string pID, double m, bool stable);
227  void addEntry(string pID, DecayTableEntry entry);
228  bool randomDecay(string pID, const TH_Channel* &decay) const;
229  const DecayTableEntry& operator[](string i) const;
230  // Retrieve width of decay channel
231  static double getWidth(const TH_Channel *ch);
232  // Print the decay table (to cout)
233  void printTable() const;
234  private:
235  unordered_map<string,DecayTableEntry> table;
236  };
237 
238 
239  // *********************************************
240  // The main decay chain class.
241  // Each link (particle) in the decay chain is an instance of this class,
242  // with pointers to its parent and child links.
243  // *********************************************
245  {
246  public:
247  // Invariant (rest) mass
248  const double m;
249  // Constructor for the base node (top particle in the decay chain).
250  ChainParticle(vec3 ipLab, const DecayTable *dc, string pID);
251  // Iteratively add random links to the decay chain by Monte Carlo to a maximum length of maxSteps or minimum energy of Emin.
252  // Use negative numbers to turn off limits
253  void generateDecayChainMC(int maxSteps, double Emin);
254  // Draw new angles for the decay products in this and all subsequent links of the decay chain
255  void reDrawAngles();
256  // Remove all subsequent links in the decay chain
257  void cutChain();
258  // Boost a given 4-momentum from this frame to the lab frame.
259  vec4 p_to_Lab(const vec4 &p) const;
260  // Calculate 4-momentum of this particle in the lab frame.
261  vec4 p_Lab() const;
262  // Calculate the energy of this particle in the lab frame. Faster than calculating the full 4-momentum.
263  double E_Lab() const;
264  // Iteratively collect decay chain endpoint states. Optional argument to only collect particles of certain types.
265  // Note that this function will not collect decay products of endpoint states. These must be extracted manually.
266  void collectEndpointStates(vector<const ChainParticle*> &endpointStates, bool includeAborted, string ipID="") const;
267  // Get number of child particles
268  int getnChildren() const {return nChildren;}
269  // Get child particle
270  const ChainParticle* operator[](int i) const;
271  // Get parent particle
272  const ChainParticle* getParent() const;
273  // Get energy in parent frame
274  double E_parentFrame() const;
275  // Get particle ID
276  string getpID() const {return pID;}
277  // Print the decay chain (to cout)
278  void printChain() const;
279  // Get weight factor (see description of the weight variable)
280  double getWeight() const {return weight;}
281  // Get boost between lab frame and CoM frame of this particle
282  void getBoost(double& gamma, double& beta) const;
283  // Get pointer to decay table
284  const DecayTable* getDecayTable() const {return decayTable;}
285  // Destructor
286  ~ChainParticle();
287  private:
288  // Helper function for printChain()
289  bool printChain(int generation, vector<int> ancestry) const;
290  // How much the decay chain (to this point) should be weighted down due to
291  // missing/disabled decay channels. Only relevant for Monte Carlo generated chains.
292  double weight;
293  // Pointer to decay table
295  // Lorentz matrices
298  // 4-momentum in parent's rest frame
300  // Particle identifier
301  string pID;
302  // How many ancestors do I have?
304  // Has this particle been kept from decaying by an energy or chain length cut?
306  // Is this particle stable, or has it decayed to a final state consisting entirely of stable particles?
307  // If so, this particle is considered a final state, and any children must be extracted manually after collecting the final states.
308  // This is to avoid having final states that can't be treated as free particles (such as quarks).
310  // Number of child particles
312  // Pointers to parent and child particles
314  vector<ChainParticle*> children;
315  // Function for updating the Lorentz boost matrices according to a new 4-momentum.
316  void update(vec4 &ip_parent);
317  // Constructor used by member functions during chain generation.
318  ChainParticle(const vec4 &pp, double m, double weight, const DecayTable *dc, ChainParticle *parent, int chainGeneration, string pID);
319  // Disable copy constructor and assignment operator. These would cause mayhem.
321  ChainParticle & operator=(const ChainParticle&);
322  };
323  typedef std::vector<const Gambit::DarkBit::DecayChain::ChainParticle*> ChainParticleVector;
324 
325  // Container for passing around ChainParticle objects.
327  {
329  ChainContainer(shared_ptr<ChainParticle> ch) : chain(ch) {}
330  shared_ptr<const ChainParticle> chain;
331  };
332 
333  }
334  }
335 }
336 
337 #endif // defined __decay_chain_hpp__
double invariantMass(const vec4 &a, const vec4 &b)
unordered_map< const TH_Channel *, bool > endpointFlags
vector< const TH_Channel * > disabledDecays
void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat)
START_MODEL beta
Definition: Axions.hpp:35
void boostMatrixParentFrame(mat4 &mat, vec4 &p_parent)
vec4 operator+(const vec4 &x, const vec4 &y)
unordered_map< string, DecayTableEntry > table
vec4 Ep4vec(const vec3 p, double m)
bool randInit
START_MODEL v0
vec4 lorentzBoost(const vec4 &inVec, const vec3 &beta_xyz)
bool useForcedTotalWidth
START_MODEL b
Definition: demo.hpp:235
vec3 operator/(const vec3 &y, double x)
Definition: decay_chain.cpp:79
vec4(double v0, double v1, double v2, double v3)
Definition: decay_chain.hpp:83
static double draw()
Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
vector< const TH_Channel * > enabledDecays
double operator[](int i) const
Definition: decay_chain.hpp:90
A container holding all annihilation and decay initial states relevant for DarkBit.
A threadsafe interface to the STL random number generators.
vec3 operator*(double x, const vec3 &y)
Definition: decay_chain.cpp:61
vec4 p_parentFrame(const vec4 &inVec, const vec4 &p_parent)
All data on a single annihilation or decay channel, e.g.
DecayTableEntry()
const DecayTable * getDecayTable() const
double totalWidth
DecayTableEntry(string pID, double m, bool stable)
double forcedTotalWidth
ChainContainer(shared_ptr< ChainParticle > ch)
double enabledWidth
double operator[](int i) const
Definition: decay_chain.hpp:64
Channel container Object containing tabularized yields for particle decay and two-body final states...
string pID
ostream & operator<<(ostream &os, const vec3 &v)
Definition: decay_chain.cpp:88
double dot(const vec3 &a, const vec3 &b)
Definition: decay_chain.cpp:93
const double m
double invisibleWidth
std::vector< const Gambit::DarkBit::DecayChain::ChainParticle * > ChainParticleVector
vector< double > randLims
shared_ptr< const ChainParticle > chain
TODO: see if we can use this one:
Definition: Analysis.hpp:33
bool stable
vec3(double v0, double v1, double v2)
Definition: decay_chain.hpp:54