1 // main83.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2013 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.
6 // This program is written by Stefan Prestel.
7 // It illustrates how to do CKKW-L merging,
8 // see the Matrix Element Merging page in the online manual.
12 using namespace Pythia8;
14 // Functions for histogramming
15 #include "fastjet/PseudoJet.hh"
16 #include "fastjet/ClusterSequence.hh"
17 #include "fastjet/CDFMidPointPlugin.hh"
18 #include "fastjet/CDFJetCluPlugin.hh"
19 #include "fastjet/D0RunIIConePlugin.hh"
21 //==========================================================================
23 // Find the Durham kT separation of the clustering from
24 // nJetMin --> nJetMin-1 jets in te input event
26 double pTfirstJet( const Event& event, int nJetMin, double Rparam) {
28 double yPartonMax = 4.;
30 // Fastjet analysis - select algorithm and parameters
31 fastjet::Strategy strategy = fastjet::Best;
32 fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
33 fastjet::JetDefinition *jetDef = NULL;
34 // For hadronic collision, use hadronic Durham kT measure
35 if(event[3].colType() != 0 || event[4].colType() != 0)
36 jetDef = new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam,
37 recombScheme, strategy);
38 // For e+e- collision, use e+e- Durham kT measure
40 jetDef = new fastjet::JetDefinition(fastjet::ee_kt_algorithm,
41 recombScheme, strategy);
43 std::vector <fastjet::PseudoJet> fjInputs;
44 // Reset Fastjet input
47 // Loop over event record to decide what to pass to FastJet
48 for (int i = 0; i < event.size(); ++i) {
49 // (Final state && coloured+photons) only!
50 if ( !event[i].isFinal()
51 || event[i].isLepton()
52 || event[i].id() == 23
53 || abs(event[i].id()) == 24
54 || abs(event[i].y()) > yPartonMax)
57 // Store as input to Fastjet
58 fjInputs.push_back( fastjet::PseudoJet (event[i].px(),
59 event[i].py(), event[i].pz(),event[i].e() ) );
62 // Do nothing for empty input
63 if (int(fjInputs.size()) == 0) {
68 // Run Fastjet algorithm
69 fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
70 // Extract kT of first clustering
71 double pTFirst = sqrt(clustSeq.exclusive_dmerge_max(nJetMin-1));
79 //==========================================================================
81 // Class for user interaction with the merging
83 class MyMergingHooks : public MergingHooks {
89 // Default constructor
94 // Functional definition of the merging scale
95 virtual double tmsDefinition( const Event& event);
97 // Function to dampen weights calculated from histories with lowest
98 // multiplicity reclustered events that do not pass the ME cuts
99 virtual double dampenIfFailCuts( const Event& inEvent );
101 // Helper function for tms definition
102 double myKTdurham(const Particle& RadAfterBranch,
103 const Particle& EmtAfterBranch, int Type, double D );
107 //--------------------------------------------------------------------------
110 MyMergingHooks::MyMergingHooks() {}
113 MyMergingHooks::~MyMergingHooks() {}
115 //--------------------------------------------------------------------------
117 double MyMergingHooks::dampenIfFailCuts( const Event& inEvent ){
119 // Get pT for pure QCD 2->2 state
121 for( int i=0; i < inEvent.size(); ++i)
122 if(inEvent[i].isFinal() && inEvent[i].colType() != 0) {
123 pT = sqrt(pow(inEvent[i].px(),2) + pow(inEvent[i].py(),2));
127 // Veto history if lowest multiplicity event does not pass ME cuts
128 if(pT < 10.) return 0.;
134 //--------------------------------------------------------------------------
136 // Definition of the merging scale
138 double MyMergingHooks::tmsDefinition( const Event& event){
140 // Cut only on QCD partons!
141 // Count particle types
142 int nFinalColoured = 0;
144 for( int i=0; i < event.size(); ++i) {
145 if(event[i].isFinal()){
146 if(event[i].id() != 23 && abs(event[i].id()) != 24)
148 if( event[i].colType() != 0)
153 // Use MergingHooks in-built functions to get information on the hard process
154 int nLeptons = nHardOutLeptons();
155 int nQuarks = nHardOutPartons();
156 int nResNow = nResInCurrent();
158 // Check if photons, electrons etc. have been produced. If so, do not veto
159 if(nFinalNow - ( (nLeptons+nQuarks)/2 - nResNow)*2 != nFinalColoured){
160 // Sometimes, Pythia detaches the decay products even though no
161 // resonance was put into the LHE file, to catch this, add another
163 if(nFinalNow != nFinalColoured) return 0.;
166 // Check that one parton has been produced. If not (e.g. in MPI), do not veto
167 int nMPI = infoPtr->nMPI();
168 if(nMPI > 1) return 0.;
170 // Declare kT algorithm parameters
173 // Declare final parton vector
174 vector <int> FinalPartPos;
175 FinalPartPos.clear();
176 // Search event record for final state partons
177 for (int i=0; i < event.size(); ++i)
178 if(event[i].isFinal() && event[i].colType() != 0)
179 FinalPartPos.push_back(i);
181 // Find minimal Durham kT in event, using own function: Check
182 // definition of separation
183 int type = (event[3].colType() == 0 && event[4].colType() == 0) ? 1 : kTtype;
185 double ktmin = event[0].e();
186 for(int i=0; i < int(FinalPartPos.size()); ++i){
188 // Compute separation to the beam axis for hadronic collisions
189 if(type == -1 || type == -2) {
190 double temp = event[FinalPartPos[i]].pT();
191 kt12 = min(kt12, temp);
193 // Compute separation to other final state jets
194 for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
195 double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
197 kt12 = min(kt12, temp);
199 // Keep the minimal Durham separation
200 ktmin = min(ktmin,kt12);
203 // Return minimal Durham kT
208 //--------------------------------------------------------------------------
210 // Function to compute durham y separation from Particle input
212 double MyMergingHooks::myKTdurham(const Particle& RadAfterBranch,
213 const Particle& EmtAfterBranch, int Type, double D ){
215 // Declare return variable
217 // Save 4-momenta of final state particles
218 Vec4 jet1 = RadAfterBranch.p();
219 Vec4 jet2 = EmtAfterBranch.p();
222 // Get angle between jets for e+e- collisions, make sure that
223 // -1 <= cos(theta) <= 1
225 if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
227 costh = costheta(jet1,jet2);
229 // Calculate kt durham separation between jets for e+e- collisions
230 ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
231 } else if( Type == -1 ){
232 // Get delta_eta and cosh(Delta_eta) for hadronic collisions
233 double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
234 double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
235 // Get delta_phi and cos(Delta_phi) for hadronic collisions
236 double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
237 double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
238 double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
239 double dPhi = acos( cosdPhi );
240 // Calculate kT durham like fastjet
241 ktdur = min( pow(pt1,2),pow(pt2,2) )
242 * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
243 } else if( Type == -2 ){
244 // Get delta_eta and cosh(Delta_eta) for hadronic collisions
245 double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
246 double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
247 double coshdEta = cosh( eta1 - eta2 );
248 // Get delta_phi and cos(Delta_phi) for hadronic collisions
249 double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
250 double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
251 double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
252 // Calculate kT durham separation "SHERPA-like"
253 ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
254 * ( coshdEta - cosdPhi ) / pow(D,2);
262 //==========================================================================
264 // Example main programm to illustrate merging
266 int main( int argc, char* argv[] ){
268 // Check that correct number of command-line arguments
270 cerr << " Unexpected number of command-line arguments. \n You are"
271 << " expected to provide the arguments \n"
272 << " 1. Input file for settings \n"
273 << " 2. Full name of the input LHE file (with path) \n"
274 << " 3. Path for output histogram files \n"
275 << " Program stopped. " << endl;
282 // 1. Input file for settings
283 // 2. Path to input LHE file
284 // 3. OUtput histogram path
285 pythia.readFile(argv[1]);
286 string iPath = string(argv[2]);
287 string oPath = string(argv[3]);
290 int nEvent = pythia.mode("Main:numberOfEvents");
292 // Construct user inut for merging
293 MergingHooks* myMergingHooks = new MyMergingHooks();
294 pythia.setMergingHooksPtr( myMergingHooks );
296 // For ISR regularisation off
297 pythia.settings.forceParm("SpaceShower:pT0Ref",0.);
299 // Declare histograms
300 Hist histPTFirst("pT of first jet",100,0.,100.);
301 Hist histPTSecond("pT of second jet",100,0.,100.);
303 // Read in ME configurations
304 pythia.init(iPath,false);
306 if(pythia.flag("Main:showChangedSettings")) {
307 pythia.settings.listChanged();
310 // Start generation loop
311 for( int iEvent=0; iEvent<nEvent; ++iEvent ){
313 // Generate next event
314 if( ! pythia.next()) continue;
316 // Get CKKWL weight of current event
317 double weight = pythia.info.mergingWeight();
319 // Fill bins with CKKWL weight
320 double pTfirst = pTfirstJet(pythia.event,1, 0.4);
321 double pTsecnd = pTfirstJet(pythia.event,2, 0.4);
322 histPTFirst.fill( pTfirst, weight);
323 histPTSecond.fill( pTsecnd, weight);
325 if(iEvent%1000 == 0) cout << iEvent << endl;
327 } // end loop over events to generate
329 // print cross section, errors
332 // Normalise histograms
334 * pythia.info.sigmaGen()
335 * 1./ double(nEvent);
338 histPTSecond *= norm;
340 // Get the number of jets in the LHE file from the file name
341 string jetsInLHEF = iPath.substr(iPath.size()-5, iPath.size());
342 jetsInLHEF = jetsInLHEF.substr(0, jetsInLHEF.size()-4);
344 // Write histograms to dat file. Use "jetsInLHEF" to label the files
345 // Once all the samples up to the maximal desired jet multiplicity from the
346 // matrix element are run, add all histograms to produce a
347 // matrix-element-merged prediction
351 suffix << jetsInLHEF << "_wv.dat";
353 // Write histograms to file
354 write.open( (char*)(oPath + "PTjet1_" + suffix.str()).c_str());
355 histPTFirst.table(write);
358 write.open( (char*)(oPath + "PTjet2_" + suffix.str()).c_str());
359 histPTSecond.table(write);
363 delete myMergingHooks;