main200
Back to index.
// main200.cc is a part of the PYTHIA event generator.
// Copyright (C) 2023 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// Simple example of the VINCIA (or DIRE) shower model(s), on Z decays at
// LEP I, with some basic event shapes, spectra, and multiplicity counts.
// Also useful as a basic test of the respective final-state showers.
// Authors:
// Peter Skands
// Keywords:
// Vincia
// Dire
// Electron‑positron
// Include Pythia8 header(s) and namespace.
#include "Pythia8/Pythia.h"
using namespace Pythia8;
// Main Program
int main() {
//************************************************************************
// Define Pythia 8 generator
Pythia pythia;
// Read user settings from file
pythia.readFile("main200.cmnd");
//************************************************************************
// Shorthands
Event& event = pythia.event;
Settings& settings = pythia.settings;
// Extract settings to be used in the main program.
int nEvent = settings.mode("Main:numberOfEvents");
int nAbort = settings.mode("Main:timesAllowErrors");
bool vinciaOn = settings.mode("PartonShowers:model") == 2;
bool helicityOn = vinciaOn && settings.flag("Vincia:helicityShower");
int iEventPri = -1;
double eCM = settings.parm("Beams:eCM");
//************************************************************************
// Initialize
if(!pythia.init()) { return EXIT_FAILURE; }
//************************************************************************
// Define a few PYTHIA utilities
Thrust Thr(1);
Sphericity SphLin(1,2);
//************************************************************************
// Define PYTHIA histograms
// Rootlink: to use ROOT instead, HIST -> TH1F ("tag",...)
// Number of quarks, partons, and charged particles
double nQsum = 0.0;
double nPsum = 0.0;
double nCsum = 0.0;
double nStrangeSum = 0.0;
double nCharmSum = 0.0;
double nBottomSum = 0.0;
double nBaryonSum = 0.0;
double nGammaSum = 0.0;
Hist histNQuarks("nQuarks", 50, -0.5, 49.5);
Hist histNPartons("nPartons", 100, -0.5, 99.5);
Hist histNCharged("nCharged", 100, -1.0, 199.0);
Hist histNGamma("nPhotons", 50, -0.5, 49.5);
Hist histXStrange("Ln(x) for strange quarks", 25, -5.0, 0.0);
Hist histXCharm("Ln(x) for charm quarks", 25, -5.0, 0.0);
Hist histXBottom("Ln(x) for bottom quarks", 25, -5.0, 0.0);
Hist histMSS("Invariant mass of s-sbar pairs",50,0.,10.);
Hist histMCC("Invariant mass of c-cbar pairs",50,0.,25.);
Hist histMBB("Invariant mass of b-bbar pairs",50,0.,50.);
Hist histX1("Ln(x) for 1st branching (QCD)",100,-7.5,0.0);
Hist histX1Gamma("Ln(x) for 1st branching (QED)",100,-7.5,0.0);
Hist histX2Gluon("Ln(x) for 2nd branching (gluons)",100,-7.5,0.0);
Hist histX2Quark("Ln(x) for 2nd branching (quarks)",100,-7.5,0.0);
// Thrust, C, and D parameters
int nBinsShapes=100;
Hist histT("1-T",nBinsShapes, 0.0, 0.5);
Hist histC("C",nBinsShapes,0.0,1.0);
Hist histD("D",nBinsShapes,0.0,1.0);
double wHistT=nBinsShapes/0.5;
double wHistCD=nBinsShapes/1.0;
//************************************************************************
// EVENT GENERATION LOOP.
// Generation, event-by-event printout, analysis, and histogramming.
// Counter for negative-weight events
double weight=1.0;
double sumWeights = 0.0;
// Begin event loop
int iAbort = 0;
for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
bool aborted = !pythia.next();
if(aborted){
event.list();
if (++iAbort < nAbort){
continue;
}
cout << " Event generation aborted prematurely, owing to error!\n";
cout<< "Event number was : "<<iEvent<<endl;
break;
}
// Check for weights
weight = pythia.info.weight();
sumWeights += weight;
// Print event with helicities
if (iEvent == 0 || iEvent == iEventPri)
if (helicityOn) event.list(helicityOn);
// Count FS charged hadrons, partons, and quarks
double nc=0;
int np=0;
int nq=0;
int nGam=0;
int nStrange=0;
int nCharm =0;
int nBottom =0;
int nStrangeBorn = (event[6].idAbs() == 3) ? 2 : 0;
int nCharmBorn = (event[6].idAbs() == 4) ? 2 : 0;
int nBottomBorn = (event[6].idAbs() == 5) ? 2 : 0;
vector<int> iStrange;
vector<int> iCharm;
vector<int> iBottom;
for (int i=1;i<event.size();i++) {
// Plot x distributions of first two emissions
if (i == 9) {
if (event[i].id() != 22)
histX1.fill(log(2*event[i].pAbs()/eCM),weight);
else
histX1Gamma.fill(log(2*event[i].pAbs()/eCM),weight);
}
if (i == 12) {
if (event[i].id() == 21)
histX2Gluon.fill(log(2*event[i].pAbs()/eCM),weight);
else if (event[i].idAbs() < 10)
histX2Quark.fill(log(2*event[i].pAbs()/eCM),weight);
}
// Count up final-state charged hadrons
if (event[i].isCharged() && event[i].status() > 80) nc++;
// Find last parton-level partons
int iDau1 = event[i].daughter1();
if (iDau1 == 0 || abs(event[iDau1].status()) > 80) {
// Count up partons and quarks
if (event[i].isQuark() || event[i].isGluon()) np++;
if (event[i].id() == 22) nGam++;
if (event[i].isQuark()) nq++;
if (event[i].idAbs() == 3) {
nStrange++;
histXStrange.fill(log(2*event[i].pAbs()/eCM),weight);
for (int j=0; j<(int)iStrange.size(); ++j) {
double m2qq = m2(event[i],event[iStrange[j]]);
histMSS.fill(sqrt(m2qq),weight);
}
iStrange.push_back(i);
}
if (event[i].idAbs() == 4) {
nCharm++;
histXCharm.fill(log(2*event[i].pAbs()/eCM),weight);
for (int j=0; j<(int)iCharm.size(); ++j) {
double m2qq = m2(event[i],event[iCharm[j]]);
histMCC.fill(sqrt(m2qq),weight);
}
iCharm.push_back(i);
}
if (event[i].idAbs() == 5) {
nBottom++;
histXBottom.fill(log(2*event[i].pAbs()/eCM),weight);
for (int j=0; j<(int)iBottom.size(); ++j) {
double m2qq = m2(event[i],event[iBottom[j]]);
histMBB.fill(sqrt(m2qq),weight);
}
iBottom.push_back(i);
}
}
if (event[i].isHadron() && event[i].isFinal()) {
int idAbs = abs(event[i].id());
if (idAbs > 1000 && idAbs < 10000) nBaryonSum++;
}
}
histNQuarks.fill(0.5*nq,weight);
histNPartons.fill(1.0*np,weight);
histNCharged.fill(1.0*nc,weight);
histNGamma.fill(1.0*nGam,weight);
nQsum += 0.5*nq * weight;
nPsum += np * weight;
nCsum += nc * weight;
nStrangeSum += (nStrange - nStrangeBorn) * weight;
nCharmSum += (nCharm - nCharmBorn) * weight;
nBottomSum += (nBottom - nBottomBorn) * weight;
nGammaSum += nGam * weight;
// Histogram thrust
Thr.analyze( event );
histT.fill(1.0-Thr.thrust(),wHistT*weight);
// Histogram Linear Sphericity values
if (np >= 2.0) {
SphLin.analyze( event );
double evC=3*(SphLin.eigenValue(1)*SphLin.eigenValue(2)
+ SphLin.eigenValue(2)*SphLin.eigenValue(3)
+ SphLin.eigenValue(3)*SphLin.eigenValue(1));
double evD=27*SphLin.eigenValue(1)*SphLin.eigenValue(2)
*SphLin.eigenValue(3);
histC.fill(evC,wHistCD*weight);
histD.fill(evD,wHistCD*weight);
}
}
//************************************************************************
// POST-RUN FINALIZATION
// Normalization, Statistics, Output.
//Normalize histograms to effective number of positive-weight events.
double normFac=1.0/sumWeights;
histT *= normFac;
histC *= normFac;
histD *= normFac;
histNQuarks *= normFac;
histNPartons *= normFac;
histNCharged *= normFac;
histNGamma *= normFac;
histXStrange *= normFac;
histXCharm *= normFac;
histXBottom *= normFac;
histX1 *= normFac;
histX1Gamma *= normFac;
histX2Gluon *= normFac;
histX2Quark *= normFac;
// Print a few histograms.
cout<<histNPartons<<endl;
if (nGammaSum > 0) cout<<histNGamma<<endl;
cout<<histNCharged<<endl;
cout<<histT<<endl;
cout<<histC<<endl;
cout<<histD<<endl;
cout<<histX1<<endl;
cout<<histX1Gamma<<endl;
cout<<histX2Gluon<<endl;
cout<<histX2Quark<<endl;
cout<<histXStrange<<endl;
cout<<histXCharm<<endl;
cout<<histXBottom<<endl;
cout<<histMSS<<endl;
cout<<histMCC<<endl;
cout<<histMBB<<endl;
// Print out end-of-run information.
pythia.stat();
cout<<endl;
cout<<fixed;
cout<< " <nGluonSplit> = "<<num2str(nQsum * normFac-1.0)<<endl;
cout<< " <nG->SS> = "<<num2str(nStrangeSum * normFac / 2.)<<endl;
cout<< " <nG->CC> = "<<num2str(nCharmSum * normFac / 2.)<<endl;
cout<< " <nG->BB> = "<<num2str(nBottomSum * normFac / 2.)<<endl;
cout<< " <nPartons> = "<<num2str(nPsum * normFac)<<endl;
cout<< " <nPhotons> = "<<num2str(nGammaSum * normFac)<<endl;
cout<< " <nCharged> = "<<num2str(nCsum * normFac)<<endl;
cout<< " <nBaryons> = "<<num2str(nBaryonSum * normFac)<<endl;
cout<<endl;
// Done.
return 0;
}