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