Warning, /include/Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // KinematicsReconstructor.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 //
0010 // This is the implementation of the non-inlined templated member
0011 // functions of the KinematicsReconstructor class.
0012 //
0013 using namespace Herwig;
0014 using namespace ThePEG;
0015
0016 namespace {
0017 /**
0018 * find showering particle for hard branchings
0019 */
0020 tShowerParticlePtr SHOWERINGPARTICLE(HardBranchingPtr a) {
0021 return a->branchingParticle();
0022 }
0023
0024 /**
0025 * find showering particle for progenitors
0026 */
0027 tShowerParticlePtr SHOWERINGPARTICLE(ShowerProgenitorPtr a) {
0028 return a->progenitor();
0029 }
0030
0031
0032 /**
0033 * Return colour line progenitor pointer for ShowerProgenitor
0034 */
0035 template<typename Value>
0036 Ptr<ThePEG::ColourLine>::transient_pointer
0037 CL(Value a, unsigned int index=0) {
0038 return const_ptr_cast<ThePEG::tColinePtr>(SHOWERINGPARTICLE(a)->colourInfo()->colourLines()[index]);
0039 }
0040
0041 /**
0042 * Return progenitor colour line size for ShowerProgenitor
0043 */
0044 template<typename Value>
0045 unsigned int CLSIZE(Value a) {
0046 return SHOWERINGPARTICLE(a)->colourInfo()->colourLines().size();
0047 }
0048
0049 /**
0050 * Return anti-colour line progenitor pointer for ShowerProgenitor
0051 */
0052 template<typename Value>
0053 Ptr<ThePEG::ColourLine>::transient_pointer
0054 ACL(Value a, unsigned int index=0) {
0055 return const_ptr_cast<ThePEG::tColinePtr>(SHOWERINGPARTICLE(a)->colourInfo()->antiColourLines()[index]);
0056 }
0057
0058 /**
0059 * Return progenitor anti-colour line size for ShowerProgenitor
0060 */
0061 template<typename Value>
0062 unsigned int ACLSIZE(Value a) {
0063 return SHOWERINGPARTICLE(a)->colourInfo()->antiColourLines().size();
0064 }
0065 }
0066
0067 template<typename Value> void KinematicsReconstructor::
0068 findPartners(Value jet,set<Value> & done,
0069 const set<Value> & jets,
0070 vector<Value> & system) const {
0071 tShowerParticlePtr part=SHOWERINGPARTICLE(jet);
0072 unsigned int partNumColourLines = part->colourInfo()-> colourLines().size();
0073 unsigned int partNumAColourLines = part->colourInfo()->antiColourLines().size();
0074 for(typename set<Value>::const_iterator cit=jets.begin();cit!=jets.end();++cit) {
0075 if(done.find(*cit)!=done.end()||!SHOWERINGPARTICLE(*cit)->coloured())
0076 continue;
0077 bool isPartner = false;
0078 // one initial one final
0079 if(part->isFinalState()!=SHOWERINGPARTICLE(*cit)->isFinalState()) {
0080 //loop over all the colours of both
0081 for(unsigned int ix=0; ix<partNumColourLines; ++ix) {
0082 for(unsigned int jx=0; jx<CLSIZE(*cit); ++jx) {
0083 if(CL(jet,ix) && CL(jet,ix)==CL(*cit,jx)) {
0084 isPartner = true;
0085 break;
0086 }
0087 }
0088 if(isPartner) break;
0089 }
0090 if(!isPartner) {
0091 //loop over anti colours of both
0092 for(unsigned int ix=0; ix<partNumAColourLines; ++ix) {
0093 for(unsigned int jx=0; jx<ACLSIZE(*cit); ++jx) {
0094 if(ACL(jet,ix) && ACL(jet,ix)==ACL(*cit,jx)) {
0095 isPartner = true;
0096 break;
0097 }
0098 }
0099 if(isPartner) break;
0100 }
0101 }
0102 }
0103 // both in either initial or final state
0104 else {
0105 // loop over the colours of the first and the anti-colours of the other
0106 if(part->colourLine()) {
0107 for(unsigned int ix=0; ix<partNumColourLines; ++ix) {
0108 for(unsigned int jx=0; jx<ACLSIZE(*cit); ++jx) {
0109 if(CL(jet,ix) && CL(jet,ix)==ACL(*cit,jx)) {
0110 isPartner = true;
0111 break;
0112 }
0113 }
0114 if(isPartner) break;
0115 }
0116 }
0117 //loop over the anti-colours of the first and the colours of the other
0118 if(part->antiColourLine()&&!isPartner) {
0119 for(unsigned int ix=0; ix<partNumAColourLines; ++ix) {
0120 for(unsigned int jx=0; jx<CLSIZE(*cit); jx++) {
0121 if(ACL(jet,ix) && ACL(jet,ix)==CL(*cit,jx)) {
0122 isPartner = true;
0123 }
0124 }
0125 if(isPartner) break;
0126 }
0127 }
0128 }
0129 if(isPartner) {
0130 system.push_back(*cit);
0131 done.insert(*cit);
0132 findPartners(*cit,done,jets,system);
0133 continue;
0134 }
0135 // special for sources/sinks
0136 if(part->colourLine()) {
0137 if(part->colourLine()->sourceNeighbours().first) {
0138 tColinePair lines = part->colourLine()->sourceNeighbours();
0139 if(lines.first == CL(*cit) || lines.first == ACL(*cit) ||
0140 lines.second == CL(*cit) || lines.second == ACL(*cit) )
0141 isPartner = true;
0142 }
0143 if(part->colourLine()->sinkNeighbours().first) {
0144 tColinePair lines = part->colourLine()->sinkNeighbours();
0145 if(lines.first == CL(*cit) || lines.first == ACL(*cit) ||
0146 lines.second == CL(*cit) || lines.second == ACL(*cit) )
0147 isPartner = true;
0148 }
0149 }
0150 if(part->antiColourLine()) {
0151 if(part->antiColourLine()->sourceNeighbours().first) {
0152 tColinePair lines = part->antiColourLine()->sourceNeighbours();
0153 if(lines.first == CL(*cit) || lines.first == ACL(*cit) ||
0154 lines.second == CL(*cit) || lines.second == ACL(*cit) )
0155 isPartner = true;
0156 }
0157 if(part->antiColourLine()->sinkNeighbours().first) {
0158 tColinePair lines = part->antiColourLine()->sinkNeighbours();
0159 if(lines.first == CL(*cit) || lines.first == ACL(*cit) ||
0160 lines.second == CL(*cit) || lines.second == ACL(*cit) )
0161 isPartner = true;
0162 }
0163 }
0164 if(isPartner) {
0165 system.push_back(*cit);
0166 done.insert(*cit);
0167 findPartners(*cit,done,jets,system);
0168 }
0169 }
0170 }
0171
0172 template<typename Value >
0173 typename Herwig::ColourSinglet<Value>::VecType KinematicsReconstructor::
0174 identifySystems(set<Value> jets,
0175 unsigned int & nnun,unsigned int & nnii,unsigned int & nnif,
0176 unsigned int & nnf ,unsigned int & nni ) const {
0177 vector<ColourSinglet<Value> > systems;
0178 set<Value> done;
0179 for(typename set<Value>::const_iterator it=jets.begin();it!=jets.end();++it) {
0180 // if not treated create new system
0181 if(done.find(*it)!=done.end()) continue;
0182 done.insert(*it);
0183 systems.push_back(ColourSinglet<Value> (UNDEFINED,*it));
0184 if(!SHOWERINGPARTICLE(*it)->coloured()) continue;
0185 findPartners(*it,done,jets,systems.back().jets);
0186 }
0187 for(unsigned int ix=0;ix<systems.size();++ix) {
0188 unsigned int ni(0),nf(0);
0189 for(unsigned int iy=0;iy<systems[ix].jets.size();++iy) {
0190 if(SHOWERINGPARTICLE(systems[ix].jets[iy])->isFinalState()) ++nf;
0191 else ++ni;
0192 }
0193 // type
0194 // initial-initial
0195 if(ni==2&&nf==0) {
0196 systems[ix].type = II;
0197 ++nnii;
0198 }
0199 // initial only
0200 else if(ni==1&&nf==0) {
0201 systems[ix].type = I;
0202 ++nni;
0203 }
0204 // initial-final
0205 else if(ni==1&&nf>0) {
0206 systems[ix].type = IF;
0207 ++nnif;
0208 }
0209 // final only
0210 else if(ni==0&&nf>0) {
0211 systems[ix].type = F;
0212 ++nnf;
0213 }
0214 // otherwise unknown
0215 else {
0216 systems[ix].type = UNDEFINED;
0217 ++nnun;
0218 }
0219 }
0220 return systems;
0221 }
0222
0223 template<typename Value >
0224 void KinematicsReconstructor::combineFinalState(vector<ColourSinglet<Value> > & systems) const {
0225 // check that 1 particle final-state systems which can be combine
0226 bool canCombine(true);
0227 for(unsigned int ix=0;ix<systems.size();++ix) {
0228 if(systems[ix].type!=F) continue;
0229 if(systems[ix].jets.size()!=1) canCombine = false;
0230 }
0231 // return if can't combine
0232 if(!canCombine) return;
0233 // otherwise combine them
0234 vector<ColourSinglet<Value> > oldsystems=systems;
0235 systems.clear();
0236 ColourSinglet<Value> finalState;
0237 finalState.type = F;
0238 for(unsigned int ix=0;ix<oldsystems.size();++ix) {
0239 if(oldsystems[ix].type==F) {
0240 for(unsigned int iy=0;iy<oldsystems[ix].jets.size();++iy)
0241 finalState.jets.push_back(oldsystems[ix].jets[iy]);
0242 }
0243 else
0244 systems.push_back(oldsystems[ix]);
0245 }
0246 systems.push_back(finalState);
0247 }
