[1] | 1 | // main22.cc is a part of the PYTHIA event generator. |
---|
| 2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
| 3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
| 4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
| 5 | |
---|
| 6 | // Simple illustration how to provide (a) your own resonance-width class, |
---|
| 7 | // and (b) your own cross-section class, with instances handed in to Pythia. |
---|
| 8 | // The hypothetical scenario is that top would have been so long-lived |
---|
| 9 | // that a toponium resonance Theta could form. Then production could |
---|
| 10 | // proceed via q qbar -> gamma*/Z* -> Theta, with decay either to |
---|
| 11 | // a fermion pair or (dominantly) to three gluons. |
---|
| 12 | // The implementation is not physically correct in any number of ways, |
---|
| 13 | // but should exemplify the strategy needed for realistic cases. |
---|
| 14 | |
---|
| 15 | #include "Pythia.h" |
---|
| 16 | |
---|
| 17 | using namespace Pythia8; |
---|
| 18 | |
---|
| 19 | //========================================================================== |
---|
| 20 | |
---|
| 21 | // The ResonanceTheta class handles a toponium resonance. |
---|
| 22 | |
---|
| 23 | class ResonanceTheta : public ResonanceWidths { |
---|
| 24 | |
---|
| 25 | public: |
---|
| 26 | |
---|
| 27 | // Constructor. |
---|
| 28 | ResonanceTheta(int idResIn) {initBasic(idResIn);} |
---|
| 29 | |
---|
| 30 | private: |
---|
| 31 | |
---|
| 32 | // Locally stored properties and couplings. |
---|
| 33 | double normTheta2qqbar, normTheta2llbar, normTheta2ggg; |
---|
| 34 | |
---|
| 35 | // Initialize constants. |
---|
| 36 | virtual void initConstants(); |
---|
| 37 | |
---|
| 38 | // Calculate various common prefactors for the current mass. |
---|
| 39 | // Superfluous here, so skipped. |
---|
| 40 | //virtual void calcPreFac(bool = false); |
---|
| 41 | |
---|
| 42 | // Calculate width for currently considered channel. |
---|
| 43 | virtual void calcWidth(bool = false); |
---|
| 44 | |
---|
| 45 | }; |
---|
| 46 | |
---|
| 47 | //-------------------------------------------------------------------------- |
---|
| 48 | |
---|
| 49 | // Initialize constants. |
---|
| 50 | |
---|
| 51 | void ResonanceTheta::initConstants() { |
---|
| 52 | |
---|
| 53 | // Dummy normalization of couplings to the allowed decay channels. |
---|
| 54 | normTheta2qqbar = 0.0001; |
---|
| 55 | normTheta2llbar = 0.0001; |
---|
| 56 | normTheta2ggg = 0.001; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | //-------------------------------------------------------------------------- |
---|
| 60 | |
---|
| 61 | // Calculate width for currently considered channel. |
---|
| 62 | |
---|
| 63 | void ResonanceTheta::calcWidth(bool) { |
---|
| 64 | |
---|
| 65 | // Expression for Theta -> q qbar (q up to b). Colour factor. |
---|
| 66 | if (id1Abs < 6) widNow = 3. * normTheta2qqbar * mHat; |
---|
| 67 | |
---|
| 68 | // Expression for Theta -> l lbar (l = e, mu, tau). |
---|
| 69 | else if (id1Abs == 11 || id1Abs == 13 || id1Abs == 15) |
---|
| 70 | widNow = normTheta2llbar * mHat; |
---|
| 71 | |
---|
| 72 | // Expression for Theta -> g g g. Colour factor. |
---|
| 73 | else if (id1Abs == 21) widNow = 8. * normTheta2ggg * mHat; |
---|
| 74 | |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | //========================================================================== |
---|
| 78 | |
---|
| 79 | // A derived class for q qbar -> Theta (toponium bound state). |
---|
| 80 | |
---|
| 81 | class Sigma1qqbar2Theta : public Sigma1Process { |
---|
| 82 | |
---|
| 83 | public: |
---|
| 84 | |
---|
| 85 | // Constructor. |
---|
| 86 | Sigma1qqbar2Theta() {} |
---|
| 87 | |
---|
| 88 | // Initialize process. |
---|
| 89 | virtual void initProc(); |
---|
| 90 | |
---|
| 91 | // Calculate flavour-independent parts of cross section. |
---|
| 92 | virtual void sigmaKin(); |
---|
| 93 | |
---|
| 94 | // Evaluate sigmaHat(sHat). Assumed flavour-independent so simple. |
---|
| 95 | virtual double sigmaHat() {return sigma;} |
---|
| 96 | |
---|
| 97 | // Select flavour, colour and anticolour. |
---|
| 98 | virtual void setIdColAcol(); |
---|
| 99 | |
---|
| 100 | // Evaluate weight for decay angles. |
---|
| 101 | virtual double weightDecay( Event& process, int iResBeg, int iResEnd); |
---|
| 102 | |
---|
| 103 | // Info on the subprocess. |
---|
| 104 | virtual string name() const {return "q qbar -> Theta";} |
---|
| 105 | virtual int code() const {return 621;} |
---|
| 106 | virtual string inFlux() const {return "qqbarSame";} |
---|
| 107 | virtual int resonanceA() const {return 663;} |
---|
| 108 | |
---|
| 109 | private: |
---|
| 110 | |
---|
| 111 | // Store flavour-specific process information and standard prefactor. |
---|
| 112 | int idTheta; |
---|
| 113 | double mRes, GammaRes, m2Res, GamMRat, normTheta2qqbar, sigma; |
---|
| 114 | |
---|
| 115 | // Pointer to properties of Theta, to access decay width. |
---|
| 116 | ParticleDataEntry* particlePtr; |
---|
| 117 | |
---|
| 118 | }; |
---|
| 119 | |
---|
| 120 | //-------------------------------------------------------------------------- |
---|
| 121 | |
---|
| 122 | // Initialize process. |
---|
| 123 | |
---|
| 124 | void Sigma1qqbar2Theta::initProc() { |
---|
| 125 | |
---|
| 126 | // Store Theta mass and width for propagator. |
---|
| 127 | idTheta = 663; |
---|
| 128 | mRes = particleDataPtr->m0(idTheta); |
---|
| 129 | GammaRes = particleDataPtr->mWidth(idTheta); |
---|
| 130 | m2Res = mRes*mRes; |
---|
| 131 | GamMRat = GammaRes / mRes; |
---|
| 132 | |
---|
| 133 | // Same normlization as in ResonanceTheta for coupling strength. |
---|
| 134 | normTheta2qqbar = 0.0001; |
---|
| 135 | |
---|
| 136 | // Set pointer to particle properties and decay table. |
---|
| 137 | particlePtr = particleDataPtr->particleDataEntryPtr(idTheta); |
---|
| 138 | |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | //-------------------------------------------------------------------------- |
---|
| 142 | |
---|
| 143 | // Evaluate sigmaHat(sHat); first step when inflavours unknown. |
---|
| 144 | |
---|
| 145 | void Sigma1qqbar2Theta::sigmaKin() { |
---|
| 146 | |
---|
| 147 | // Incoming width with colour factor. |
---|
| 148 | double widthIn = normTheta2qqbar * mH / 3.; |
---|
| 149 | |
---|
| 150 | // Breit-Wigner, including some (guessed) spin factors. |
---|
| 151 | double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
| 152 | |
---|
| 153 | // Outgoing width: only includes channels left open. |
---|
| 154 | double widthOut = particlePtr->resWidthOpen(663, mH); |
---|
| 155 | |
---|
| 156 | // Total answer. |
---|
| 157 | sigma = widthIn * sigBW * widthOut; |
---|
| 158 | |
---|
| 159 | } |
---|
| 160 | |
---|
| 161 | //-------------------------------------------------------------------------- |
---|
| 162 | |
---|
| 163 | // Select identity, colour and anticolour. |
---|
| 164 | |
---|
| 165 | void Sigma1qqbar2Theta::setIdColAcol() { |
---|
| 166 | |
---|
| 167 | // Flavours trivial. |
---|
| 168 | setId( id1, id2, idTheta); |
---|
| 169 | |
---|
| 170 | // Colour flow topologies. Swap when antiquarks. |
---|
| 171 | setColAcol( 1, 0, 0, 1, 0, 0); |
---|
| 172 | if (id1 < 0) swapColAcol(); |
---|
| 173 | |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | //-------------------------------------------------------------------------- |
---|
| 177 | |
---|
| 178 | // Evaluate weight for Theta -> g g g. |
---|
| 179 | |
---|
| 180 | double Sigma1qqbar2Theta::weightDecay( Event& process, int iResBeg, |
---|
| 181 | int iResEnd) { |
---|
| 182 | |
---|
| 183 | // Should be Theta decay. (This is only option here, so overkill.) |
---|
| 184 | if (iResEnd != iResBeg || process[iResBeg].idAbs() != idTheta) |
---|
| 185 | return 1.; |
---|
| 186 | |
---|
| 187 | // Should be decay to three gluons. |
---|
| 188 | int i1 = process[iResBeg].daughter1(); |
---|
| 189 | int i2 = i1 + 1; |
---|
| 190 | int i3 = i2 + 1; |
---|
| 191 | if (i3 != process[iResBeg].daughter2() || process[i1].id() != 21) |
---|
| 192 | return 1.; |
---|
| 193 | |
---|
| 194 | // Energy fractions x_i = 2 E_i/m_Theta of gluons in Theta rest frame. |
---|
| 195 | double x1 = 2. * process[i1].p() * process[iResBeg].p() |
---|
| 196 | / process[iResBeg].m2(); |
---|
| 197 | double x2 = 2. * process[i2].p() * process[iResBeg].p() |
---|
| 198 | / process[iResBeg].m2(); |
---|
| 199 | double x3 = 2. * process[i3].p() * process[iResBeg].p() |
---|
| 200 | / process[iResBeg].m2(); |
---|
| 201 | |
---|
| 202 | // Matrix-element expression for Theta -> g g g. |
---|
| 203 | double wtME = pow2( (1. - x1) / (x2 * x3) ) |
---|
| 204 | + pow2( (1. - x2) / (x1 * x3) ) + pow2( (1. - x3) / (x1 * x2) ); |
---|
| 205 | double wtMEmax = 2.; |
---|
| 206 | return wtME / wtMEmax; |
---|
| 207 | |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | //========================================================================== |
---|
| 211 | |
---|
| 212 | int main() { |
---|
| 213 | |
---|
| 214 | // Number of events to generate. Max number of errors. |
---|
| 215 | // Warning: generation of complete events is much slower than if you use |
---|
| 216 | // PartonLevel:all = off to only get cross sections, so adjust nEvent. |
---|
| 217 | int nEvent = 1000; |
---|
| 218 | int nAbort = 5; |
---|
| 219 | |
---|
| 220 | // Pythia generator. |
---|
| 221 | Pythia pythia; |
---|
| 222 | |
---|
| 223 | // Create the toponium resonance and a few production/decay channels. |
---|
| 224 | // Warning: many more exist, e.g. weak ones of one top quark. |
---|
| 225 | // Note: to obtain the correct width for the Breit-Wigner you must |
---|
| 226 | // include all channels, but you only need leave those on that you |
---|
| 227 | // want to study. |
---|
| 228 | pythia.readString("663:new = Theta void 3 0 0 342.0 0.2 300. 400. 0."); |
---|
| 229 | pythia.readString("663:addChannel = 1 0. 0 1 -1"); |
---|
| 230 | pythia.readString("663:addChannel = 1 0. 0 2 -2"); |
---|
| 231 | pythia.readString("663:addChannel = 1 0. 0 3 -3"); |
---|
| 232 | pythia.readString("663:addChannel = 1 0. 0 4 -4"); |
---|
| 233 | pythia.readString("663:addChannel = 1 0. 0 5 -5"); |
---|
| 234 | pythia.readString("663:addChannel = 1 0. 0 11 -11"); |
---|
| 235 | pythia.readString("663:addChannel = 1 0. 0 13 -13"); |
---|
| 236 | pythia.readString("663:addChannel = 1 0. 0 15 -15"); |
---|
| 237 | pythia.readString("663:addChannel = 1 0. 0 21 21 21"); |
---|
| 238 | |
---|
| 239 | // Create instance of a class to calculate the width of Theta to the |
---|
| 240 | // above channels. Hand in pointer to Pythia. |
---|
| 241 | ResonanceWidths* resonanceTheta = new ResonanceTheta(663); |
---|
| 242 | pythia.setResonancePtr(resonanceTheta); |
---|
| 243 | |
---|
| 244 | // Create instance of a class to generate the q qbar -> Theta process |
---|
| 245 | // from an external matrix element. Hand in pointer to Pythia. |
---|
| 246 | SigmaProcess* sigma1Theta = new Sigma1qqbar2Theta(); |
---|
| 247 | pythia.setSigmaPtr(sigma1Theta); |
---|
| 248 | |
---|
| 249 | // Optionally only compare cross sections. |
---|
| 250 | //pythia.readString("PartonLevel:all = off"); |
---|
| 251 | pythia.readString("Check:nErrList = 2"); |
---|
| 252 | |
---|
| 253 | // Initialization for LHC. |
---|
| 254 | pythia.init(); |
---|
| 255 | |
---|
| 256 | // Book histogram. |
---|
| 257 | Hist mTheta("Theta mass", 100, 300., 400.); |
---|
| 258 | |
---|
| 259 | // Begin event loop. |
---|
| 260 | int iAbort = 0; |
---|
| 261 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
| 262 | |
---|
| 263 | // Generate events. Quit if many failures. |
---|
| 264 | if (!pythia.next()) { |
---|
| 265 | if (++iAbort < nAbort) continue; |
---|
| 266 | cout << " Event generation aborted prematurely, owing to error!\n"; |
---|
| 267 | break; |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | // Fill Theta mass. End of event loop. |
---|
| 271 | mTheta.fill( pythia.process[5].m() ); |
---|
| 272 | } |
---|
| 273 | |
---|
| 274 | // Final statistics. Print histogram. |
---|
| 275 | pythia.stat(); |
---|
| 276 | cout << mTheta; |
---|
| 277 | |
---|
| 278 | // Done. |
---|
| 279 | return 0; |
---|
| 280 | } |
---|