32 using std::ostringstream;
90 os << v[0] <<
", " << v[1] <<
", " << v[2];
95 return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
142 "(" << v[0] <<
", " << v[1] <<
", " << v[2] <<
", " << v[3] <<
")";
147 return a[0]*b[0]-a[1]*b[1]-a[2]*b[2]-a[3]*b[3];
151 double E = sqrt(
dot(p,p)+m*m);
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)
216 out[i] += m.
vals[i][j]*v[j];
228 for(
int k=0; k<4; k++)
242 os << m.
vals[i][j] <<
", ";
268 while(r1*r1+r2*r2 >=1.0);
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);
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);
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];
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);
305 return lorentz*inVec;
309 vec3 beta_xyz = -p_parent.
xyz()/p_parent[0];
314 vec3 beta_xyz = -p_parent.
xyz()/p_parent[0];
315 double gamma = p_parent[0]/m;
321 return sqrt(
dot(tmp,tmp));
334 "Generating rand table at runtime. This should have happened\n" 335 "during initialization, and might not be threadsafe here.");
338 vector<double>::const_iterator pos = upper_bound(randLims.begin(),
339 randLims.end(),pick);
340 return pos - randLims.begin();
344 if(!hasEnabledDecays())
return false;
346 int idx = findChannelIdx(pick);
347 decay = enabledDecays[idx];
354 for(vector<const TH_Channel*>::const_iterator
355 it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
358 randLims.push_back(tmp);
366 for(vector<const TH_Channel*>::const_iterator
367 it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
372 for(vector<const TH_Channel*>::const_iterator
373 it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
377 totalWidth += invisibleWidth;
378 if(randInit) generateRandTable();
382 for(vector<const TH_Channel*>::const_iterator
383 it = enabledDecays.begin(); it != enabledDecays.end(); ++it)
385 if((*it) == in)
return true;
391 for(vector<const TH_Channel*>::const_iterator
392 it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
394 if((*it) == in)
return true;
400 if(isEnabled(in))
return true;
401 if(isDisabled(in))
return true;
409 "Trying to add decay with >2 final states to DecayTableEntry.\n" 410 "Channel added as disabled.");
414 enabledDecays.push_back(in);
420 disabledDecays.push_back(in);
425 invisibleWidth = width;
434 "Trying to enable decay with >2 final states in decay chain.\n" 439 for(vector<const TH_Channel*>::iterator
440 it = disabledDecays.begin(); it != disabledDecays.end(); ++it)
445 disabledDecays.erase(it);
448 if(!found)
return false;
449 enabledDecays.push_back(in);
452 if(randInit) generateRandTable();
458 for(vector<const TH_Channel*>::iterator it = enabledDecays.begin();
459 it != enabledDecays.end(); ++it)
464 enabledDecays.erase(it);
467 if(!found)
return false;
468 disabledDecays.push_back(in);
471 if(randInit) generateRandTable();
476 return enabledWidth/getTotalWidth();
480 useForcedTotalWidth = enabled;
481 forcedTotalWidth = width;
485 return useForcedTotalWidth ? forcedTotalWidth : totalWidth;
489 return (enabledDecays.size()>0);
500 std::cout <<
"Importing CascadeMC DecayTable from process catalog..." << std::endl;
502 set<string> finalStates;
504 for(vector<TH_Process>::const_iterator it = cat.
processList.begin();
508 if(it->isAnnihilation)
continue;
510 string pID = it->particle1ID;
512 bool stable = ((it->channelList).size()<1);
514 std::cout <<
"Importing " << pID << std::endl;
516 if(disabledList.count(pID)==1) stable =
true;
528 for(vector<TH_Channel>::const_iterator it2 = (
529 it->channelList).begin(); it2 != (it->channelList).end(); ++it2)
533 if((it2->nFinalStates) > 2)
541 for(vector<string>::const_iterator
542 it3 = (it2->finalStateIDs).begin();
543 it3 != (it2->finalStateIDs).end(); ++it3)
545 finalStates.insert(*it3);
552 "Unstable particle "+pID+
" has no available decay channels. Treating it as stable in cascade decays.");
559 for(unordered_map<string,DecayTableEntry>::iterator it = table.begin();
560 it != table.end(); ++it)
563 for(vector<const TH_Channel*>::const_iterator
564 it2 = (it->second.enabledDecays).begin();
565 it2 != (it->second.enabledDecays).end(); ++it2)
567 bool isEndpoint =
true;
569 for(vector<string>::const_iterator
570 it3 = ((*it2)->finalStateIDs).begin();
571 it3 != ((*it2)->finalStateIDs).end(); ++it3)
577 if(!table[*it3].stable)
584 it->second.endpointFlags[*it2] = isEndpoint;
586 it->second.generateRandTable();
590 for(set<string>::iterator it = finalStates.begin();
591 it != finalStates.end(); ++it)
596 addEntry(*it,m,
true);
600 std::cout <<
"...done" << std::endl;
605 return table.find(index) != table.end();
614 table.insert ( pair<string,DecayTableEntry>(pID,entry) );
621 ans=(table.at(pID)).randomDecay(decay);
623 catch (std::out_of_range& e)
626 "No partcile "+pID+
" in decay table."));
637 catch (std::out_of_range& e)
640 "No partcile "+i+
" in decay table."));
646 return ch->
genRate->bind()->eval();
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)
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)
671 for(vector<string>::const_iterator
672 it3 = ((*it2)->finalStateIDs).begin();
673 it3 != ((*it2)->finalStateIDs).end(); ++it3)
675 std::cout << *it3 <<
", ";
677 std::cout <<
"Width: " << getWidth(*it2) << endl;
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)
685 for(vector<string>::const_iterator
686 it3 = ((*it2)->finalStateIDs).begin();
687 it3 != ((*it2)->finalStateIDs).end(); ++it3)
689 std::cout << *it3 <<
", ";
691 std::cout <<
"Width: " << getWidth(*it2) << endl;
693 std::cout << std::endl;
705 m((*dc)[pID].m), weight(1), decayTable(dc), pID(pID),
706 chainGeneration(0), abortedDecay(false), isEndpoint(false),
707 nChildren(0), parent(NULL)
718 "Overwriting existing decay in decay chain.");
728 and ((Emin < 0) or (
E_Lab()> Emin)) )
735 "Unable to pick allowed decay for "+
pID+
". Keeping particle stable.");
743 err =
"Invalid decay channel in decay table.\n";
744 err+=
"N!=2 body decays are currently not supported in cascade decays.\n";
756 "Kinematically impossible decay in decay chain:\n" <<
759 "Please check your process catalog." << endl;
760 err <<
"Relevant particle masses: " <<
m <<
" -> " << m1 <<
" + " << m2;
763 const double &Etot =
m;
764 double E1 = 0.5*(Etot*Etot+m1*m1-m2*m2)/Etot;
766 double abs_p = sqrt(E1*E1-m1*m1);
769 vec4 p1(E1, abs_p*dir);
770 vec4 p2(E2,-abs_p*dir);
772 double wt =
weight*(*decayTable)[
pID].getEnabledBranching();
787 children[i]->generateDecayChainMC(maxSteps, Emin);
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);
839 &endpointStates,
bool includeAborted,
string ipID)
const 843 if(includeAborted && ((ipID==
"") || (ipID==
pID)))
844 endpointStates.push_back(
this);
850 for(vector<ChainParticle*>::const_iterator it=
children.begin();
853 (*it)->collectEndpointStates(endpointStates,includeAborted,
pID);
858 endpointStates.push_back(
this);
869 "Trying to access non-existing decay chain entry."));
883 std::cout <<
"*********************" << endl;
884 std::cout <<
"Decay chain printout:" << endl;
885 std::cout <<
"---------------------" << endl;
887 std::cout <<
"0 " <<
pID <<
", p = " <<
p_Lab() <<
888 ", Weight: " <<
weight << endl;
889 std::cout <<
"---------------------" << endl;
896 std::cout <<
"Generation " << gen <<
":" << endl;
900 vector<int> ancestry;
901 ancestry.push_back(0);
902 ancestry.push_back(i);
903 bool more =
children[i]->printChain(gen,ancestry);
907 std::cout <<
"---------------------" << endl;
920 vector<int> ancestry2 = ancestry;
921 ancestry2.push_back(i);
926 for(vector<int>::const_iterator it=ancestry.begin();
927 it!=ancestry.end(); ++it)
929 std::cout << *it <<
" ";
931 std::cout <<
pID <<
", p = " <<
p_Lab() <<
", Weight: " <<
weight << endl;
938 gamma = b.
vals[0][0];
double invariantMass(const vec4 &a, const vec4 &b)
vector< ChainParticle * > children
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
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.
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)
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)
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.
const DecayTable * decayTable
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)
void update(vec4 &ip_parent)
A container holding all annihilation and decay initial states relevant for DarkBit.
vec3 operator*(double x, const vec3 &y)
bool hasEntry(string) const
static double getWidth(const TH_Channel *ch)
vec4 p_to_Lab(const vec4 &p) const
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)
double dot(const vec3 &a, const vec3 &b)
double E_parentFrame() const
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:
double E(const double x, const double y)
const ChainParticle * getParent() const
const DecayTableEntry & operator[](string i) const