runPythia.h

/*******************************************************************************

Example of a module which sees the JetscapeEvent and can read/modify it.

*******************************************************************************/

#include <string>
#include "Pythia8/Pythia.h"
using namespace Pythia8;


void runPythia(JetscapeEvent * js) {

  // Number of events, generated and listed ones.
  int nEvent    = 1;
  int nListJets = 5;

  // Generator. LHC process and output selection. Initialization.
  Pythia pythia;
  char numstr[21]; // enough to hold all numbers up to 64-bits
  sprintf(numstr, "%d", js->parameterSet.comEnergy);
  std::string comE = "Beams:eCM = ";
  comE += numstr;
  comE += ".";
  pythia.readString(comE);
  pythia.readString("HardQCD:all = on");
  pythia.readString("PhaseSpace:pTHatMin = 200.");
  pythia.readString("Next:numberShowInfo = 0");
  pythia.readString("Next:numberShowProcess = 0");
  pythia.readString("Next:numberShowEvent = 0");
  pythia.init();

  // Common parameters for the two jet finders.
  double etaMax   = 4.;
  double radius   = 0.7;
  double pTjetMin = 10.;
  // Exclude neutrinos (and other invisible) from study.
  int    nSel     = 2;
  // Range and granularity of CellJet jet finder.
  int    nEta     = 80;
  int    nPhi     = 64;

  // Set up SlowJet jet finder, with anti-kT clustering
  // and pion mass assumed for non-photons..
  SlowJet slowJet( -1, radius, pTjetMin, etaMax, nSel, 1);

  // Set up CellJet jet finder.
  CellJet cellJet( etaMax, nEta, nPhi, nSel);

  // Histograms. Note similarity in names, even when the two jet finders
  // do not calculate identically the same property (pT vs. ET, y vs. eta).
  Hist nJetsD("number of jets, CellJet - SlowJet", 45, -22.5, 22.5);

  // Begin event loop. Generate event. Skip if error.
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
    if (!pythia.next()) continue;

    // Analyze Slowet jet properties. List first few.
    slowJet. analyze( pythia.event );
    if (iEvent < nListJets) slowJet.list();

    // Fill SlowJet inclusive jet distributions.
    for (int i = 0; i < slowJet.sizeJet(); ++i) {
      js->jetCollection.push_back(Jet(0,slowJet.pT(i),slowJet.y(i),slowJet.phi(i)));
    }


    // Compare number of jets for the two finders.
    nJetsD.fill( cellJet.size() - slowJet.sizeJet() );

  // End of event loop. Statistics. Histograms.
  }
  pythia.stat();
  cout << nJetsD ;

  // Done.
  std::cout<<js<<std::endl;
}