StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main10.cc
1 // main10.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 // Example how you can use UserHooks to trace pT spectrum through program,
7 // and veto undesirable jet multiplicities.
8 
9 #include "Pythia8/Pythia.h"
10 using namespace Pythia8;
11 
12 //==========================================================================
13 
14 // Put histograms here to make them global, so they can be used both
15 // in MyUserHooks and in the main program.
16 
17 Hist pTtrial("trial pT spectrum", 100, 0., 400.);
18 Hist pTselect("selected pT spectrum (before veto)", 100, 0., 400.);
19 Hist pTaccept("accepted pT spectrum (after veto)", 100, 0., 400.);
20 Hist nPartonsB("number of partons before veto", 20, -0.5, 19.5);
21 Hist nJets("number of jets before veto", 20, -0.5, 19.5);
22 Hist nPartonsA("number of partons after veto", 20, -0.5, 19.5);
23 Hist nFSRatISR("number of FSR emissions at first ISR emission",
24  20, -0.5, 19.5);
25 
26 //==========================================================================
27 
28 // Write own derived UserHooks class.
29 
30 class MyUserHooks : public UserHooks {
31 
32 public:
33 
34  // Constructor creates anti-kT jet finder with (-1, R, pTmin, etaMax).
35  MyUserHooks() { slowJet = new SlowJet(-1, 0.7, 10., 5.); }
36 
37  // Destructor deletes anti-kT jet finder.
38  ~MyUserHooks() {delete slowJet;}
39 
40  // Allow process cross section to be modified...
41  virtual bool canModifySigma() {return true;}
42 
43  // ...which gives access to the event at the trial level, before selection.
44  virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr,
45  const PhaseSpace* phaseSpacePtr, bool inEvent) {
46 
47  // All events should be 2 -> 2, but kill them if not.
48  if (sigmaProcessPtr->nFinal() != 2) return 0.;
49 
50  // Extract the pT for 2 -> 2 processes in the event generation chain
51  // (inEvent = false for initialization).
52  if (inEvent) {
53  pTHat = phaseSpacePtr->pTHat();
54  // Fill histogram of pT spectrum.
55  pTtrial.fill( pTHat );
56  }
57 
58  // Here we do not modify 2 -> 2 cross sections.
59  return 1.;
60  }
61 
62  // Allow a veto for the interleaved evolution in pT.
63  virtual bool canVetoPT() {return true;}
64 
65  // Do the veto test at a pT scale of 5 GeV.
66  virtual double scaleVetoPT() {return 5.;}
67 
68  // Access the event in the interleaved evolution.
69  virtual bool doVetoPT(int iPos, const Event& event) {
70 
71  // iPos <= 3 for interleaved evolution; skip others.
72  if (iPos > 3) return false;
73 
74  // Fill histogram of pT spectrum at this stage.
75  pTselect.fill(pTHat);
76 
77  // Extract a copy of the partons in the hardest system.
78  subEvent(event);
79  nPartonsB.fill( workEvent.size() );
80 
81  // Find number of jets with given conditions.
82  slowJet->analyze(event);
83  int nJet = slowJet->sizeJet();
84  nJets.fill( nJet );
85 
86  // Veto events which do not have exactly three jets.
87  if (nJet != 3) return true;
88 
89  // Statistics of survivors.
90  nPartonsA.fill( workEvent.size() );
91  pTaccept.fill(pTHat);
92 
93  // Do not veto events that got this far.
94  return false;
95 
96  }
97 
98  // Allow a veto after (by default) first step.
99  virtual bool canVetoStep() {return true;}
100 
101  // Access the event in the interleaved evolution after first step.
102  virtual bool doVetoStep( int iPos, int nISR, int nFSR, const Event& ) {
103 
104  // Only want to study what happens at first ISR emission
105  if (iPos == 2 && nISR == 1) nFSRatISR.fill( nFSR );
106 
107  // Not intending to veto any events here.
108  return false;
109  }
110 
111 private:
112 
113  // The anti-kT (or kT, C/A) jet finder.
114  SlowJet* slowJet;
115 
116  // Save the pThat scale.
117  double pTHat;
118 
119 };
120 
121 //==========================================================================
122 
123 int main() {
124 
125  // Generator.
126  Pythia pythia;
127 
128  // Process selection. No need to study hadron level.
129  pythia.readString("HardQCD:all = on");
130  pythia.readString("PhaseSpace:pTHatMin = 50.");
131  pythia.readString("HadronLevel:all = off");
132 
133  // Set up to do a user veto and send it in.
134  MyUserHooks* myUserHooks = new MyUserHooks();
135  pythia.setUserHooksPtr( myUserHooks);
136 
137  // Tevatron initialization.
138  pythia.readString("Beams:idB = -2212");
139  pythia.readString("Beams:eCM = 1960.");
140  pythia.init();
141 
142  // Begin event loop.
143  for (int iEvent = 0; iEvent < 1000; ++iEvent) {
144 
145  // Generate events.
146  pythia.next();
147 
148  // End of event loop.
149  }
150 
151  // Statistics. Histograms.
152  pythia.stat();
153  cout << pTtrial << pTselect << pTaccept
154  << nPartonsB << nJets << nPartonsA
155  << nFSRatISR;
156 
157  // Done.
158  delete myUserHooks;
159  return 0;
160 }