PYTHIA8/pythia8145/examples/main29.cc

   1 // main29.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 // Illustration how to generate various two-body channels from
   7 // astroparticle processes, e.g. neutralino annihilation or decay.
   8 // To this end a "blob" of energy is created with unit cross section,
   9 // from the fictitious collision of two non-radiating incoming e+e-.
  10 // In the accompanying main29.cmnd file the decay channels of this
  11 // blob can be set up. Furthermore, only gamma, e+-, p/pbar and
  12 // neutrinos are stable, everything else is set to decay.
  13
  14 #include "Pythia.h"
  15
  16 using namespace Pythia8;
  17
  18 //==========================================================================
  19
  20 // A derived class for (e+ e- ->) GenericResonance -> various final states.
  21
  22 class Sigma1GenRes : public Sigma1Process {
  23
  24 public:
  25
  26   // Constructor.
  27   Sigma1GenRes() {}
  28
  29   // Evaluate sigmaHat(sHat): dummy unit cross section.
  30   virtual double sigmaHat() {return 1.;}
  31
  32   // Select flavour. No colour or anticolour.
  33   virtual void setIdColAcol() {setId( -11, 11, 999999);
  34     setColAcol( 0, 0, 0, 0, 0, 0);}
  35
  36   // Info on the subprocess.
  37   virtual string name()    const {return "GenericResonance";}
  38   virtual int    code()    const {return 9001;}
  39   virtual string inFlux()  const {return "ffbarSame";}
  40
  41 };
  42
  43 //==========================================================================
  44
  45 int main() {
  46
  47   // Pythia generator.
  48   Pythia pythia;
  49
  50   // A class to generate the fictitious resonance initial state.
  51   SigmaProcess* sigma1GenRes = new Sigma1GenRes();
  52
  53   // Hand pointer to Pythia.
  54   pythia.setSigmaPtr( sigma1GenRes);
  55
  56   // Read in the rest of the settings and data from a separate file.
  57   pythia.readFile("main29.cmnd");
  58
  59   // Initialization.
  60   pythia.init();
  61
  62   // Extract settings to be used in the main program.
  63   int    nEvent  = pythia.mode("Main:numberOfEvents");
  64   int    nList   = pythia.mode("Main:numberToList");
  65   int    nShow   = pythia.mode("Main:timesToShow");
  66   int    nAbort  = pythia.mode("Main:timesAllowErrors");
  67   bool   showCS  = pythia.flag("Main:showChangedSettings");
  68   bool   showCPD = pythia.flag("Main:showChangedParticleData");
  69
  70   // List changes.
  71   if (showCS)  pythia.settings.listChanged();
  72   if (showCPD) pythia.particleData.listChanged();
  73
  74   // Histogram particle spectra.
  75   Hist eGamma("energy spectrum of photons",        100, 0., 250.);
  76   Hist eE(    "energy spectrum of e+ and e-",      100, 0., 250.);
  77   Hist eP(    "energy spectrum of p and pbar",     100, 0., 250.);
  78   Hist eNu(   "energy spectrum of neutrinos",      100, 0., 250.);
  79   Hist eRest( "energy spectrum of rest particles", 100, 0., 250.);
  80
  81   // Begin event loop.
  82   int nPace  = max(1, nEvent / max(1, nShow) );
  83   int iAbort = 0;
  84   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
  85     if (nShow > 0 && iEvent%nPace == 0)
  86       cout << " Now begin event " << iEvent << endl;
  87
  88     // Generate events. Quit if many failures.
  89     if (!pythia.next()) {
  90       if (++iAbort < nAbort) continue;
  91       cout << " Event generation aborted prematurely, owing to error!\n";
  92       break;
  93     }
  94
  95     // List first few events.
  96     if (iEvent < nList) {
  97       pythia.info.list();
  98       pythia.event.list();
  99     }
 100
 101     // Loop over all particles and analyze the final-state ones.
 102     for (int i = 0; i < pythia.event.size(); ++i)
 103     if (pythia.event[i].isFinal()) {
 104       int idAbs = pythia.event[i].idAbs();
 105       double eI = pythia.event[i].e();
 106       if (idAbs == 22) eGamma.fill(eI);
 107       else if (idAbs == 11) eE.fill(eI);
 108       else if (idAbs == 2212) eP.fill(eI);
 109       else if (idAbs == 12 || idAbs == 14 || idAbs == 16) eNu.fill(eI);
 110       else {
 111         eRest.fill(eI);
 112         cout << " Error: stable id = " << pythia.event[i].id() << endl;
 113       }
 114     }
 115
 116   // End of event loop.
 117   }
 118
 119   // Final statistics and histograms.
 120   pythia.statistics();
 121   cout << eGamma << eE << eP << eNu << eRest;
 122
 123   // Done.
 124   return 0;
 125 }