PYTHIA8/pythia8145/examples/main25.cc

   1 // main25.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 your own cross-section class,
   7 // with an instance handed in to Pythia for internal generation.
   8 // The class is a lighly modified copy of an existing Pythia class,
   9 // and cross sections are compared with the internal implementation.
  10
  11 #include "Pythia.h"
  12
  13 using namespace Pythia8;
  14
  15 //==========================================================================
  16
  17 // A derived class for b g -> H0 b (Standard Model Higgs).
  18
  19 class Sigma2bg2Hb : public Sigma2Process {
  20
  21 public:
  22
  23   // Constructor.
  24   Sigma2bg2Hb() {}
  25
  26   // Initialize process.
  27   virtual void initProc();
  28
  29   // Calculate flavour-independent parts of cross section.
  30   virtual void sigmaKin();
  31
  32   // Evaluate sigmaHat(sHat).
  33   virtual double sigmaHat();
  34
  35   // Select flavour, colour and anticolour.
  36   virtual void setIdColAcol();
  37
  38   // Evaluate weight for decay angles.
  39   virtual double weightDecay( Event& process, int iResBeg, int iResEnd);
  40
  41   // Info on the subprocess.
  42   virtual string name()    const {return "b g -> H b";}
  43   virtual int    code()    const {return 10001;}
  44   virtual string inFlux()  const {return "qg";}
  45   virtual int    id3Mass() const {return 25;}
  46   virtual int    id4Mass() const {return 5;}
  47
  48 private:
  49
  50   // Store flavour-specific process information and standard prefactor.
  51   double m2b, m2W, thetaWRat, openFrac, sigma;
  52
  53 };
  54
  55 //--------------------------------------------------------------------------
  56
  57 // Initialize process.
  58
  59 void Sigma2bg2Hb::initProc() {
  60
  61   // Masses and couplings.
  62   m2b       = pow2( particleDataPtr->m0(5) );
  63   m2W       = pow2( particleDataPtr->m0(24) );
  64   thetaWRat = 1. / (24. * couplingsPtr->sin2thetaW());
  65
  66   // Suppression from secondary widths.
  67   openFrac = particleDataPtr->resOpenFrac(25);
  68
  69 }
  70
  71 //--------------------------------------------------------------------------
  72
  73 // Evaluate sigmaHat(sHat); first step when inflavours unknown.
  74
  75 void Sigma2bg2Hb::sigmaKin() {
  76
  77   // Initial values. Running mass provides coupling.
  78   double mHat   = sqrt(sH);
  79   double m2bRun = pow2( particleDataPtr->mRun(5, mHat) );
  80
  81   // Cross section.
  82   sigma = (M_PI / sH2) * alpS * alpEM * thetaWRat * (m2bRun/m2W)
  83     * ( sH / (m2b - uH) + 2. * m2b * (s3 - uH) / pow2(m2b - uH)
  84     + (m2b - uH) / sH - 2. * m2b / (m2b - uH)
  85     + 2. * (s3 - uH)  * (s3 - m2b - sH) / ((m2b - uH) * sH) );
  86   sigma *= openFrac;
  87
  88 }
  89
  90 //--------------------------------------------------------------------------
  91
  92 // Evaluate sigmaHat(sHat); second step for given inflavours.
  93
  94 double Sigma2bg2Hb::sigmaHat() {
  95
  96   // Cross section vanishing except for incoming b(bar) quark.
  97   if (abs(id1) != 5 && abs(id2) != 5) return 0.;
  98
  99   // Cross section as already calculated.
 100   return sigma;
 101
 102 }
 103
 104 //--------------------------------------------------------------------------
 105
 106 // Select identity, colour and anticolour.
 107
 108 void Sigma2bg2Hb::setIdColAcol() {
 109
 110   // Flavour set up for b g -> H0 q.
 111   int idq = (id2 == 21) ? id1 : id2;
 112   setId( id1, id2, 25, idq);
 113
 114   // tH defined between b_in and b_out: must swap tHat <-> uHat if b g in.
 115   swapTU = (id2 == 21);
 116
 117   // Colour flow topologies. Swap when antiquarks.
 118   if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
 119   else           setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
 120   if (idq < 0) swapColAcol();
 121
 122 }
 123
 124 //--------------------------------------------------------------------------
 125
 126 // Evaluate weight for Z0 Z0 or W+W- decay angles in Higgs decay.
 127
 128 double Sigma2bg2Hb::weightDecay( Event& process, int iResBeg,
 129   int iResEnd) {
 130
 131   // For Higgs decay hand over to standard routine, else done.
 132   if (process[process[iResBeg].mother1()].idAbs() == 25)
 133        return weightHiggsDecay( process, iResBeg, iResEnd);
 134   else return 1.;
 135
 136 }
 137
 138 //==========================================================================
 139
 140 int main() {
 141
 142   // Number of events to generate and to list. Max number of errors.
 143   // Warning: generation of complete events is much slower than if you use
 144   // PartonLevel:all = off to only get cross sections, so adjust nEvent.
 145   int nEvent = 1000;
 146   int nList  = 1;
 147   int nAbort = 5;
 148
 149   // Pythia generator.
 150   Pythia pythia;
 151
 152   // A class to generate the b g -> H b process from external matrix element.
 153   SigmaProcess* sigma2bg2Hb = new Sigma2bg2Hb();
 154
 155   // Hand pointer to Pythia.
 156   pythia.setSigmaPtr( sigma2bg2Hb);
 157
 158   // Optionally compare with (almost) same process implemented internally.
 159   pythia.readString("HiggsSM:qg2Hq = on");
 160
 161   // Phase space cut on pThat.
 162   pythia.readString("PhaseSpace:pTHatMin = 20.");
 163
 164   // Optionally only want to study H0 -> gamma gamma channel.
 165   pythia.readString("25:onMode = off");
 166   pythia.readString("25:onIfMatch = 22 22");
 167
 168   // Optionally only compare cross sections.
 169   //pythia.readString("PartonLevel:all = off");
 170
 171   // Initialization for LHC.
 172   pythia.init( 2212, 2212, 14000.);
 173
 174   // List changes.
 175   pythia.settings.listChanged();
 176   pythia.particleData.listChanged();
 177
 178   // Begin event loop.
 179   int iAbort = 0;
 180   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
 181     if (iEvent%(max(1,nEvent/20)) == 0) cout << " Now begin event "
 182       << iEvent << "\n";
 183
 184     // Generate events. Quit if many failures.
 185     if (!pythia.next()) {
 186       if (++iAbort < nAbort) continue;
 187       cout << " Event generation aborted prematurely, owing to error!\n";
 188       break;
 189     }
 190
 191     // List first few events.
 192     if (iEvent < nList) {
 193       pythia.info.list();
 194       pythia.process.list();
 195       pythia.event.list();
 196     }
 197
 198   // End of event loop.
 199   }
 200
 201   // Final statistics.
 202   pythia.statistics();
 203
 204   // Done.
 205   return 0;
 206 }