Warning, /include/ThePEG/Utilities/UtilityBase.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // UtilityBase.tcc is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG 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 namespace ThePEG {
0011
0012 template <typename PType>
0013 LorentzRotation UtilityBase::boostToCM(const pair<PType,PType> & pp) {
0014 typedef ParticleTraits<PType> Traits;
0015 LorentzRotation boost(-(Traits::momentum(pp.first) +
0016 Traits::momentum(pp.second)).boostVector());
0017 Traits::transform(pp.first, boost);
0018 Traits::transform(pp.second, boost);
0019 LorentzRotation rot;
0020 rot.rotateZ(-Traits::momentum(pp.first).phi());
0021 rot.rotateY(-Traits::momentum(pp.first).theta());
0022 rot.rotateZ(Traits::momentum(pp.first).phi());
0023 Traits::transform(pp.first, rot);
0024 Traits::transform(pp.second, rot);
0025 return rot*boost;
0026 }
0027
0028 template <typename PType>
0029 LorentzRotation UtilityBase::boostToCM(const Triplet<PType,PType,PType> & pt) {
0030 typedef ParticleTraits<PType> Traits;
0031 LorentzRotation rot = getBoostToCM(pt);
0032 Traits::transform(pt.first, rot);
0033 Traits::transform(pt.second, rot);
0034 Traits::transform(pt.third, rot);
0035 return rot;
0036 }
0037
0038 template <typename PType>
0039 LorentzRotation UtilityBase::getBoostToCM(const pair<PType,PType> & pp) {
0040 typedef ParticleTraits<PType> Traits;
0041 LorentzMomentum p1 = Traits::momentum(pp.first);
0042 Boost b = (p1 + Traits::momentum(pp.second)).boostVector();
0043 p1.boost(-b);
0044 LorentzRotation rot(-b);
0045 rot.rotateZ(-p1.phi());
0046 rot.rotateY(-p1.theta());
0047 rot.rotateZ(p1.phi());
0048 return rot;
0049 }
0050
0051 template <typename PType>
0052 LorentzRotation UtilityBase::
0053 getBoostToCM(const Triplet<PType,PType,PType> & pt) {
0054 typedef ParticleTraits<PType> Traits;
0055 LorentzMomentum p1 = Traits::momentum(pt.first);
0056 LorentzMomentum p2 = Traits::momentum(pt.second);
0057 Boost b = (p1 + p2 + Traits::momentum(pt.third)).boostVector();
0058 p1.boost(-b);
0059 LorentzRotation rot(-b);
0060 rot.rotateZ(-p1.phi());
0061 rot.rotateY(-p1.theta());
0062 p2.transform(rot);
0063 rot.rotateZ(-p2.phi());
0064 return rot;
0065 }
0066
0067 template <typename PType>
0068 LorentzRotation UtilityBase::
0069 getBoostFromCM(const Triplet<PType,PType,PType> & pt) {
0070 return getBoostToCM(pt).inverse();
0071 }
0072
0073 template <typename PType>
0074 LorentzRotation UtilityBase::getBoostFromCM(const pair<PType,PType> & pp) {
0075 typedef ParticleTraits<PType> Traits;
0076 LorentzMomentum p1 = Traits::momentum(pp.first);
0077 Boost b = (p1 + Traits::momentum(pp.second)).boostVector();
0078 p1.boost(-b);
0079 LorentzRotation rot;
0080 rot.rotateZ(-p1.phi());
0081 rot.rotateY(p1.theta());
0082 rot.rotateZ(p1.phi());
0083 rot.boost(b);
0084 return rot;
0085 }
0086
0087 template <typename Iterator>
0088 LorentzRotation UtilityBase::boostToCM(Iterator first, Iterator last,
0089 Iterator zAxis, Iterator xzPlane) {
0090 if ( first == last ) return LorentzRotation();
0091 LorentzRotation boost(-sumMomentum(first, last).boostVector());
0092 typedef typename std::iterator_traits<Iterator>::value_type PType;
0093 typedef ParticleTraits<PType> Traits;
0094
0095 transform(first, last, boost);
0096 if ( zAxis != last ) {
0097 LorentzRotation rot;
0098 rot.rotateZ(-Traits::momentum(*zAxis).phi());
0099 rot.rotateY(-Traits::momentum(*zAxis).theta());
0100 if ( xzPlane == last ) rot.rotateZ(Traits::momentum(*zAxis).phi());
0101 transform(first, last, rot);
0102 boost = rot*boost;
0103 }
0104 if ( xzPlane != last ) {
0105 LorentzRotation rot;
0106 rot.rotateZ(-Traits::momentum(*xzPlane).phi());
0107 transform(first, last, rot);
0108 boost = rot*boost;
0109 }
0110
0111 return boost;
0112 }
0113
0114 template <typename LV>
0115 LorentzRotation UtilityBase::
0116 transformFromCMS(const LV & p) {
0117 if(p.rho()/p.e()>=1.)
0118 throw Exception() << "Boost >=1. in "
0119 << "UtilityBase::transformFromCMS()"
0120 << Exception::eventerror;
0121 LorentzRotation r(0.0, 0.0, p.rho()/p.e());
0122 r.rotateZ(-p.phi());
0123 r.rotateY(p.theta());
0124 r.rotateZ(p.phi());
0125 return r;
0126 }
0127
0128 template <typename LV>
0129 LorentzRotation UtilityBase::
0130 transformFromCMS(const LV & sum, LV zAxis) {
0131 LorentzRotation r;
0132 Boost bz(sum.boostVector());
0133 zAxis.boost(-bz);
0134 r.rotateZ(-zAxis.phi());
0135 r.rotateY(zAxis.theta());
0136 r.rotateZ(zAxis.phi());
0137 r.boost(bz);
0138 return r;
0139 }
0140
0141 template <typename LV>
0142 LorentzRotation UtilityBase::
0143 transformFromCMS(const LV & sum, const LV & zAxis, LV xyPlane) {
0144 LorentzRotation r;
0145 LorentzRotation r1 = transformFromCMS(sum, zAxis);
0146 xyPlane.transform(r1.inverse());
0147 r.rotateZ(xyPlane.phi());
0148 r.transform(r1);
0149 return r;
0150 }
0151
0152 template <typename LV>
0153 LorentzRotation UtilityBase::
0154 transformToCMS(const LV & p) {
0155 LorentzRotation r;
0156 r.rotateZ(-p.phi());
0157 r.rotateY(-p.theta());
0158 r.rotateZ(p.phi());
0159 if(p.rho()/p.e()>1.) throw Exception() << "Boost >=1. in "
0160 << "UtilityBase::transformToCMS "
0161 << Exception::eventerror;
0162 r.boostZ(-p.rho()/p.e());
0163 return r;
0164 }
0165
0166 template <typename LV>
0167 LorentzRotation UtilityBase::
0168 transformToCMS(const LV & sum, LV zAxis) {
0169 LorentzRotation r(-sum.boostVector());
0170 zAxis.transform(r);
0171 r.rotateZ(-zAxis.phi());
0172 r.rotateY(-zAxis.theta());
0173 r.rotateZ(zAxis.phi());
0174 return r;
0175 }
0176
0177 template <typename LV>
0178 LorentzRotation UtilityBase::
0179 transformToCMS(const LV & sum, const LV & zAxis, LV xyPlane) {
0180 LorentzRotation r(transformToCMS(sum, zAxis));
0181 xyPlane.transform(r);
0182 r.rotateZ(-xyPlane.phi());
0183 return r;
0184 }
0185
0186 template <typename PType>
0187 void UtilityBase::setMomentum(PType & p, const Momentum3 & q) {
0188 typedef ParticleTraits<PType> Traits;
0189 if ( Traits::momentum(p).m2() <= ZERO ||
0190 Traits::mass(p) <= ZERO ) {
0191 Traits::set3Momentum(p, q);
0192 } else {
0193 Traits::transform(p, transformToCMS(Traits::momentum(p)));
0194 setMomentumFromCMS(p, q);
0195 }
0196 }
0197
0198 template <typename PType>
0199 void UtilityBase::setMomentumFromCMS(PType & p, const Momentum3 & q) {
0200 typedef ParticleTraits<PType> Traits;
0201 LorentzMomentum q4(q, sqrt(q.mag2() + sqr(Traits::mass(p))));
0202 if ( Traits::momentum(p) == LorentzMomentum() )
0203 Traits::set5Momentum(p, q4);
0204 else
0205 Traits::transform(p, transformFromCMS(q4));
0206 return;
0207 }
0208
0209 template <typename Iter>
0210 void UtilityBase::setMomentum(Iter first, Iter last, const Momentum3 & q) {
0211 if ( first == last ) return;
0212 typedef typename std::iterator_traits<Iter>::value_type PType;
0213 typedef ParticleTraits<PType> Traits;
0214 Iter second = first;
0215 if ( ++second == last ) return setMomentum(Traits::ref(*first), q);
0216 LorentzMomentum sum = sumMomentum(first, last);
0217 transform(first, last, transformToCMS(sum));
0218 setMomentumFromCMS(first, last, sum.m2(), q);
0219 }
0220
0221 template <typename Iter>
0222 void UtilityBase::
0223 setMomentum(Iter first, Iter last, const Momentum3 & q, double eps) {
0224 if ( first == last ) return;
0225 typedef typename std::iterator_traits<Iter>::value_type PType;
0226 typedef ParticleTraits<PType> Traits;
0227 Iter second = first;
0228 if ( ++second == last ) return setMomentum(Traits::ref(*first), q);
0229 LorentzRotation r;
0230 LorentzMomentum sum = sumMomentum(first, last);
0231 if ( sum.rho2() > 1.0e-12*MeV2 ) {
0232 r.rotateZ(-sum.phi());
0233 r.rotateY(-sum.theta());
0234 r.rotateZ(sum.phi());
0235 }
0236 if((sum.vect()-q).mag2()< 1.0e-12*MeV2) return;
0237 Energy2 ppo = sqr(sum.rho() + sum.e());
0238 Energy2 ppn = sqr(q.mag() + sqrt(q.mag2() + sum.m2()));
0239 r.boost(0.0, 0.0, (ppn - ppo)/(ppn + ppo));
0240 r.rotateZ(-q.phi());
0241 r.rotateY(q.theta());
0242 r.rotateZ(q.phi());
0243 transform(first, last, r);
0244 sum = sumMomentum(first, last);
0245 if ( abs(sum.x() - q.x()) > eps*sum.e() ||
0246 abs(sum.y() - q.y()) > eps*sum.e() ||
0247 abs(sum.z() - q.z()) > eps*sum.e() ) setMomentum(first, last, q, eps);
0248 }
0249
0250 template <typename Iter>
0251 void UtilityBase::
0252 setMomentumFromCMS(Iter first, Iter last, Energy2 m2, const Momentum3 & q) {
0253 if ( first == last ) return;
0254 typedef typename std::iterator_traits<Iter>::value_type PType;
0255 typedef ParticleTraits<PType> Traits;
0256 Iter second = first;
0257 if ( ++second == last ) return setMomentumFromCMS(Traits::ref(*first), q);
0258 transform(first, last,
0259 transformFromCMS(LorentzMomentum(q, sqrt(q.mag2() + m2))));
0260 }
0261
0262 template <typename Cont1, typename Cont2>
0263 void UtilityBase::add(Cont1 & c1, const Cont2 & c2) {
0264 typename Cont2::const_iterator i = c2.begin();
0265 while ( i != c2.end() ) c1.insert(c1.end(), *i++);
0266 }
0267
0268 }
