gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
decay_chain.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
23 
24 //#define DARKBIT_DEBUG
25 
26 namespace Gambit
27 {
28  namespace DarkBit
29  {
30  namespace DecayChain
31  {
32  using std::ostringstream;
33  using std::set;
34  using std::endl;
35  using std::pair;
36 
37  // *********************************************
38  // 3-vector related
39  // *********************************************
40 
41  double vec3::length() const
42  {
43  return sqrt(vals[0]*vals[0]+vals[1]*vals[1]+vals[2]*vals[2]);
44  }
46  {
47  double l = length();
48  for(int i=0;i<3;i++)
49  {
50  vals[i] /= l;
51  }
52  }
53  void vec3::normalize(const double len)
54  {
55  double l = length();
56  for(int i=0;i<3;i++)
57  {
58  vals[i] *= len/l;
59  }
60  }
61  vec3 operator* (double x, const vec3 &y)
62  {
63  vec3 retVec = y;
64  for(int i=0;i<3;i++)
65  {
66  retVec[i] *= x;
67  }
68  return retVec;
69  }
70  vec3 operator* (const vec3 &y, double x)
71  {
72  vec3 retVec = y;
73  for(int i=0;i<3;i++)
74  {
75  retVec[i] *= x;
76  }
77  return retVec;
78  }
79  vec3 operator/ (const vec3 &y, double x)
80  {
81  vec3 retVec = y;
82  for(int i=0;i<3;i++)
83  {
84  retVec[i] /= x;
85  }
86  return retVec;
87  }
88  std::ostream& operator<<(std::ostream& os, const vec3& v)
89  {
90  os << v[0] << ", " << v[1] << ", " << v[2];
91  return os;
92  }
93  double dot(const vec3 &a, const vec3 &b)
94  {
95  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
96  }
97 
98 
99  // *********************************************
100  // 4-vector related
101  // *********************************************
102 
103  vec4 operator* (double x, const vec4 &y)
104  {
105  vec4 retVec = y;
106  for(int i=0;i<4;i++)
107  {
108  retVec[i] *= x;
109  }
110  return retVec;
111  }
112  vec4 operator* (const vec4 &y, double x)
113  {
114  vec4 retVec = y;
115  for(int i=0;i<4;i++)
116  {
117  retVec[i] *= x;
118  }
119  return retVec;
120  }
121  vec4 operator+ (const vec4 &x, const vec4 &y)
122  {
123  vec4 retVec = x;
124  for(int i=0;i<4;i++)
125  {
126  retVec[i] += y[i];
127  }
128  return retVec;
129  }
130  vec4 operator- (const vec4 &x, const vec4 &y)
131  {
132  vec4 retVec = x;
133  for(int i=0;i<4;i++)
134  {
135  retVec[i] -= y[i];
136  }
137  return retVec;
138  }
139  std::ostream& operator<<(std::ostream& os, const vec4& v)
140  {
141  os <<
142  "(" << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] <<")";
143  return os;
144  }
145  double dot(const vec4 &a, const vec4 &b)
146  {
147  return a[0]*b[0]-a[1]*b[1]-a[2]*b[2]-a[3]*b[3];
148  }
149  vec4 Ep4vec(const vec3 p, double m)
150  {
151  double E = sqrt(dot(p,p)+m*m);
152  return vec4(E,p);
153  }
154 
155 
156  // *********************************************
157  // 4x4 matrix related
158  // *********************************************
159 
161  double v00, double v01, double v02, double v03,
162  double v10, double v11, double v12, double v13,
163  double v20, double v21, double v22, double v23,
164  double v30, double v31, double v32, double v33)
165  {
166  vals[0][0] = v00;
167  vals[0][1] = v01;
168  vals[0][2] = v02;
169  vals[0][3] = v03;
170  vals[1][0] = v10;
171  vals[1][1] = v11;
172  vals[1][2] = v12;
173  vals[1][3] = v13;
174  vals[2][0] = v20;
175  vals[2][1] = v21;
176  vals[2][2] = v22;
177  vals[2][3] = v23;
178  vals[3][0] = v30;
179  vals[3][1] = v31;
180  vals[3][2] = v32;
181  vals[3][3] = v33;
182  }
183  mat4::mat4(double v)
184  {
185  vals[0][0] = v;
186  vals[0][1] = v;
187  vals[0][2] = v;
188  vals[0][3] = v;
189  vals[1][0] = v;
190  vals[1][1] = v;
191  vals[1][2] = v;
192  vals[1][3] = v;
193  vals[2][0] = v;
194  vals[2][1] = v;
195  vals[2][2] = v;
196  vals[2][3] = v;
197  vals[3][0] = v;
198  vals[3][1] = v;
199  vals[3][2] = v;
200  vals[3][3] = v;
201  }
203  {
204  return mat4(1,0,0,0,
205  0,1,0,0,
206  0,0,1,0,
207  0,0,0,1);
208  }
209  vec4 operator* (const mat4 &m, const vec4 &v)
210  {
211  vec4 out(0);
212  for(int i=0;i<4;i++)
213  {
214  for(int j=0;j<4;j++)
215  {
216  out[i] += m.vals[i][j]*v[j];
217  }
218  }
219  return out;
220  }
221  mat4 operator* (const mat4 &m1, const mat4 &m2)
222  {
223  mat4 out(0);
224  for(int i=0;i<4;i++)
225  {
226  for(int j=0;j<4;j++)
227  {
228  for(int k=0; k<4; k++)
229  {
230  out.vals[i][j] += m1.vals[i][k]*m2.vals[k][j];
231  }
232  }
233  }
234  return out;
235  }
236  std::ostream& operator<<(std::ostream& os, const mat4& m)
237  {
238  for(int i=0;i<4;i++)
239  {
240  for(int j=0;j<4;j++)
241  {
242  os << m.vals[i][j] << ", ";
243  }
244  if(i!=3)
245  os << endl;
246  }
247  return os;
248  }
249 
250 
251  // *********************************************
252  // Utility functions
253  // *********************************************
254 
255  double rand_m1_1()
256  {
257  return -1.0 + 2.0 * rand_0_1();
258  }
260  {
261  // Marsaglia 1972 rejection method
262  double r1,r2;
263  do
264  {
265  r1 = rand_m1_1();
266  r2 = rand_m1_1();
267  }
268  while(r1*r1+r2*r2 >=1.0);
269  vec3 v;
270  v[0] = 2.0*r1*sqrt(1-r1*r1-r2*r2);
271  v[1] = 2.0*r2*sqrt(1-r1*r1-r2*r2);
272  v[2] = 1.0-2.0*(r1*r1+r2*r2);
273  return v;
274  }
275  void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat)
276  {
277  double b = beta_xyz.length();
278  double bm2 = b==0 ? 0 : 1.0/(b*b);
279  const double &bx = beta_xyz[0];
280  const double &by = beta_xyz[1];
281  const double &bz = beta_xyz[2];
282  double g = 1.0/sqrt(1-b*b);
283  mat = mat4( g, -g*bx, -g*by, -g*bz,
284  -g*bx, 1+(g-1)*bx*bx*bm2, (g-1)*bx*by*bm2, (g-1)*bx*bz*bm2,
285  -g*by, (g-1)*by*bx*bm2, 1+(g-1)*by*by*bm2, (g-1)*by*bz*bm2,
286  -g*bz, (g-1)*bz*bx*bm2, (g-1)*bz*by*bm2, 1+(g-1)*bz*bz*bm2);
287  }
288  void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat, double gamma)
289  {
290  double b = beta_xyz.length();
291  double bm2 = b==0 ? 0 : 1.0/(b*b);
292  const double &bx = beta_xyz[0];
293  const double &by = beta_xyz[1];
294  const double &bz = beta_xyz[2];
295  double &g = gamma;
296  mat = mat4( g, -g*bx, -g*by, -g*bz,
297  -g*bx, 1+(g-1)*bx*bx*bm2, (g-1)*bx*by*bm2, (g-1)*bx*bz*bm2,
298  -g*by, (g-1)*by*bx*bm2, 1+(g-1)*by*by*bm2, (g-1)*by*bz*bm2,
299  -g*bz, (g-1)*bz*bx*bm2, (g-1)*bz*by*bm2, 1+(g-1)*bz*bz*bm2);
300  }
301  vec4 lorentzBoost(const vec4 &inVec, const vec3 &beta_xyz)
302  {
303  mat4 lorentz;
304  lorentzMatrix(beta_xyz, lorentz);
305  return lorentz*inVec;
306  }
307  vec4 p_parentFrame(const vec4 &inVec, const vec4 &p_parent)
308  {
309  vec3 beta_xyz = -p_parent.xyz()/p_parent[0];
310  return lorentzBoost(inVec, beta_xyz);
311  }
312  void boostMatrixParentFrame(mat4 &mat, vec4 &p_parent, double m)
313  {
314  vec3 beta_xyz = -p_parent.xyz()/p_parent[0];
315  double gamma = p_parent[0]/m;
316  lorentzMatrix(beta_xyz, mat,gamma);
317  }
318  double invariantMass(const vec4 &a, const vec4 &b)
319  {
320  vec4 tmp = a+b;
321  return sqrt(dot(tmp,tmp));
322  }
323 
324 
325  // *********************************************
326  // DecayTableEntry functions
327  // *********************************************
328 
329  int DecayTableEntry::findChannelIdx(double pick) const
330  {
331  if(!randInit)
332  {
334  "Generating rand table at runtime. This should have happened\n"
335  "during initialization, and might not be threadsafe here.");
336  generateRandTable();
337  }
338  vector<double>::const_iterator pos = upper_bound(randLims.begin(),
339  randLims.end(),pick);
340  return pos - randLims.begin();
341  }
342  bool DecayTableEntry::randomDecay(const TH_Channel* &decay) const
343  {
344  if(!hasEnabledDecays()) return false;
345  double pick = rand_0_1();
346  int idx = findChannelIdx(pick);
347  decay = enabledDecays[idx];
348  return true;
349  }
351  {
352  randLims.clear();
353  double tmp=0;
354  for(vector<const TH_Channel*>::const_iterator
355  it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
356  {
357  tmp += DecayTable::getWidth(*it)/enabledWidth;
358  randLims.push_back(tmp);
359  }
360  randInit=true;
361  }
363  {
364  totalWidth = 0;
365  enabledWidth = 0;
366  for(vector<const TH_Channel*>::const_iterator
367  it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
368  {
369  enabledWidth += DecayTable::getWidth(*it);
370  totalWidth += DecayTable::getWidth(*it);
371  }
372  for(vector<const TH_Channel*>::const_iterator
373  it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
374  {
375  totalWidth += DecayTable::getWidth(*it);
376  }
377  totalWidth += invisibleWidth;
378  if(randInit) generateRandTable();
379  }
381  {
382  for(vector<const TH_Channel*>::const_iterator
383  it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
384  {
385  if((*it) == in) return true;
386  }
387  return false;
388  }
390  {
391  for(vector<const TH_Channel*>::const_iterator
392  it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
393  {
394  if((*it) == in) return true;
395  }
396  return false;
397  }
399  {
400  if(isEnabled(in)) return true;
401  if(isDisabled(in)) return true;
402  return false;
403  }
405  {
406  if(in->nFinalStates > 2)
407  {
409  "Trying to add decay with >2 final states to DecayTableEntry.\n"
410  "Channel added as disabled.");
411  addDisabled(in);
412  return;
413  }
414  enabledDecays.push_back(in);
415  totalWidth += DecayTable::getWidth(in);
416  enabledWidth += DecayTable::getWidth(in);
417  }
419  {
420  disabledDecays.push_back(in);
421  totalWidth += DecayTable::getWidth(in);
422  }
424  {
425  invisibleWidth = width;
426  totalWidth += width;
427  }
429  {
430  // N>2 body decays currently not possible.
431  if(in->nFinalStates > 2)
432  {
434  "Trying to enable decay with >2 final states in decay chain.\n"
435  "Request ignored.");
436  return false;
437  }
438  bool found = false;
439  for(vector<const TH_Channel*>::iterator
440  it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
441  {
442  if((*it) == in)
443  {
444  found = true;
445  disabledDecays.erase(it);
446  }
447  }
448  if(!found) return false;
449  enabledDecays.push_back(in);
450  enabledWidth += DecayTable::getWidth(in);
451  // Re-generate Monte Carlo table if necessary
452  if(randInit) generateRandTable();
453  return true;
454  }
456  {
457  bool found = false;
458  for(vector<const TH_Channel*>::iterator it = enabledDecays.begin();
459  it != enabledDecays.end(); ++it)
460  {
461  if((*it) == in)
462  {
463  found = true;
464  enabledDecays.erase(it);
465  }
466  }
467  if(!found) return false;
468  disabledDecays.push_back(in);
469  enabledWidth -= DecayTable::getWidth(in);
470  // Re-generate Monte Carlo table if necessary
471  if(randInit) generateRandTable();
472  return true;
473  }
475  {
476  return enabledWidth/getTotalWidth();
477  }
478  void DecayTableEntry::forceTotalWidth(bool enabled, double width)
479  {
480  useForcedTotalWidth = enabled;
481  forcedTotalWidth = width;
482  }
484  {
485  return useForcedTotalWidth ? forcedTotalWidth : totalWidth;
486  }
488  {
489  return (enabledDecays.size()>0);
490  }
491 
492  // *********************************************
493  // DecayTable functions
494  // *********************************************
495 
497  const SimYieldTable &tab, set<string> disabledList)
498  {
499 #ifdef DARKBIT_DEBUG
500  std::cout << "Importing CascadeMC DecayTable from process catalog..." << std::endl;
501 #endif
502  set<string> finalStates;
503  // Register all decaying particles and their decays
504  for(vector<TH_Process>::const_iterator it = cat.processList.begin();
505  it != cat.processList.end(); ++it)
506  {
507  // Only interested in decay processes
508  if(it->isAnnihilation) continue;
509 
510  string pID = it->particle1ID;
511  double m = cat.getParticleProperty(pID).mass;
512  bool stable = ((it->channelList).size()<1);
513 #ifdef DARKBIT_DEBUG
514  std::cout << "Importing " << pID << std::endl;
515 #endif
516  if(disabledList.count(pID)==1) stable = true;
517  // If tabulated spectra exist for decays of this particle, consider
518  // it stable for the purpose of decay chain generation.
519  if(tab.hasAnyChannel(pID)) stable = true;
520 // if(!stable and (width <=0.0))
521 // {
522 // piped_warnings.request(LOCAL_INFO,
523 // "Unstable particle "+pID+" with zero width in decay table. Treating it as stable in cascade decays.");
524 // stable = true;
525 // }
526  // Create DecayTableEntry and insert decay channels
527  DecayTableEntry entry(pID,m,stable);
528  for(vector<TH_Channel>::const_iterator it2 = (
529  it->channelList).begin(); it2 != (it->channelList).end(); ++it2)
530  {
531  // N>2 body decays currently not supported, but should be included
532  // as disabled to get the correct total width.
533  if((it2->nFinalStates) > 2)
534  {
535  entry.addDisabled(&(*it2));
536  }
537  else
538  {
539  entry.addChannel(&(*it2));
540  }
541  for(vector<string>::const_iterator
542  it3 = (it2->finalStateIDs).begin();
543  it3 != (it2->finalStateIDs).end(); ++it3)
544  {
545  finalStates.insert(*it3);
546  }
547  }
548  entry.setInvisibleWidth(it->genRateMisc->bind()->eval());
549  if(!stable and entry.enabledDecays.size() == 0)
550  {
552  "Unstable particle "+pID+" has no available decay channels. Treating it as stable in cascade decays.");
553  entry.stable = true;
554  }
555  addEntry(pID,entry);
556  }
557  // Flag channels where all final final states are stable as endpoints.
558  // Loop over all particles
559  for(unordered_map<string,DecayTableEntry>::iterator it = table.begin();
560  it != table.end(); ++it)
561  {
562  // Loop over all decays
563  for(vector<const TH_Channel*>::const_iterator
564  it2 = (it->second.enabledDecays).begin();
565  it2 != (it->second.enabledDecays).end(); ++it2)
566  {
567  bool isEndpoint = true;
568  // Loop over all final states
569  for(vector<string>::const_iterator
570  it3 = ((*it2)->finalStateIDs).begin();
571  it3 != ((*it2)->finalStateIDs).end(); ++it3)
572  {
573  // All (and only) decaying particles are registered so far. If
574  // particle is registered, it's not stable.
575  if(hasEntry(*it3))
576  {
577  if(!table[*it3].stable)
578  {
579  isEndpoint = false;
580  break;
581  }
582  }
583  }
584  it->second.endpointFlags[*it2] = isEndpoint;
585  }
586  it->second.generateRandTable();
587  }
588  // Iterate over final states and register particles that have not
589  // already been registered (because they are considered stable).
590  for(set<string>::iterator it = finalStates.begin();
591  it != finalStates.end(); ++it)
592  {
593  if(!hasEntry(*it))
594  {
595  double m = cat.getParticleProperty(*it).mass;
596  addEntry(*it,m,true);
597  }
598  }
599 #ifdef DARKBIT_DEBUG
600  std::cout << "...done" << std::endl;
601 #endif
602  }
603  bool DecayTable::hasEntry(string index) const
604  {
605  return table.find(index) != table.end();
606  }
607  void DecayTable::addEntry(string pID, double m, bool stable)
608  {
609  table.insert (
610  pair<string,DecayTableEntry>(pID,DecayTableEntry(pID,m,stable)) );
611  }
612  void DecayTable::addEntry(string pID, DecayTableEntry entry)
613  {
614  table.insert ( pair<string,DecayTableEntry>(pID,entry) );
615  }
616  bool DecayTable::randomDecay(string pID, const TH_Channel* &decay) const
617  {
618  bool ans=false;
619  try
620  {
621  ans=(table.at(pID)).randomDecay(decay);
622  }
623  catch (std::out_of_range& e)
624  {
626  "No partcile "+pID+" in decay table."));
627  }
628  return ans;
629  }
631  {
632  const DecayTableEntry *ent = NULL;
633  try
634  {
635  ent= &table.at(i);
636  }
637  catch (std::out_of_range& e)
638  {
640  "No partcile "+i+" in decay table."));
641  }
642  return *ent;
643  }
645  {
646  return ch->genRate->bind()->eval();
647  }
649  {
650 #ifdef DARKBIT_DEBUG
651  std::cout << std::endl;
652  std::cout << "***********************" << endl;
653  std::cout << "CMC DecayTable printout" << endl;
654  std::cout << "***********************" << endl;
655  std::cout << std::endl;
656  for(unordered_map<string,DecayTableEntry>::const_iterator
657  it = table.begin(); it != table.end(); ++it)
658  {
659  std::cout << "Particle: " <<(it->first) << endl;
660  std::cout << "Set stable: " << (it->second).stable << endl;
661  std::cout << "Mass: " <<(it->second).m << endl;
662  std::cout << "Total width: " << (it->second.getTotalWidth())<< endl;
663  std::cout << "Enabled branching ratio: "
664  << (it->second.getEnabledBranching()) << endl;
665  std::cout << "Enabled decays:" << endl;
666  for(vector<const TH_Channel*>::const_iterator
667  it2 = (it->second.enabledDecays).begin();
668  it2 != (it->second.enabledDecays).end(); ++it2)
669  {
670  std::cout << " ";
671  for(vector<string>::const_iterator
672  it3 = ((*it2)->finalStateIDs).begin();
673  it3 != ((*it2)->finalStateIDs).end(); ++it3)
674  {
675  std::cout << *it3 << ", ";
676  }
677  std::cout << "Width: " << getWidth(*it2) << endl;
678  }
679  std::cout << "Disabled decays:" << endl;
680  for(vector<const TH_Channel*>::const_iterator
681  it2 = (it->second.disabledDecays).begin();
682  it2 != (it->second.disabledDecays).end(); ++it2)
683  {
684  std::cout << " ";
685  for(vector<string>::const_iterator
686  it3 = ((*it2)->finalStateIDs).begin();
687  it3 != ((*it2)->finalStateIDs).end(); ++it3)
688  {
689  std::cout << *it3 << ", ";
690  }
691  std::cout << "Width: " << getWidth(*it2) << endl;
692  }
693  std::cout << std::endl;
694  }
695 #endif
696  }
697 
698 
699  // *********************************************
700  // ChainParticle functions
701  // *********************************************
702 
704  vec3 ipLab, const DecayTable *dc, string pID) :
705  m((*dc)[pID].m), weight(1), decayTable(dc), pID(pID),
706  chainGeneration(0), abortedDecay(false), isEndpoint(false),
707  nChildren(0), parent(NULL)
708  {
709  p_parent=Ep4vec(ipLab,m);
712  }
713  void ChainParticle::generateDecayChainMC(int maxSteps, double Emin)
714  {
715  if(nChildren!=0)
716  {
718  "Overwriting existing decay in decay chain.");
719  cutChain();
720  }
721  // Stable particles flagged as endpoints
722  if((*decayTable)[pID].stable)
723  {
724  isEndpoint = true;
725  }
726  // Check whether or not to proceed with decay
727  else if( ((maxSteps < 0) or (int(chainGeneration) < maxSteps))
728  and ((Emin < 0) or (E_Lab()> Emin)) )
729  {
730  const TH_Channel *chn;
731  bool canDecay = decayTable->randomDecay(pID, chn);
732  if(!canDecay)
733  {
735  "Unable to pick allowed decay for "+ pID+". Keeping particle stable.");
736  abortedDecay = true;
737  return;
738  }
739  // Only 2-body decays are currently allowed
740  if((chn->nFinalStates) != 2)
741  {
742  string err;
743  err = "Invalid decay channel in decay table.\n";
744  err+= "N!=2 body decays are currently not supported in cascade decays.\n";
746  return;
747  }
748  nChildren = 2; // chn->nFinalStates;
749  // Kinematics for 2-body decays
750  double m1 = (*decayTable)[(chn->finalStateIDs)[0]].m;
751  double m2 = (*decayTable)[(chn->finalStateIDs)[1]].m;
752  if(m1+m2>m)
753  {
754  ostringstream err;
755  err <<
756  "Kinematically impossible decay in decay chain:\n" <<
757  pID << "-> " <<
758  ((chn->finalStateIDs)[0]) << ", " << ((chn->finalStateIDs)[1]) << "\n" <<
759  "Please check your process catalog." << endl;
760  err << "Relevant particle masses: " << m << " -> " << m1 << " + " << m2;
761  throw(Piped_exceptions::description(LOCAL_INFO, err.str()));
762  }
763  const double &Etot = m;
764  double E1 = 0.5*(Etot*Etot+m1*m1-m2*m2)/Etot;
765  double E2 = Etot-E1;
766  double abs_p = sqrt(E1*E1-m1*m1);
767  // Assume isotropic decay in CM frame (no spin correlations).
768  vec3 dir = randOnSphere();
769  vec4 p1(E1, abs_p*dir);
770  vec4 p2(E2,-abs_p*dir);
771  // Weight from not including all possible decay channels
772  double wt = weight*(*decayTable)[pID].getEnabledBranching();
773  children.push_back(new ChainParticle(p1, m1, wt, decayTable, this,
774  chainGeneration+1, chn->finalStateIDs[0]));
775  children.push_back(new ChainParticle(p2, m2, wt, decayTable, this,
776  chainGeneration+1, chn->finalStateIDs[1]));
777  // Reached chain endpoint. Don't attempt further decays
778  if((*decayTable)[pID].endpointFlags.at(chn))
779  {
780  isEndpoint = true;
781  }
782  // Continue chain from child links.
783  else
784  {
785  for(int i=0;i<nChildren;i++)
786  {
787  children[i]->generateDecayChainMC(maxSteps, Emin);
788  }
789  }
790  }
791  else
792  {
793  abortedDecay = true;
794  }
795  }
797  {
798  if(nChildren==2)
799  {
800  double m1 = children[0]->m;
801  vec3 dir = randOnSphere();
802  double E1 = children[0]->p_parent[0];
803  double E2 = children[1]->p_parent[0];
804  double abs_p = sqrt(E1*E1-m1*m1);
805  vec4 p1(E1, abs_p*dir);
806  vec4 p2(E2,-abs_p*dir);
807  children[0]->update(p1);
808  children[1]->update(p2);
809  children[0]->reDrawAngles();
810  children[1]->reDrawAngles();
811  }
812  }
814  {
815  for(int i=0;i<nChildren; i++) delete children[i];
816  children.clear();
817  nChildren = 0;
818  }
820  {
821  return boostToLabFrame*p;
822  }
824  {
825  if(chainGeneration==0) return p_parent;
827  }
828  double ChainParticle::E_Lab() const
829  {
830  if(chainGeneration==0) return p_parent[0];
831  double E = 0;
832  for(int j=0;j<4;j++)
833  {
834  E += (parent->boostToLabFrame).vals[0][j]*p_parent[j];
835  }
836  return E;
837  }
838  void ChainParticle::collectEndpointStates(vector<const ChainParticle*>
839  &endpointStates, bool includeAborted, string ipID) const
840  {
841  if(abortedDecay)
842  {
843  if(includeAborted && ((ipID=="") || (ipID==pID)))
844  endpointStates.push_back(this);
845  }
846  else
847  {
848  if(nChildren!=0 and !isEndpoint)
849  {
850  for(vector<ChainParticle*>::const_iterator it=children.begin();
851  it!=children.end(); ++it)
852  {
853  (*it)->collectEndpointStates(endpointStates,includeAborted,pID);
854  }
855  }
856  else if((ipID=="") or (ipID==pID) or isEndpoint)
857  {
858  endpointStates.push_back(this);
859  }
860  }
861  }
863  {
864  if(i<nChildren)
865  return children[i];
866  else
867  {
869  "Trying to access non-existing decay chain entry."));
870  return this;
871  }
872  }
874  {
875  return parent;
876  }
878  {
879  return p_parent[0];
880  }
882  {
883  std::cout << "*********************" << endl;
884  std::cout << "Decay chain printout:" << endl;
885  std::cout << "---------------------" << endl;
886  std::cout << "Generation " << chainGeneration << ":" << endl;
887  std::cout << "0 " << pID << ", p = " << p_Lab() <<
888  ", Weight: " << weight << endl;
889  std::cout << "---------------------" << endl;
890  if(nChildren>0)
891  {
892  bool run = false;
893  int gen = chainGeneration+1;
894  do
895  {
896  std::cout << "Generation " << gen <<":" << endl;
897  run= false;
898  for(int i=0;i<nChildren;i++)
899  {
900  vector<int> ancestry;
901  ancestry.push_back(0);
902  ancestry.push_back(i);
903  bool more = children[i]->printChain(gen,ancestry);
904  run = more || run;
905  }
906  gen++;
907  std::cout << "---------------------" << endl;
908  }
909  while(run);
910  }
911  }
912  bool ChainParticle::printChain(int generation, vector<int> ancestry) const
913  {
914  if(generation < chainGeneration) return false;
915  if(generation > chainGeneration)
916  {
917  bool more = false;
918  for(int i=0;i<nChildren;i++)
919  {
920  vector<int> ancestry2 = ancestry;
921  ancestry2.push_back(i);
922  if(children[i]->printChain(generation,ancestry2)) more = true;
923  }
924  return more;
925  }
926  for(vector<int>::const_iterator it=ancestry.begin();
927  it!=ancestry.end(); ++it)
928  {
929  std::cout << *it << " ";
930  }
931  std::cout << pID << ", p = " << p_Lab() << ", Weight: " << weight << endl;
932  if(nChildren>0) return true;
933  return false;
934  }
935  void ChainParticle::getBoost(double& gamma, double& beta) const
936  {
937  const mat4& b=boostToLabFrame;
938  gamma = b.vals[0][0];
939  beta = sqrt(b.vals[0][1]*b.vals[0][1]+b.vals[0][2]*b.vals[0][2]+
940  b.vals[0][3]*b.vals[0][3])/gamma;
941  }
943  {
944  for(int i=0;i<nChildren; i++) delete children[i];
945  }
946  void ChainParticle::update(vec4 &ip_parent)
947  {
948  p_parent = ip_parent;
951  }
952  ChainParticle::ChainParticle(const vec4 &pp, double m, double weight,
954  string pID) :
955  m(m), weight(weight), decayTable(dc), p_parent(pp), pID(pID),
956  chainGeneration(chainGeneration), abortedDecay(false),
957  isEndpoint(false), nChildren(0), parent(parent)
958  {
961  }
962 
963  } // namespace DecayChain
964  } // namespace DarkBit
965 } // namespace Gambit
double invariantMass(const vec4 &a, const vec4 &b)
double getTotalWidth() const
TH_ParticleProperty getParticleProperty(str) const
Retrieve properties of a given particle involved in one or more processes in this catalog...
void addDisabled(const TH_Channel *in)
void setInvisibleWidth(double width)
void addEntry(string pID, double m, bool stable)
void lorentzMatrix(const vec3 &beta_xyz, mat4 &mat)
std::pair< std::string, std::string > description
Definition: exceptions.hpp:306
START_MODEL beta
Definition: Axions.hpp:36
void collectEndpointStates(vector< const ChainParticle *> &endpointStates, bool includeAborted, string ipID="") const
double getEnabledBranching() const
void boostMatrixParentFrame(mat4 &mat, vec4 &p_parent)
vec4 operator+(const vec4 &x, const vec4 &y)
void getBoost(double &gamma, double &beta) const
vec4 Ep4vec(const vec3 p, double m)
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:547
#define LOCAL_INFO
Definition: local_info.hpp:34
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void run(bool, bool, bool, int, double, double, int, int, int, int, int, double, char[], int, int[], bool, bool, bool, bool, double, int, double(*)(double *, int, int, void *), void(*)(int, int, int, double *, double *, double *, double, double, double, void *), void *)
vec4 lorentzBoost(const vec4 &inVec, const vec3 &beta_xyz)
START_MODEL b
Definition: demo.hpp:270
bool randomDecay(string pID, const TH_Channel *&decay) const
void forceTotalWidth(bool enabled, double width)
bool randomDecay(const TH_Channel *&decay) const
vec3 operator/(const vec3 &y, double x)
Definition: decay_chain.cpp:79
int findChannelIdx(double pick) const
bool hasEnabledDecays() const
double mass
Particle mass (GeV)
bool disableDecay(const TH_Channel *in)
unsigned int nFinalStates
Number of final state particles in this channel.
bool isDisabled(const TH_Channel *in) const
daFunk::Funk genRate
Energy dependence of final state particles. Includes v_rel ("v") as last argument in case of annihila...
vector< const TH_Channel * > enabledDecays
ChainParticle(vec3 ipLab, const DecayTable *dc, string pID)
bool enableDecay(const TH_Channel *in)
A container holding all annihilation and decay initial states relevant for DarkBit.
Header file that includes all GAMBIT headers required for a module source file.
vec3 operator*(double x, const vec3 &y)
Definition: decay_chain.cpp:61
static double getWidth(const TH_Channel *ch)
vec4 p_parentFrame(const vec4 &inVec, const vec4 &p_parent)
All data on a single annihilation or decay channel, e.g.
std::vector< std::string > finalStateIDs
Final state identifiers.
Rollcall header for module DarkBit.
bool hasAnyChannel(std::string p1) const
void addChannel(const TH_Channel *in)
void generateRandTable() const
Channel container Object containing tabularized yields for particle decay and two-body final states...
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 ChainParticle * operator[](int i) const
bool isRegistered(const TH_Channel *in) const
bool isEnabled(const TH_Channel *in) const
std::vector< TH_Process > processList
Vector of all processes in this catalog.
void generateDecayChainMC(int maxSteps, double Emin)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
void update()
bool stable
const ChainParticle * getParent() const
const DecayTableEntry & operator[](string i) const