PYTHIA8/pythia8145/examples/main10.cc

   1 // main10.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 // Example how you can use UserHooks to trace pT spectrum through program,
   7 // and veto undesirable jet multiplicities.
   8
   9 #include "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 cone-jet finder with (etaMax, nEta, nPhi).
  35   MyUserHooks() { coneJet = new CellJet(5., 100, 64); }
  36
  37   // Destructor deletes cone-jet finder.
  38   ~MyUserHooks() {delete coneJet;}
  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 eTjet > 10 for R = 0.7.
  82     coneJet->analyze(event, 10., 0.7);
  83     int nJet = coneJet->size();
  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 cone-jet finder.
 114   CellJet* coneJet;
 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.init( 2212, -2212, 1960.);
 139
 140   // Begin event loop.
 141   for (int iEvent = 0; iEvent < 1000; ++iEvent) {
 142
 143     // Generate events.
 144     pythia.next();
 145
 146   // End of event loop.
 147   }
 148
 149   // Statistics. Histograms.
 150   pythia.statistics();
 151   cout << pTtrial << pTselect << pTaccept
 152        << nPartonsB << nJets << nPartonsA
 153        << nFSRatISR;
 154
 155   // Done.
 156   return 0;
 157 }