/usr/share/pythia8-examples/examples/main83.cc is in pythia8-examples 8.1.86-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 | // main83.cc is a part of the PYTHIA event generator.
// Copyright (C) 2014 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// This program is written by Stefan Prestel.
// It illustrates how to do CKKW-L merging,
// see the Matrix Element Merging page in the online manual.
#include "Pythia8/Pythia.h"
using namespace Pythia8;
// Functions for histogramming
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/CDFMidPointPlugin.hh"
#include "fastjet/CDFJetCluPlugin.hh"
#include "fastjet/D0RunIIConePlugin.hh"
//==========================================================================
// Find the Durham kT separation of the clustering from
// nJetMin --> nJetMin-1 jets in te input event
double pTfirstJet( const Event& event, int nJetMin, double Rparam) {
double yPartonMax = 4.;
// Fastjet analysis - select algorithm and parameters
fastjet::Strategy strategy = fastjet::Best;
fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
fastjet::JetDefinition *jetDef = NULL;
// For hadronic collision, use hadronic Durham kT measure
if(event[3].colType() != 0 || event[4].colType() != 0)
jetDef = new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam,
recombScheme, strategy);
// For e+e- collision, use e+e- Durham kT measure
else
jetDef = new fastjet::JetDefinition(fastjet::ee_kt_algorithm,
recombScheme, strategy);
// Fastjet input
std::vector <fastjet::PseudoJet> fjInputs;
// Reset Fastjet input
fjInputs.resize(0);
// Loop over event record to decide what to pass to FastJet
for (int i = 0; i < event.size(); ++i) {
// (Final state && coloured+photons) only!
if ( !event[i].isFinal()
|| event[i].isLepton()
|| event[i].id() == 23
|| abs(event[i].id()) == 24
|| abs(event[i].y()) > yPartonMax)
continue;
// Store as input to Fastjet
fjInputs.push_back( fastjet::PseudoJet (event[i].px(),
event[i].py(), event[i].pz(),event[i].e() ) );
}
// Do nothing for empty input
if (int(fjInputs.size()) == 0) {
delete jetDef;
return 0.0;
}
// Run Fastjet algorithm
fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
// Extract kT of first clustering
double pTFirst = sqrt(clustSeq.exclusive_dmerge_max(nJetMin-1));
delete jetDef;
// Return kT
return pTFirst;
}
//==========================================================================
// Class for user interaction with the merging
class MyMergingHooks : public MergingHooks {
private:
public:
// Default constructor
MyMergingHooks();
// Destructor
~MyMergingHooks();
// Functional definition of the merging scale
virtual double tmsDefinition( const Event& event);
// Function to dampen weights calculated from histories with lowest
// multiplicity reclustered events that do not pass the ME cuts
virtual double dampenIfFailCuts( const Event& inEvent );
// Helper function for tms definition
double myKTdurham(const Particle& RadAfterBranch,
const Particle& EmtAfterBranch, int Type, double D );
};
//--------------------------------------------------------------------------
// Constructor
MyMergingHooks::MyMergingHooks() {}
// Desctructor
MyMergingHooks::~MyMergingHooks() {}
//--------------------------------------------------------------------------
double MyMergingHooks::dampenIfFailCuts( const Event& inEvent ){
// Get pT for pure QCD 2->2 state
double pT = 0.;
for( int i=0; i < inEvent.size(); ++i)
if(inEvent[i].isFinal() && inEvent[i].colType() != 0) {
pT = sqrt(pow(inEvent[i].px(),2) + pow(inEvent[i].py(),2));
break;
}
// Veto history if lowest multiplicity event does not pass ME cuts
if(pT < 10.) return 0.;
return 1.;
}
//--------------------------------------------------------------------------
// Definition of the merging scale
double MyMergingHooks::tmsDefinition( const Event& event){
// Cut only on QCD partons!
// Count particle types
int nFinalColoured = 0;
int nFinalNow =0;
for( int i=0; i < event.size(); ++i) {
if(event[i].isFinal()){
if(event[i].id() != 23 && abs(event[i].id()) != 24)
nFinalNow++;
if( event[i].colType() != 0)
nFinalColoured++;
}
}
// Use MergingHooks in-built functions to get information on the hard process
int nLeptons = nHardOutLeptons();
int nQuarks = nHardOutPartons();
int nResNow = nResInCurrent();
// Check if photons, electrons etc. have been produced. If so, do not veto
if(nFinalNow - ( (nLeptons+nQuarks)/2 - nResNow)*2 != nFinalColoured){
// Sometimes, Pythia detaches the decay products even though no
// resonance was put into the LHE file, to catch this, add another
// if statement
if(nFinalNow != nFinalColoured) return 0.;
}
// Check that one parton has been produced. If not (e.g. in MPI), do not veto
int nMPI = infoPtr->nMPI();
if(nMPI > 1) return 0.;
// Declare kT algorithm parameters
double Dparam = 0.4;
int kTtype = -1;
// Declare final parton vector
vector <int> FinalPartPos;
FinalPartPos.clear();
// Search event record for final state partons
for (int i=0; i < event.size(); ++i)
if(event[i].isFinal() && event[i].colType() != 0)
FinalPartPos.push_back(i);
// Find minimal Durham kT in event, using own function: Check
// definition of separation
int type = (event[3].colType() == 0 && event[4].colType() == 0) ? 1 : kTtype;
// Find minimal kT
double ktmin = event[0].e();
for(int i=0; i < int(FinalPartPos.size()); ++i){
double kt12 = ktmin;
// Compute separation to the beam axis for hadronic collisions
if(type == -1 || type == -2) {
double temp = event[FinalPartPos[i]].pT();
kt12 = min(kt12, temp);
}
// Compute separation to other final state jets
for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
type, Dparam);
kt12 = min(kt12, temp);
}
// Keep the minimal Durham separation
ktmin = min(ktmin,kt12);
}
// Return minimal Durham kT
return ktmin;
}
//--------------------------------------------------------------------------
// Function to compute durham y separation from Particle input
double MyMergingHooks::myKTdurham(const Particle& RadAfterBranch,
const Particle& EmtAfterBranch, int Type, double D ){
// Declare return variable
double ktdur;
// Save 4-momenta of final state particles
Vec4 jet1 = RadAfterBranch.p();
Vec4 jet2 = EmtAfterBranch.p();
if( Type == 1) {
// Get angle between jets for e+e- collisions, make sure that
// -1 <= cos(theta) <= 1
double costh;
if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
else {
costh = costheta(jet1,jet2);
}
// Calculate kt durham separation between jets for e+e- collisions
ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
} else if( Type == -1 ){
// Get delta_eta and cosh(Delta_eta) for hadronic collisions
double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
// Get delta_phi and cos(Delta_phi) for hadronic collisions
double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
double dPhi = acos( cosdPhi );
// Calculate kT durham like fastjet
ktdur = min( pow(pt1,2),pow(pt2,2) )
* ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
} else if( Type == -2 ){
// Get delta_eta and cosh(Delta_eta) for hadronic collisions
double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
double coshdEta = cosh( eta1 - eta2 );
// Get delta_phi and cos(Delta_phi) for hadronic collisions
double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
// Calculate kT durham separation "SHERPA-like"
ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
* ( coshdEta - cosdPhi ) / pow(D,2);
} else {
ktdur = 0.0;
}
// Return kT
return sqrt(ktdur);
}
//==========================================================================
// Example main programm to illustrate merging
int main( int argc, char* argv[] ){
// Check that correct number of command-line arguments
if (argc != 4) {
cerr << " Unexpected number of command-line arguments. \n You are"
<< " expected to provide the arguments \n"
<< " 1. Input file for settings \n"
<< " 2. Full name of the input LHE file (with path) \n"
<< " 3. Path for output histogram files \n"
<< " Program stopped. " << endl;
return 1;
}
Pythia pythia;
// Input parameters:
// 1. Input file for settings
// 2. Path to input LHE file
// 3. OUtput histogram path
pythia.readFile(argv[1]);
string iPath = string(argv[2]);
string oPath = string(argv[3]);
// Number of events
int nEvent = pythia.mode("Main:numberOfEvents");
// Construct user inut for merging
MergingHooks* myMergingHooks = new MyMergingHooks();
pythia.setMergingHooksPtr( myMergingHooks );
// For ISR regularisation off
pythia.settings.forceParm("SpaceShower:pT0Ref",0.);
// Declare histograms
Hist histPTFirst("pT of first jet",100,0.,100.);
Hist histPTSecond("pT of second jet",100,0.,100.);
// Read in ME configurations
pythia.init(iPath,false);
if(pythia.flag("Main:showChangedSettings")) {
pythia.settings.listChanged();
}
// Start generation loop
for( int iEvent=0; iEvent<nEvent; ++iEvent ){
// Generate next event
if( ! pythia.next()) continue;
// Get CKKWL weight of current event
double weight = pythia.info.mergingWeight();
// Fill bins with CKKWL weight
double pTfirst = pTfirstJet(pythia.event,1, 0.4);
double pTsecnd = pTfirstJet(pythia.event,2, 0.4);
histPTFirst.fill( pTfirst, weight);
histPTSecond.fill( pTsecnd, weight);
if(iEvent%1000 == 0) cout << iEvent << endl;
} // end loop over events to generate
// print cross section, errors
pythia.stat();
// Normalise histograms
double norm = 1.
* pythia.info.sigmaGen()
* 1./ double(nEvent);
histPTFirst *= norm;
histPTSecond *= norm;
// Get the number of jets in the LHE file from the file name
string jetsInLHEF = iPath.substr(iPath.size()-5, iPath.size());
jetsInLHEF = jetsInLHEF.substr(0, jetsInLHEF.size()-4);
// Write histograms to dat file. Use "jetsInLHEF" to label the files
// Once all the samples up to the maximal desired jet multiplicity from the
// matrix element are run, add all histograms to produce a
// matrix-element-merged prediction
ofstream write;
stringstream suffix;
suffix << jetsInLHEF << "_wv.dat";
// Write histograms to file
write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str());
histPTFirst.table(write);
write.close();
write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str());
histPTSecond.table(write);
write.close();
delete myMergingHooks;
return 0;
// Done
}
|