PYTHIA8/pythia8145/examples/main51.cc

   1 // main51.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 illustrating how to link in Pythia 6.4
   7 // for the generation of hard processes. That part is now considered
   8 // obsolete, but for some debug work this example has been collected
   9 // from various code pieces that were separately available up until
  10 // Pythia 8.125. In addition to modifying the below code to fit your
  11 // needs, you also have to modify examples/Makefile to link properly
  12 // for main51, including the Pythia 6.4xx library to be used (where
  13 // xx is the current subversion number).
  14
  15 // All hard PYTHIA 6.4 processes should be available for full generation
  16 // in PYTHIA 8, at least to the extent that they are defined for p p,
  17 // pbar p or e+ e-collisions. Soft processes, i.e. elastic and diffractive
  18 // scattering, as well as minimum-bias events, require a different
  19 // kinematics machinery, and can only be generated with the internal
  20 // PYTHIA 8 processes.
  21
  22 // PYTHIA 6.4 does its own rejection of events internally, according to
  23 // the strategy option 3. However, the current runtime interface does not
  24 // take cross-section information from PYTHIA 6.4. This means that both
  25 // the initial maxima and the final cross sections printed by the PYTHIA 8
  26 // routines are irrelevant in this case. Instead you have to study the
  27 // standard PYTHIA 6.4 initialization printout and call on pystat(...)
  28 // at the end of the run to obtain this information. It also means that
  29 // you cannot mix with internally generated PYTHIA 8 events.
  30
  31 //==========================================================================
  32
  33 #include "Pythia.h"
  34 #include "LHAFortran.h"
  35
  36 using namespace Pythia8;
  37
  38 //==========================================================================
  39
  40 // Declare the Fortran subroutines that may be used.
  41 // This code section is generic.
  42
  43 #ifdef _WIN32
  44   #define pygive_ PYGIVE
  45   #define pyinit_ PYINIT
  46   #define pyupin_ PYUPIN
  47   #define pyupev_ PYUPEV
  48   #define pylist_ PYLIST
  49   #define pystat_ PYSTAT
  50 #endif
  51
  52 extern "C" {
  53 #ifdef _WIN32
  54   extern void pyinit_(const char*, int, const char*, int, const char*, int, double&);
  55 #else
  56   extern void pyinit_(const char*, const char*, const char*, double&, int, int, int);
  57 #endif
  58 }
  59
  60 extern "C" {
  61   extern void pygive_(const char*, int);
  62   extern void pyupin_();
  63   extern void pyupev_();
  64   extern void pylist_(int&);
  65   extern void pystat_(int&);
  66 }
  67
  68 //==========================================================================
  69
  70 // Interfaces to the above routines, to make the C++ calls similar to Fortran.
  71 // This code section is generic.
  72
  73 class Pythia6Interface {
  74
  75 public:
  76
  77   // Give in a command to change a setting.
  78   static void pygive(const string cmnd) {
  79     const char* cstring = cmnd.c_str(); int len = cmnd.length();
  80     pygive_(cstring, len);
  81   }
  82
  83   // Initialize the generation for the given beam confiuration.
  84   static void pyinit(const string frame, const string beam,
  85     const string target, double wIn) {
  86     const char* cframe = frame.c_str(); int lenframe = frame.length();
  87     const char* cbeam = beam.c_str(); int lenbeam = beam.length();
  88     const char* ctarget = target.c_str(); int lentarget = target.length();
  89 #ifdef _WIN32
  90     pyinit_(cframe, lenframe, cbeam, lenbeam, ctarget, lentarget, wIn);
  91 #else
  92     pyinit_(cframe, cbeam, ctarget, wIn, lenframe, lenbeam, lentarget);
  93 #endif
  94   }
  95
  96   // Fill the initialization information in the HEPRUP commonblock.
  97   static void pyupin() {pyupin_();}
  98
  99   // Generate the next hard process and
 100   // fill the event information in the HEPEUP commonblock
 101   static void pyupev() {pyupev_();}
 102
 103   // List the event at the process level.
 104   static void pylist(int mode) {pylist_(mode);}
 105
 106   // Print statistics on the event generation process.
 107   static void pystat(int mode) {pystat_(mode);}
 108
 109 };
 110
 111 //==========================================================================
 112
 113 // Implement initialization fillHepRup method for Pythia6 example.
 114 // This code section is specific to the kind of precesses you generate.
 115
 116 // Of all parameters that could be set with pygive, only those that
 117 // influence the generation of the hard processes have any impact,
 118 // since this is the only part of the Fortran code that is used.
 119
 120 bool LHAupFortran::fillHepRup() {
 121
 122   // Set process to generate.
 123   // Example 1: QCD production; must set pTmin.
 124   //Pythia6Interface::pygive("msel = 1");
 125   //Pythia6Interface::pygive("ckin(3) = 20.");
 126   // Example 2: t-tbar production.
 127   //Pythia6Interface::pygive("msel = 6");
 128   // Example 3: Z0 production; must set mMin.
 129   Pythia6Interface::pygive("msel = 11");
 130   Pythia6Interface::pygive("ckin(1) = 50.");
 131
 132   // Speed up initialization: multiple interactions only in C++ code.
 133   Pythia6Interface::pygive("mstp(81)=0");
 134
 135   // Initialize for 14 TeV pp collider.
 136   Pythia6Interface::pyinit("cms","p","p",14000.);
 137
 138   // Fill initialization information in HEPRUP.
 139   Pythia6Interface::pyupin();
 140
 141   // Done.
 142   return true;
 143
 144 }
 145
 146 //==========================================================================
 147
 148 // Implement event generation fillHepEup method for Pythia 6 example.
 149 // This code section is generic.
 150
 151 bool LHAupFortran::fillHepEup() {
 152
 153   // Generate and fill the next Pythia6 event in HEPEUP.
 154   Pythia6Interface::pyupev();
 155
 156   // Done.
 157   return true;
 158
 159 }
 160
 161 //==========================================================================
 162
 163 // The main program.
 164 // This code section is specific to the physics study you want to do.
 165
 166 int main() {
 167
 168   // Generator. Shorthand for the event and for settings.
 169   Pythia pythia;
 170   Event& event = pythia.event;
 171   Settings& settings = pythia.settings;
 172
 173   // Set Pythia8 generation aspects.
 174   pythia.readString("BeamRemnants:primordialKThard = 2.");
 175   pythia.readString("MultipleInteractions:bProfile = 3");
 176
 177   // Initialize to access Pythia6 generator by Les Houches interface.
 178   LHAupFortran pythia6;
 179   pythia.init(&pythia6);
 180
 181   // Set some generation values.
 182   int nEvent = 100;
 183   int nList = 1;
 184   int nAbort = 10;
 185
 186   // List changed settings data.
 187   settings.listChanged();
 188
 189   // Histograms.
 190   double eCM = 14000.;
 191   double epTol = 1e-7 * eCM;
 192   Hist epCons("deviation from energy-momentum conservation",100,0.,epTol);
 193   Hist nFinal("final particle multiplicity",100,-0.5,1599.5);
 194   Hist nChg("final charged multiplicity",100,-0.5,799.5);
 195
 196   // Begin event loop.
 197   int iAbort = 0;
 198   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
 199
 200     // Generate events. Quit if too many failures.
 201     if (!pythia.next()) {
 202       if (++iAbort < nAbort) continue;
 203       cout << " Event generation aborted prematurely, owing to error!\n";
 204       break;
 205     }
 206
 207     // List first few events, both hard process and complete events.
 208     if (iEvent < nList) {
 209       pythia.info.list();
 210       // This call to Pythia6 is superfluous, but shows it can be done.
 211       Pythia6Interface::pylist(1);
 212       pythia.process.list();
 213       event.list();
 214     }
 215
 216     // Reset quantities to be summed over event.
 217     int nfin = 0;
 218     int nch = 0;
 219     Vec4 pSum = - (event[1].p() + event[2].p());
 220
 221     // Loop over final particles in the event.
 222     for (int i = 0; i < event.size(); ++i)
 223     if (event[i].isFinal()) {
 224       ++nfin;
 225       if (event[i].isCharged()) ++nch;
 226       pSum += event[i].p();
 227     }
 228
 229     // Fill summed quantities.
 230     double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
 231       + abs(pSum.pz());
 232     epCons.fill(epDev);
 233     nFinal.fill(nfin);
 234     nChg.fill(nch);
 235
 236   // End of event loop.
 237   }
 238
 239   // Final statistics. Must do call to Pythia6 explicitly.
 240   pythia.statistics();
 241   Pythia6Interface::pystat(1);
 242
 243   // Histogram output.
 244   cout << epCons << nFinal<< nChg;
 245
 246   // Done.
 247   return 0;
 248 }