PYTHIA8/pythia8145/examples/main15.cc

   1 // main15.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 // This is a simple test program.
   7 // It illustrates how either
   8 // (a) B decays (sections marked "Repeated decays:), or
   9 // (b) all hadronization (sections marked "Repeated hadronization:")
  10 // could be repeated a number of times for each event,
  11 // to improve statistics when this could be a problem.
  12 // Option (a) is faster than (b), but less generic.
  13
  14 #include "Pythia.h"
  15 using namespace Pythia8;
  16
  17 int main() {
  18
  19   // Main switchs: redo B decays only or redo all hadronization, but not both.
  20   bool redoBDecays = false;
  21   bool redoHadrons = true;
  22   if (redoHadrons) redoBDecays = false;
  23
  24   // Number of events, generated and listed ones.
  25   int nEvent = 100;
  26   int nList = 1;
  27
  28   // Number of times decays/hadronization should be redone for each event.
  29   int nRepeat = 10;
  30   if (!redoBDecays && !redoHadrons) nRepeat = 1;
  31
  32   // Generator. Shorthand for event.
  33   Pythia pythia;
  34   Event& event = pythia.event;
  35
  36   // Simulate b production above given pTmin scale.
  37   // Warning: these processes do not catch all possible production modes.
  38   // You would need to use HardQCD:all or even SoftQCD:minBias for that.
  39   pythia.readString("HardQCD:gg2bbbar = on");
  40   pythia.readString("HardQCD:qqbar2bbbar = on");
  41   pythia.readString("PhaseSpace:pTHatMin = 50.");
  42
  43   // Repeated decays: list of weakly decaying B hadrons.
  44   // Note: this list is overkill; some will never be produced.
  45   int bCodes[28] = {511, 521, 531, 541, 5122, 5132, 5142, 5232, 5242,
  46     5332, 5342, 5412, 5414, 5422, 5424, 5432, 5434, 5442, 5444, 5512,
  47     5514, 5522, 5524, 5532, 5534, 5542, 5544, 5544 };
  48   int nCodes = 28;
  49
  50   // Repeated decays: location of B handrons.
  51   vector<int> iBHad;
  52   int nBHad = 0;
  53
  54   // Repeated hadronization: spare copy of event.
  55   Event savedEvent;
  56
  57   // Repeated hadronization: switch off normal HadronLevel call.
  58   if (redoHadrons) pythia.readString("HadronLevel:all = off");
  59
  60   // Initialize for LHC energies.
  61   pythia.init( 2212, 2212, 14000.);
  62
  63   // Show changed settings and particleData.
  64   pythia.settings.listChanged();
  65   pythia.particleData.listChanged();
  66
  67   // Histogram invariant mass of muon pairs.
  68   Hist nBperEvent("number of b quarks in an event", 10, -0.5, 9.5);
  69   Hist nSameEvent("number of times same event is used", 10, -0.5, 9.5);
  70   Hist oppSignMass("mass of opposite-sign muon pair", 100, 0.0, 100.0);
  71   Hist sameSignMass("mass of same-sign muon pair", 100, 0.0, 100.0);
  72
  73   // Begin event loop.
  74   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
  75
  76     // Repeated decays: switch off decays of weakly decaying B hadrons.
  77     // (More compact solution than repeated readString(..).)
  78     if (redoBDecays) for (int iC = 0; iC < nCodes; ++iC)
  79       pythia.particleData.mayDecay( bCodes[iC], false);
  80
  81     // Generate event. Skip it if error.
  82     if (!pythia.next()) continue;
  83
  84     // List first few events.
  85     if (iEvent < nList) {
  86       pythia.info.list();
  87       pythia.process.list();
  88       event.list();
  89     }
  90
  91     // Find and histogram number of b quarks.
  92     int nBquark = 0;
  93     int stat;
  94     for (int i = 0; i < event.size(); ++i) {
  95       stat = event[i].statusAbs();
  96       if (event[i].idAbs() == 5 && (stat == 62 || stat == 63)) ++nBquark;
  97     }
  98     nBperEvent.fill( nBquark );
  99
 100     // Repeated decays: find all locations where B hadrons are stored.
 101     if (redoBDecays) {
 102       iBHad.resize(0);
 103       for (int i = 0; i < event.size(); ++i) {
 104         int idAbs = event[i].idAbs();
 105         for (int iC = 0; iC < 28; ++iC)
 106         if (idAbs == bCodes[iC]) {
 107           iBHad.push_back(i);
 108           break;
 109         }
 110       }
 111
 112       // Repeated decays: check that #b = #B.
 113       nBHad = iBHad.size();
 114       if (nBquark != nBHad) cout << " Warning: " << nBquark
 115         << " b quarks but " << nBHad << " B hadrons" << endl;
 116
 117       // Repeated decays: store size of current event.
 118       event.saveSize();
 119
 120       // Repeated decays: switch back on weakly decaying B hadrons.
 121       for (int iC = 0; iC < nCodes; ++iC)
 122         pythia.particleData.mayDecay( bCodes[iC], true);
 123
 124     //  Repeated hadronization: copy event into spare position.
 125     } else if (redoHadrons) {
 126       savedEvent = event;
 127     }
 128
 129     // Begin loop over rounds of decays / hadronization for same event.
 130     int nWithPair = 0;
 131     for (int iRepeat = 0; iRepeat < nRepeat; ++iRepeat) {
 132
 133       // Repeated decays: remove B decay products from previous round.
 134       if (redoBDecays) {
 135         if (iRepeat > 0) {
 136           event.restoreSize();
 137
 138           // Repeated decays: mark decayed B hadrons as undecayed.
 139           for (int iB = 0; iB < nBHad; ++iB) event[ iBHad[iB] ].statusPos();
 140         }
 141
 142         // Repeated decays: do decays of B hadrons, sequentially for products.
 143         // Note: modeDecays does not work for bottomonium (or heavier) states,
 144         // since there decays like Upsilon -> g g g also need hadronization.
 145         // Also, there is no provision for Bose-Einstein effects.
 146         if (!pythia.moreDecays()) continue;
 147
 148
 149       // Repeated hadronization: restore saved event record.
 150       } else if (redoHadrons) {
 151         if (iRepeat > 0) event = savedEvent;
 152
 153         // Repeated hadronization: do HadronLevel (repeatedly).
 154         if (!pythia.forceHadronLevel()) continue;
 155       }
 156
 157       // List last repetition of first few events.
 158       if ( (redoBDecays || redoHadrons) && iEvent < nList
 159         && iRepeat == nRepeat - 1) event.list();
 160
 161       // Look for muons among decay products (also from charm/tau/...).
 162       vector<int> iMuNeg, iMuPos;
 163       for (int i = 0; i < event.size(); ++i) {
 164         int id = event[i].id();
 165         if (id ==  13) iMuNeg.push_back(i);
 166         if (id == -13) iMuPos.push_back(i);
 167       }
 168
 169       // Check whether pair(s) present.
 170       int nMuNeg = iMuNeg.size();
 171       int nMuPos = iMuPos.size();
 172       if (nMuNeg + nMuPos > 1) {
 173         ++nWithPair;
 174
 175         // Fill masses of opposite-sign pairs.
 176         for (int iN = 0; iN < nMuNeg; ++iN)
 177         for (int iP = 0; iP < nMuPos; ++iP)
 178           oppSignMass.fill(
 179             (event[iMuNeg[iN]].p() + event[iMuPos[iP]].p()).mCalc() );
 180
 181         // Fill masses of same-sign pairs.
 182         for (int i1 = 0; i1 < nMuNeg - 1; ++i1)
 183         for (int i2 = i1 + 1; i2 < nMuNeg; ++i2)
 184           sameSignMass.fill(
 185             (event[iMuNeg[i1]].p() + event[iMuNeg[i2]].p()).mCalc() );
 186         for (int i1 = 0; i1 < nMuPos - 1; ++i1)
 187         for (int i2 = i1 + 1; i2 < nMuPos; ++i2)
 188           sameSignMass.fill(
 189             (event[iMuPos[i1]].p() + event[iMuPos[i2]].p()).mCalc() );
 190
 191       // Finished analysis of current round.
 192       }
 193
 194     // End of loop over many rounds. fill number of rounds with pairs.
 195     }
 196     nSameEvent.fill( nWithPair );
 197
 198   // End of event loop.
 199   }
 200
 201   // Statistics. Histograms.
 202   pythia.statistics();
 203   cout << nBperEvent << nSameEvent << oppSignMass << sameSignMass << endl;
 204
 205   // Done.
 206   return 0;
 207 }