/usr/share/pythia8-examples/examples/main17.cc is in pythia8-examples 8.1.80-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 | // main17.cc is a part of the PYTHIA event generator.
// Copyright (C) 2013 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// This is a simple test program.
// It illustrates
// (a) how to use UserHooks to regularize onium cross section for pT -> 0,
// (b) how decays could be handled externally.
#include "Pythia8/Pythia.h"
using namespace Pythia8;
//==========================================================================
// A derived class to do J/psi decays.
class JpsiDecay : public DecayHandler {
public:
// Constructor.
JpsiDecay(ParticleData* pdtPtrIn, Rndm* rndmPtrIn) {times = 0;
pdtPtr = pdtPtrIn; rndmPtr = rndmPtrIn;}
// Routine for doing the decay.
bool decay(vector<int>& idProd, vector<double>& mProd,
vector<Vec4>& pProd, int iDec, const Event& event);
private:
// Count number of times JpsiDecay is called.
int times;
// Pointer to the particle data table.
ParticleData* pdtPtr;
// Pointer to the random number generator.
Rndm* rndmPtr;
};
//--------------------------------------------------------------------------
// The actual J/psi decay routine.
// Not intended for realism, just to illustrate the principles.
bool JpsiDecay::decay(vector<int>& idProd, vector<double>& mProd,
vector<Vec4>& pProd, int iDec, const Event& event) {
// Always do decay J/psi -> mu+ mu-; store the muons.
idProd.push_back(-13);
idProd.push_back(13);
// Muon mass(es), here from Pythia tables, also stored.
double mMuon = pdtPtr->m0(13);
mProd.push_back(mMuon);
mProd.push_back(mMuon);
// Calculate muon energy and momentum in J/psi rest frame.
double eMuon = 0.5 * mProd[0];
double pAbsMuon = sqrt(eMuon * eMuon - mMuon * mMuon);
// Assume decay angles isotropic in rest frame.
double cosTheta = 2. * rndmPtr->flat() - 1.;
double sinTheta = sqrt(max(0., 1. - cosTheta * cosTheta));
double phi = 2. * M_PI * rndmPtr->flat();
double pxMuon = pAbsMuon * sinTheta * cos(phi);
double pyMuon = pAbsMuon * sinTheta * sin(phi);
double pzMuon = pAbsMuon * cosTheta;
// Define mu+ and mu- four-vectors in the J/psi rest frame.
Vec4 pMuPlus( pxMuon, pyMuon, pzMuon, eMuon);
Vec4 pMuMinus( -pxMuon, -pyMuon, -pzMuon, eMuon);
// Boost them by velocity vector of the J/psi mother and store.
pMuPlus.bst(pProd[0]);
pMuMinus.bst(pProd[0]);
pProd.push_back(pMuPlus);
pProd.push_back(pMuMinus);
// Print message the first few times, to show that it works.
if (times++ < 10) {
int iMother = event[iDec].mother1();
int idMother = event[iMother].id();
cout << "\n J/psi decay performed, J/psi in line " << iDec
<< ", mother id = " << idMother << "\n";
}
// Done
return true;
}
//==========================================================================
int main() {
// Number of events to generate and to list. Max number of errors.
int nEvent = 2000;
int nList = 2;
int nAbort = 5;
// Pythia generator.
Pythia pythia;
// Initialization for charmonium (singlet+octet) production at the LHC.
pythia.readString("Charmonium:all = on");
pythia.readString("Beams:eCM = 7000.");
// Normally cutoff at pTHat = 1, but push it lower combined with dampening.
pythia.readString("PhaseSpace:pTHatMin = 0.5");
pythia.readString("PhaseSpace:pTHatMinDiverge = 0.5");
// Set up to do a user veto and send it in.
// First argument: multiplies the pT0 of multiparton interactions
// to define the pT dampeing scale.
// Second argument: howe many powers of alpha_strong to
// reweight with new (larger) argument.
// Third argument: choice of process scale two different ways;
// probably does not make much difference.
// See "User Hooks" in manual for detail on SuppressSmallPT.
UserHooks* oniumUserHook = new SuppressSmallPT( 1., 3, false);
pythia.setUserHooksPtr( oniumUserHook);
// A class to do J/psi decays externally.
DecayHandler* handleDecays = new JpsiDecay(&pythia.particleData,
&pythia.rndm);
// The list of particles the class can handle.
vector<int> handledParticles;
handledParticles.push_back(443);
// Hand pointer and list to Pythia.
pythia.setDecayPtr( handleDecays, handledParticles);
// Switch off automatic event listing in favour of manual.
pythia.readString("Next:numberShowInfo = 0");
pythia.readString("Next:numberShowProcess = 0");
pythia.readString("Next:numberShowEvent = 0");
// Initialization.
pythia.init();
// Book histograms.
Hist pThard("pTHat of hard subprocess", 100, 0., 50.);
Hist pTJPsi("pT of J/Psi", 100, 0., 50.);
// Begin event loop.
int iList = 0;
int iAbort = 0;
for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
// Generate events. Quit if many failures.
if (!pythia.next()) {
if (++iAbort < nAbort) continue;
cout << " Event generation aborted prematurely, owing to error!\n";
break;
}
// Histogram pThard spectrum of process.
double pTHat = pythia.info.pTHat();
pThard.fill( pTHat );
// Look for event with externally handled decays.
bool externalDecay = false;
for (int i = 0; i < pythia.event.size(); ++i) {
int status = pythia.event[i].statusAbs();
if (status == 93 || status == 94) {externalDecay = true; break;}
}
// List first few events with external decay.
if (externalDecay && ++iList <= nList) {
pythia.process.list();
pythia.event.list();
}
// Histogram pT spectrum of J/Psi.
for (int i = 0; i < pythia.event.size(); ++i)
if (pythia.event[i].id() == 443) pTJPsi.fill( pythia.event[i].pT() );
// End of event loop.
}
// Final statistics. Print histograms.
pythia.stat();
cout << pThard << pTJPsi;
// Done.
delete handleDecays;
delete oniumUserHook;
return 0;
}
|