source: HiSusy/trunk/Pythia8/pythia8170/src/MergingHooks.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 89.5 KB
Line 
1// MergingHooks.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2012 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// This file is written by Stefan Prestel.
7// Function definitions (not found in the header) for the HardProcess and
8// MergingHooks classes.
9
10#include "MergingHooks.h"
11
12namespace Pythia8 {
13 
14//==========================================================================
15
16// The HardProcess class.
17
18//--------------------------------------------------------------------------
19
20// Declaration of hard process class
21// This class holds information on the desired hard 2->2 process to be merged
22// This class is a container class for History class use.
23
24// Initialisation on the process string
25
26void HardProcess::initOnProcess( string process, ParticleData* particleData) {
27  state.init("(hard process)", particleData);
28  translateProcessString(process);
29}
30
31//--------------------------------------------------------------------------
32
33// Initialisation on the path to LHE file
34
35void HardProcess::initOnLHEF( string LHEfile, ParticleData* particleData) {
36  state.init("(hard process)", particleData);
37  translateLHEFString(LHEfile);
38}
39
40//--------------------------------------------------------------------------
41
42// Function to access the LHE file and read relevant information.
43// The merging scale will be read from the +1 jet sample, called
44//   LHEpath_1.lhe
45// while the hard process will be read from
46//   LHEpath_0.lhe
47// Currently, only read from MadEvent- or Sherpa-generated LHE files
48// is automatic, else the user is asked to supply the necessary
49// information.
50
51void HardProcess::translateLHEFString( string LHEpath){
52
53  // Open path to LHEF and extract merging scale
54  ifstream infile;
55  infile.open( (char*)( LHEpath +"_0.lhe").c_str());
56
57  // Check with ME generator has been used
58  int iLine =0;
59  int nLinesMax = 200;
60  string lineGenerator;
61  while( iLine < nLinesMax
62        && lineGenerator.find("SHERPA", 0) == string::npos
63        && lineGenerator.find("POWHEG-BOX", 0) == string::npos
64        && lineGenerator.find("Pythia8", 0) == string::npos
65        && lineGenerator.find("MadGraph", 0) == string::npos){
66    iLine++;
67    lineGenerator = " ";
68    getline(infile,lineGenerator);
69  }
70  infile.close();
71
72  vector <int> incom;
73  vector <int> inter;
74  vector <int> outgo;
75  // Particle identifiers, ordered in such a way that e.g. the "u"
76  // in a mu is not mistaken for an u quark
77  int inParticleNumbers[] = {
78        // Leptons
79        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
80        // Jet container
81        2212,2212,0,0,0,0,
82        // Quarks
83        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6};
84
85  string inParticleNamesSH[] = {
86        // Leptons
87        "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
88        // Jet container
89        "-93","93","-90","90","-91","91",
90        // Quarks
91        "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6"};
92  string inParticleNamesMG[] =  {
93        // Leptons
94        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
95        // Jet container
96        "p~","p","l+","l-","vl~","vl",
97        // Quarks
98        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};
99
100  // Declare intermediate particle identifiers
101  int interParticleNumbers[] = {
102         // Electroweak gauge bosons
103         22,23,-24,24,25,2400,
104         // Top quarks
105         -6,6,
106         // Dummy index as back-up
107         0,
108         // All squarks
109        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
110        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
111        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
112  // Declare names of intermediate particles
113  string interParticleNamesMG[] = {
114        // Electroweak gauge bosons
115        "a","z","w-","w+","h","W",
116         // Top quarks
117         "t~","t",
118        // Dummy index as back-up
119        "xx",
120         // All squarks
121         "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
122         "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"}; 
123
124  // Declare final state particle identifiers
125  int outParticleNumbers[] = {
126        // Leptons
127        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
128        // Jet container and lepton containers
129        2212,2212,0,0,0,0,1200,1100,5000,
130        // Quarks
131        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6,
132        // SM uncoloured bosons
133        22,23,-24,24,25,2400,
134        // Neutralino in SUSY
135        1000022,
136        // All squarks
137        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
138        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
139        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
140  // Declare names of final state particles
141  string outParticleNamesMG[] =  {
142        // Leptons
143        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
144        // Jet container and lepton containers
145        "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
146        // Quarks
147        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
148        // SM uncoloured bosons
149        "a","z","w-","w+","h","W",
150        // Neutralino in SUSY
151        "n1",
152        // All squarks
153        "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
154        "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};
155
156  string outParticleNamesSH[] = {
157        // Leptons
158        "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
159        // Jet container and lepton containers
160        "-93","93","-90","90","-91","91","0","0","0",
161        // Quarks
162        "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6",
163        // SM uncoloured bosons
164        "22","23","-24","24","25","0",
165        // Neutralino in SUSY
166        "1000022",
167        // All squarks
168        "-1000001","1000001","-1000002","1000002","-1000003","1000003",
169                   "-1000004","1000004",
170        "-1000005","1000005","-1000006","1000006","-2000001","2000001",
171                   "-2000002","2000002",
172        "-2000003","2000003","-2000004","2000004","-2000005","2000005",
173                   "-2000006","2000006"};
174
175  // Declare size of particle name arrays
176  int nIn   = 30;
177  int nInt  = 33;
178  int nOut  = 64;
179
180  // Save type of the generator, in order to be able to extract
181  // the tms definition
182  meGenType = (lineGenerator.find("MadGraph", 0) != string::npos)     ? -1
183              : (lineGenerator.find("SHERPA", 0) != string::npos)     ? -2
184              : (lineGenerator.find("POWHEG-BOX", 0) != string::npos) ? -3
185              : (lineGenerator.find("Pythia8", 0) != string::npos)    ? -4
186              : 0;
187
188  if (meGenType == -2){
189    // Now read merging scale
190    // Open path to LHEF and extract merging scale
191    infile.open( (char*)( LHEpath +"_1.lhe").c_str());
192    string lineTMS;
193    while(lineTMS.find("NJetFinder ", 0) == string::npos){
194      lineTMS = " ";
195      getline(infile,lineTMS);
196    }
197    infile.close();
198    lineTMS = lineTMS.substr(0,lineTMS.find(" 0.0 ",0));
199    lineTMS = lineTMS.substr(lineTMS.find(" ", 0)+3,lineTMS.size());
200    // Remove whitespaces
201    while(lineTMS.find(" ", 0) != string::npos)
202      lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
203    // Replace d with e
204    if ( lineTMS.find("d", 0) != string::npos)
205      lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
206    tms = atof((char*)lineTMS.c_str());
207
208    // Now read hard process
209    // Open path to LHEF and extract hard process
210    infile.open( (char*)( LHEpath +"_0.lhe").c_str());
211    string line;
212    while(line.find("Process", 0) == string::npos){
213      line = " ";
214      getline(infile,line);
215    }
216    infile.close();
217    line = line.substr(line.find(" ",0),line.size());
218
219    // Cut string into incoming and outgoing pieces
220    vector <string> pieces;
221    pieces.push_back( line.substr(0,line.find("->", 0)) );
222    // Do not count additional final jets
223    int end = (line.find("{", 0) != string::npos) ? line.find("{", 0)-2
224            : line.size();
225    pieces.push_back( line.substr(line.find("->", 0)+2, end) );
226
227    // Get incoming particles
228    for(int i=0; i < nIn; ++i) {
229      for(int n = pieces[0].find(inParticleNamesSH[i], 0);
230             n != int(string::npos);
231             n = pieces[0].find(inParticleNamesSH[i], n)) {
232        incom.push_back(inParticleNumbers[i]);
233        pieces[0].erase(pieces[0].begin()+n,
234                        pieces[0].begin()+n+inParticleNamesSH[i].size());
235        n=0;
236      }
237    }
238    // Get intermediate particles
239    // If intermediates are still empty, fill intermediate with default value
240    inter.push_back(0);
241    // Get final particles
242    for(int i=0; i < nOut; ++i) {
243      for(int n = pieces[1].find(outParticleNamesSH[i], 0);
244             n != int(string::npos);
245             n = pieces[1].find(outParticleNamesSH[i], n)) {
246        outgo.push_back(outParticleNumbers[i]);
247        pieces[1].erase(pieces[1].begin()+n,
248                        pieces[1].begin()+n+outParticleNamesSH[i].size());
249        n=0;
250      }
251    }
252
253  } else if (meGenType == -1 || meGenType == -3 || meGenType == -4){
254
255    // Now read merging scale
256    string lineTMS;
257
258    if (meGenType == -1) {
259      // Open path to LHEF and extract merging scale
260      infile.open( (char*)( LHEpath +"_1.lhe").c_str());
261      while(lineTMS.find("ktdurham", 0) == string::npos){
262        lineTMS = " ";
263        getline(infile,lineTMS);
264      }
265      lineTMS = lineTMS.substr(0,lineTMS.find("=",0));
266      infile.close();
267    } else {
268      lineTMS = "30.";
269    }
270
271    // Remove whitespaces
272    while(lineTMS.find(" ", 0) != string::npos)
273      lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
274    // Replace d with e
275    if ( lineTMS.find("d", 0) != string::npos)
276      lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
277    tms = atof((char*)lineTMS.c_str());
278
279    // Now read hard process
280    // Open path to LHEF and extract hard process
281    infile.open( (char*)( LHEpath +"_0.lhe").c_str());
282    string line;
283    while(line.find("@1", 0) == string::npos){
284      line = " ";
285      getline(infile,line);
286    }
287    infile.close();
288    line = line.substr(0,line.find("@",0));
289
290    // Count number of resonances
291    int appearances = 0;
292    for(int n = line.find("(", 0); n != int(string::npos);
293            n = line.find("(", n)) {
294      appearances++;
295      n++;
296    }
297
298    // Cut string in incoming, resonance+decay and outgoing pieces
299    vector <string> pieces;
300    for(int i =0; i < appearances;++i) {
301      int n = line.find("(", 0);
302      pieces.push_back(line.substr(0,n));
303      line = line.substr(n+1,line.size());
304    }
305    // Cut last resonance from rest
306    if (appearances > 0) {
307      pieces.push_back( line.substr(0,line.find(")",0)) );
308      pieces.push_back( line.substr(line.find(")",0)+1,line.size()) );
309    }
310
311    // If the string was not cut into pieces, i.e. no resonance was
312    // required, cut string using '>' as delimiter
313    if (pieces.empty() ){
314      appearances = 0;
315      for(int n = line.find(">", 0); n != int(string::npos);
316              n = line.find(">", n)) {
317        appearances++;
318        n++;
319      }
320
321      // Cut string in incoming and outgoing pieces
322      for(int i =0; i < appearances;++i) {
323        int n = line.find(">", 0);
324        pieces.push_back(line.substr(0,n));
325        line = line.substr(n+1,line.size());
326      }
327
328      if (appearances == 1) pieces.push_back(line);
329      if (appearances > 1) {
330        pieces.push_back( line.substr(0,line.find(">",0)) );
331        pieces.push_back( line.substr(line.find(">",0)+1,line.size()) );
332      }
333    }
334
335    // Get incoming particles
336    for(int i=0; i < nIn; ++i) {
337      for(int n = pieces[0].find(inParticleNamesMG[i], 0);
338             n != int(string::npos);
339             n = pieces[0].find(inParticleNamesMG[i], n)) {
340        incom.push_back(inParticleNumbers[i]);
341        pieces[0].erase(pieces[0].begin()+n,
342                        pieces[0].begin()+n+inParticleNamesMG[i].size());
343        n=0;
344      }
345    }
346
347    // Check intermediate resonances and decay products
348    for(int i =1; i < int(pieces.size()); ++i){
349      // Seperate strings into intermediate and outgoing, if not already done
350      int k = pieces[i].find(">", 0);
351
352      string intermediate = (pieces[i].find(">", 0) != string::npos) ? 
353                             pieces[i].substr(0,k) : "";
354      string outgoing = (pieces[i].find(">", 0) != string::npos) ? 
355                         pieces[i].substr(k+1,pieces[i].size()) : pieces[i];
356
357      // Get intermediate particles
358      for(int j=0; j < nInt; ++j) {
359        for(int n = intermediate.find(interParticleNamesMG[j], 0);
360               n != int(string::npos);
361               n = intermediate.find(interParticleNamesMG[j], n)) {
362          inter.push_back(interParticleNumbers[j]);
363          intermediate.erase(intermediate.begin()+n,
364                      intermediate.begin()+n+interParticleNamesMG[j].size());
365          n=0;
366        }
367      }
368
369      // Get outgoing particles
370      for(int j=0; j < nOut; ++j) {
371        for(int n = outgoing.find(outParticleNamesMG[j], 0);
372               n != int(string::npos);
373               n = outgoing.find(outParticleNamesMG[j], n)) {
374          outgo.push_back(outParticleNumbers[j]);
375          outgoing.erase(outgoing.begin()+n,
376                         outgoing.begin()+n+outParticleNamesMG[j].size());
377          n=0;
378        }
379      }
380
381      // For arbitrary or non-existing intermediate, remember zero for each
382      // two outgoing particles, without bosons.
383      if (inter.empty()) {
384
385        // For final state bosons, bookkeep the final state boson as
386        // intermediate as well.
387        int nBosons = 0;
388        for(int l=0; l < int(outgo.size()); ++l)
389          if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
390            nBosons++;
391
392        int nZeros = (outgo.size() - nBosons)/2;
393        for(int l=0; l < nZeros; ++l)
394          inter.push_back(0);
395      }
396
397      // For final state bosons, bookkeep the final state boson as
398      // intermediate as well.
399      for(int l=0; l < int(outgo.size()); ++l)
400        if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
401          inter.push_back(outgo[l]);
402
403    }
404
405  } else {
406
407    cout << "Reading of tms and hard process information from LHEF currently"
408         << " only automated for MadEvent- or SHERPA-produced LHEF" << endl;
409    int tempInt      = 0;
410    cout << "Use default process pp -> e+ve + jets? (0:no / 1:yes): ";
411    cin >> tempInt;
412    cout << endl;
413
414    if (tempInt == 0){
415      tempInt = 0;
416      double tempDouble  = 0.0;
417      cout << "Please specify merging scale (kT Durham, in GeV): ";
418      cin >> tempDouble;
419      tms = tempDouble;
420      meGenType = -1;
421      cout << endl;
422      cout << "Please specify first incoming particle ";
423      cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
424      cin >> tempInt;
425      incom.push_back(tempInt);
426      tempInt = 0;
427      cout << endl;
428      cout << "Please specify second incoming particle ";
429      cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
430      cin >> tempInt;
431      incom.push_back(tempInt);
432      tempInt = 0;
433      cout << endl;
434      cout << "Please specify intermediate particle, if any ";
435      cout << "(0 for none, else PDG code): ";
436      cin >> tempInt;
437      inter.push_back(tempInt);
438      cout << endl;
439      do {
440        tempInt = 0;
441        cout << "Please specify outgoing particle ";
442        cout << "(jet=2212, else PDG code, exit with 99): ";
443        cin >> tempInt;
444        if (tempInt != 99) outgo.push_back(tempInt);
445      } while(tempInt != 99);
446      cout << endl;
447    } else {
448      cout << "LHE file not produced by SHERPA or MG/ME - ";
449      cout << "Using default process and tms" << endl;
450      incom.push_back(2212);
451      incom.push_back(2212);
452      inter.push_back(24);
453      outgo.push_back(-11);
454      outgo.push_back(12);
455      tms = 10.;
456      meGenType = -1;
457    }
458  }
459
460  // Now store incoming, intermediate and outgoing
461  // Set intermediate tags
462  for(int i=0; i < int(inter.size()); ++i)
463    hardIntermediate.push_back(inter[i]);
464
465  // Set the incoming particle tags
466  if (incom.size() != 2)
467    cout << "Only two incoming particles allowed" << endl;
468  else {
469    hardIncoming1 = incom[0];
470    hardIncoming2 = incom[1];
471  }
472
473  // Remember final state bosons
474  int nBosons = 0;
475  for(int i=0; i < int(outgo.size()); ++i)
476    if ( (abs(outgo[i]) > 20 && abs(outgo[i]) <= 25) || outgo[i] == 2400)
477      nBosons++;
478  // Remember b-quark container
479  int nBQuarks = 0;
480  for(int i=0; i < int(outgo.size()); ++i)
481    if ( outgo[i] == 5000)
482      nBQuarks++;
483  // Remember jet container
484  int nJets = 0;
485  for(int i=0; i < int(outgo.size()); ++i)
486    if ( outgo[i] == 2212)
487      nJets++;
488  // Remember lepton container
489  int nLeptons = 0;
490  for(int i=0; i < int(outgo.size()); ++i)
491    if ( outgo[i] == 1100)
492      nLeptons++;
493  // Remember lepton container
494  int nNeutrinos = 0;
495  for(int i=0; i < int(outgo.size()); ++i)
496    if ( outgo[i] == 1200)
497      nNeutrinos++;
498  int nContainers = nLeptons + nNeutrinos + nJets + nBQuarks;
499
500  // Set final particle identifiers
501  if ( (outgo.size() - nBosons - nContainers)%2 == 1) {
502    cout << "Only even number of outgoing particles allowed" << endl; 
503    for(int i=0; i < int(outgo.size()); ++i)
504      cout << outgo[i] << endl;
505  } else {
506
507    // Push back particles / antiparticles
508    for(int i=0; i < int(outgo.size()); ++i)
509      if (outgo[i] > 0
510        && outgo[i] != 2212
511        && outgo[i] != 5000
512        && outgo[i] != 1100
513        && outgo[i] != 1200
514        && outgo[i] != 2400
515        && outgo[i] != 1000022)
516        hardOutgoing2.push_back( outgo[i]);
517      else if (outgo[i] < 0)
518        hardOutgoing1.push_back( outgo[i]);
519
520    // Save final state W-boson container as particle
521    for(int i=0; i < int(outgo.size()); ++i)
522      if ( outgo[i] == 2400)
523        hardOutgoing2.push_back( outgo[i]);
524
525    // Push back jets, distribute evenly amongst particles / antiparticles
526    // Push back majorana particles, distribute evenly
527    int iNow = 0;
528    for(int i=0; i < int(outgo.size()); ++i)
529      if ( (outgo[i] == 2212
530        || outgo[i] == 5000
531        || outgo[i] == 1200
532        || outgo[i] == 1000022)
533        && iNow%2 == 0 ){
534        hardOutgoing2.push_back( outgo[i]);
535        iNow++;
536      } else if ( (outgo[i] == 2212
537               || outgo[i] == 5000
538               || outgo[i] == 1100
539               || outgo[i] == 1000022)
540               && iNow%2 == 1 ){
541        hardOutgoing1.push_back( outgo[i]);
542        iNow++;
543      }
544  }
545
546  // Done
547}
548
549//--------------------------------------------------------------------------
550
551// Function to translate a string specitying the core process into the
552// internal notation
553// Currently, the input string has to be in MadEvent notation
554
555void HardProcess::translateProcessString( string process){
556
557  vector <int> incom;
558  vector <int> inter;
559  vector <int> outgo;
560  // Particle identifiers, ordered in such a way that e.g. the "u"
561  // in a mu is not mistaken for an u quark
562  int inParticleNumbers[] = {
563        // Leptons
564        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
565        // Jet container
566        2212,2212,0,0,0,0,
567        // Quarks
568        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6};
569  string inParticleNamesMG[] =  {
570        // Leptons
571        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
572        // Jet container
573        "p~","p","l+","l-","vl~","vl",
574        // Quarks
575        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};
576
577  // Declare intermediate particle identifiers
578  int interParticleNumbers[] = {
579         // Electroweak gauge bosons
580         22,23,-24,24,25,2400,
581         // Top quarks
582         -6,6,
583         // Dummy index as back-up
584         0,
585         // All squarks
586        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
587        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
588        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
589  // Declare names of intermediate particles
590  string interParticleNamesMG[] = {
591        // Electroweak gauge bosons
592        "a","z","w-","w+","h","W",
593         // Top quarks
594         "t~","t",
595        // Dummy index as back-up
596        "xx",
597         // All squarks
598         "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
599         "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"}; 
600
601  // Declare final state particle identifiers
602  int outParticleNumbers[] = {
603        // Leptons
604        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
605        // Jet container and lepton containers
606        2212,2212,0,0,0,0,1200,1100,5000,
607        // Quarks
608        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6,
609        // SM uncoloured bosons
610        22,23,-24,24,25,2400,
611        // Neutralino in SUSY
612        1000022,
613        // All squarks
614        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
615        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
616        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
617  // Declare names of final state particles
618  string outParticleNamesMG[] =  {
619        // Leptons
620        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
621        // Jet container and lepton containers
622        "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
623        // Quarks
624        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
625        // SM uncoloured bosons
626        "a","z","w-","w+","h","W",
627        // Neutralino in SUSY
628        "n1",
629        // All squarks
630        "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
631        "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};
632
633  // Declare size of particle name arrays
634  int nIn   = 30;
635  int nInt  = 33;
636  int nOut  = 64;
637
638  // Start mapping user-defined particles onto particle ids.
639  //string fullProc = "pp>{blaa,124}LEPTONS,NEUTRINOS";
640  string fullProc = process;
641
642  // Find user-defined hard process content
643  // Count number of user particles
644  int nUserParticles = 0;
645  for(int n = fullProc.find("{", 0); n != int(string::npos);
646          n = fullProc.find("{", n)) {
647    nUserParticles++;
648    n++;
649  }
650  // Cut user-defined particles from remaining process
651  vector <string> userParticleStrings;
652  for(int i =0; i < nUserParticles;++i) {
653    int n = fullProc.find("{", 0);
654    userParticleStrings.push_back(fullProc.substr(0,n));
655    fullProc = fullProc.substr(n+1,fullProc.size());
656  }
657  // Cut remaining process string from rest
658  if (nUserParticles > 0)
659    userParticleStrings.push_back(
660      fullProc.substr( 0, fullProc.find("}",0) ) );
661  // Remove curly brackets and whitespace
662  for(int i =0; i < int(userParticleStrings.size());++i) {
663    while(userParticleStrings[i].find("{", 0) != string::npos)
664      userParticleStrings[i].erase(userParticleStrings[i].begin()
665                                  +userParticleStrings[i].find("{", 0));
666    while(userParticleStrings[i].find("}", 0) != string::npos)
667      userParticleStrings[i].erase(userParticleStrings[i].begin()
668                                  +userParticleStrings[i].find("}", 0));
669    while(userParticleStrings[i].find(" ", 0) != string::npos)
670      userParticleStrings[i].erase(userParticleStrings[i].begin()
671                                  +userParticleStrings[i].find(" ", 0));
672  }
673  // Convert particle numbers in user particle to integers
674  vector<int>userParticleNumbers;
675  if ( int(userParticleStrings.size()) > 1) {
676    for( int i = 1; i < int(userParticleStrings.size()); ++i) {
677      userParticleNumbers.push_back( 
678        atoi((char*)userParticleStrings[i].substr(
679          userParticleStrings[i].find(",",0)+1, 
680          userParticleStrings[i].size()).c_str() ) );
681    }
682  }
683
684  // Save remaining process string
685  if (nUserParticles > 0)
686    userParticleStrings.push_back(
687      fullProc.substr(
688        fullProc.find("}",0)+1, fullProc.size() ) );
689  // Remove curly brackets and whitespace
690  for( int i = 0; i < int(userParticleStrings.size()); ++i) {
691    while(userParticleStrings[i].find("{", 0) != string::npos)
692      userParticleStrings[i].erase(userParticleStrings[i].begin()
693                                  +userParticleStrings[i].find("{", 0));
694    while(userParticleStrings[i].find("}", 0) != string::npos)
695      userParticleStrings[i].erase(userParticleStrings[i].begin()
696                                  +userParticleStrings[i].find("}", 0));
697    while(userParticleStrings[i].find(" ", 0) != string::npos)
698      userParticleStrings[i].erase(userParticleStrings[i].begin()
699                                  +userParticleStrings[i].find(" ", 0));
700  }
701
702  // Start mapping residual process string onto particle IDs.
703  // Declare leftover process after user-defined particles have been converted
704  string residualProc;
705  if ( int(userParticleStrings.size()) > 1 )
706    residualProc = userParticleStrings.front() + userParticleStrings.back();
707  else
708    residualProc = fullProc;
709
710  // Remove comma separation
711  while(residualProc.find(",", 0) != string::npos)
712    residualProc.erase(residualProc.begin()+residualProc.find(",",0));
713
714  // Count number of resonances
715  int appearances = 0;
716  for(int n = residualProc.find("(", 0); n != int(string::npos);
717          n = residualProc.find("(", n)) {
718    appearances++;
719    n++;
720  }
721
722  // Cut string in incoming, resonance+decay and outgoing pieces
723  vector <string> pieces;
724  for(int i =0; i < appearances;++i) {
725    int n = residualProc.find("(", 0);
726    pieces.push_back(residualProc.substr(0,n));
727    residualProc = residualProc.substr(n+1,residualProc.size());
728  }
729  // Cut last resonance from rest
730  if (appearances > 0) {
731    pieces.push_back( residualProc.substr(0,residualProc.find(")",0)) );
732    pieces.push_back( residualProc.substr(
733      residualProc.find(")",0)+1, residualProc.size()) );
734  }
735
736  // If the string was not cut into pieces, i.e. no resonance was
737  // required, cut string using '>' as delimiter
738  if (pieces.empty() ){
739    appearances = 0;
740    for(int n = residualProc.find(">", 0); n != int(string::npos);
741            n = residualProc.find(">", n)) {
742      appearances++;
743      n++;
744    }
745
746    // Cut string in incoming and outgoing pieces
747    for(int i =0; i < appearances;++i) {
748      int n = residualProc.find(">", 0);
749      pieces.push_back(residualProc.substr(0,n));
750      residualProc = residualProc.substr(n+1,residualProc.size());
751    }
752
753    if (appearances == 1) pieces.push_back(residualProc);
754    if (appearances > 1) {
755      pieces.push_back( residualProc.substr(0,residualProc.find(">",0)) );
756      pieces.push_back( residualProc.substr(
757        residualProc.find(">",0)+1, residualProc.size()) );
758    }
759  }
760
761  // Get incoming particles
762  for(int i=0; i < nIn; ++i) {
763    for(int n = pieces[0].find(inParticleNamesMG[i], 0);
764           n != int(string::npos);
765           n = pieces[0].find(inParticleNamesMG[i], n)) {
766      incom.push_back(inParticleNumbers[i]);
767      pieces[0].erase(pieces[0].begin()+n,
768                      pieces[0].begin()+n+inParticleNamesMG[i].size());
769      n=0;
770    }
771  }
772
773  // Check intermediate resonances and decay products
774  for(int i =1; i < int(pieces.size()); ++i){
775    // Seperate strings into intermediate and outgoing, if not already done
776    int k = pieces[i].find(">", 0);
777
778    string intermediate = (pieces[i].find(">", 0) != string::npos) ? 
779                           pieces[i].substr(0,k) : "";
780    string outgoing = (pieces[i].find(">", 0) != string::npos) ? 
781                       pieces[i].substr(k+1,pieces[i].size()) : pieces[i];
782
783    // Get intermediate particles
784    for(int j=0; j < nInt; ++j) {
785      for(int n = intermediate.find(interParticleNamesMG[j], 0);
786             n != int(string::npos);
787             n = intermediate.find(interParticleNamesMG[j], n)) {
788        inter.push_back(interParticleNumbers[j]);
789        intermediate.erase(intermediate.begin()+n,
790                    intermediate.begin()+n+interParticleNamesMG[j].size());
791        n=0;
792      }
793    }
794
795    // Get outgoing particles
796    for(int j=0; j < nOut; ++j) {
797      for(int n = outgoing.find(outParticleNamesMG[j], 0);
798             n != int(string::npos);
799             n = outgoing.find(outParticleNamesMG[j], n)) {
800        outgo.push_back(outParticleNumbers[j]);
801        outgoing.erase(outgoing.begin()+n,
802                       outgoing.begin()+n+outParticleNamesMG[j].size());
803        n=0;
804      }
805    }
806
807    // For arbitrary or non-existing intermediate, remember zero for each
808    // two outgoing particles, without bosons.
809    if (inter.empty()) {
810
811      // For final state bosons, bookkeep the final state boson as
812      // intermediate as well.
813      int nBosons = 0;
814      for(int l=0; l < int(outgo.size()); ++l)
815        if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
816          nBosons++;
817
818      int nZeros = (outgo.size() - nBosons)/2;
819      for(int l=0; l < nZeros; ++l)
820        inter.push_back(0);
821    }
822
823    // For final state bosons, bookkeep the final state boson as
824    // intermediate as well.
825    for(int l=0; l < int(outgo.size()); ++l)
826      if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
827        inter.push_back(outgo[l]);
828
829  }
830
831  // Now store incoming, intermediate and outgoing
832  // Set intermediate tags
833  for(int i=0; i < int(inter.size()); ++i)
834    hardIntermediate.push_back(inter[i]);
835
836  // Set the incoming particle tags
837  if (incom.size() != 2)
838    cout << "Only two incoming particles allowed" << endl;
839  else {
840    hardIncoming1 = incom[0];
841    hardIncoming2 = incom[1];
842  }
843
844  // Remember final state bosons
845  int nBosons = 0;
846  for(int i=0; i < int(outgo.size()); ++i)
847    if ( (abs(outgo[i]) > 20 && abs(outgo[i]) <= 25) || outgo[i] == 2400)
848      nBosons++;
849  // Remember b-quark container
850  int nBQuarks = 0;
851  for(int i=0; i < int(outgo.size()); ++i)
852    if ( outgo[i] == 5000)
853      nBQuarks++;
854  // Remember jet container
855  int nJets = 0;
856  for(int i=0; i < int(outgo.size()); ++i)
857    if ( outgo[i] == 2212)
858      nJets++;
859  // Remember lepton container
860  int nLeptons = 0;
861  for(int i=0; i < int(outgo.size()); ++i)
862    if ( outgo[i] == 1100)
863      nLeptons++;
864  // Remember lepton container
865  int nNeutrinos = 0;
866  for(int i=0; i < int(outgo.size()); ++i)
867    if ( outgo[i] == 1200)
868      nNeutrinos++;
869  int nContainers = nLeptons + nNeutrinos + nJets + nBQuarks;
870
871  // Set final particle identifiers
872  if ( (outgo.size() - nBosons - nContainers)%2 == 1) {
873    cout << "Only even number of outgoing particles allowed" << endl; 
874    for(int i=0; i < int(outgo.size()); ++i)
875      cout << outgo[i] << endl;
876  } else {
877
878    // Start with user-defined particles.
879    for( int i = 0; i < int(userParticleNumbers.size()); ++i)
880      if (userParticleNumbers[i] > 0) {
881        hardOutgoing2.push_back( userParticleNumbers[i]);
882        hardIntermediate.push_back(0);
883        // For non-existing intermediate, remember zero.
884      } else if (userParticleNumbers[i] < 0) {
885        hardOutgoing1.push_back( userParticleNumbers[i]);
886        // For non-existing intermediate, remember zero.
887        hardIntermediate.push_back(0);
888      }
889
890    // Push back particles / antiparticles
891    for(int i=0; i < int(outgo.size()); ++i)
892      if (outgo[i] > 0
893        && outgo[i] != 2212
894        && outgo[i] != 5000
895        && outgo[i] != 1100
896        && outgo[i] != 1200
897        && outgo[i] != 2400
898        && outgo[i] != 1000022)
899        hardOutgoing2.push_back( outgo[i]);
900      else if (outgo[i] < 0)
901        hardOutgoing1.push_back( outgo[i]);
902
903    // Save final state W-boson container as particle
904    for(int i=0; i < int(outgo.size()); ++i)
905      if ( outgo[i] == 2400)
906        hardOutgoing2.push_back( outgo[i]);
907
908    // Push back jets, distribute evenly among particles / antiparticles
909    // Push back majorana particles, distribute evenly
910    int iNow = 0;
911    for(int i=0; i < int(outgo.size()); ++i)
912      if ( (outgo[i] == 2212
913        || outgo[i] == 5000
914        || outgo[i] == 1200
915        || outgo[i] == 1000022)
916        && iNow%2 == 0 ){
917        hardOutgoing2.push_back( outgo[i]);
918        iNow++;
919      } else if ( (outgo[i] == 2212
920               || outgo[i] == 5000
921               || outgo[i] == 1100
922               || outgo[i] == 1000022)
923               && iNow%2 == 1 ){
924        hardOutgoing1.push_back( outgo[i]);
925        iNow++;
926      }
927  }
928
929  // Done
930}
931
932//--------------------------------------------------------------------------
933
934// Function to check if the candidates stored in Pos1 and Pos2, together with
935// a proposed candidate iPos are allowed.
936
937bool HardProcess::allowCandidates(int iPos, vector<int> Pos1, 
938  vector<int> Pos2, const Event& event){
939
940  bool allowed = true;
941
942  // Find colour-partner of new candidate
943  int type = (event[iPos].col() > 0) ? 1 : (event[iPos].acol() > 0) ? -1 : 0;
944
945  if (type == 0) return true;
946
947  if (type == 1){
948    int col = event[iPos].col();
949    int iPartner = 0;
950    for(int i=0; i < int(event.size()); ++i)
951      if ( i != iPos
952        && (( event[i].isFinal() && event[i].acol() == col)
953          ||( event[i].status() == -21 && event[i].col() == col) ))
954      iPartner = i;
955
956    vector<int> partners;
957    for(int i=0; i < int(event.size()); ++i)
958      for(int j=0; j < int(Pos1.size()); ++j)
959        if ( Pos1[j] != 0 && i != Pos1[j] && event[Pos1[j]].colType() != 0
960        && (( event[i].isFinal() && event[i].col() == event[Pos1[j]].acol())
961          ||( event[i].status() == -21
962           && event[i].acol() == event[Pos1[j]].acol()) ))
963         partners.push_back(i);
964
965    // Never allow equal initial partners!
966    if (event[iPartner].status() == -21){
967      for(int i=0; i < int(partners.size()); ++i)
968        if ( partners[i] == iPartner)
969          allowed = false;
970    }
971
972  } else {
973    int col = event[iPos].acol();
974    int iPartner = 0;
975    for(int i=0; i < int(event.size()); ++i)
976      if ( i != iPos
977        && (( event[i].isFinal() && event[i].col()  == col)
978          ||(!event[i].isFinal() && event[i].acol() == col) ))
979      iPartner = i;
980
981    vector<int> partners;
982    for(int i=0; i < int(event.size()); ++i)
983      for(int j=0; j < int(Pos2.size()); ++j)
984        if ( Pos2[j] != 0 && i != Pos2[j] && event[Pos2[j]].colType() != 0
985        && (( event[i].isFinal() && event[i].acol() == event[Pos2[j]].col())
986          ||( event[i].status() == -21
987           && event[i].col() == event[Pos2[j]].col()) ))
988         partners.push_back(i);
989
990
991    // Never allow equal initial partners!
992    if (event[iPartner].status() == -21){
993      for(int i=0; i < int(partners.size()); ++i){
994        if ( partners[i] == iPartner)
995          allowed = false;
996      }
997    }
998
999  }
1000
1001  return allowed;
1002
1003}
1004
1005//--------------------------------------------------------------------------
1006
1007// Function to identify the hard subprocess in the current event
1008
1009void HardProcess::storeCandidates( const Event& event, string process){
1010
1011  // Store the reference event
1012  state.clear();
1013  state = event;
1014
1015  // Local copy of intermediate bosons
1016  vector<int> intermediates;
1017  for(int i =0; i < int(hardIntermediate.size());++i)
1018    intermediates.push_back( hardIntermediate[i]);
1019
1020  // Local copy of outpoing partons
1021  vector<int> outgoing1;
1022  for(int i =0; i < int(hardOutgoing1.size());++i)
1023    outgoing1.push_back( hardOutgoing1[i]);
1024  vector<int> outgoing2;
1025  for(int i =0; i < int(hardOutgoing2.size());++i)
1026    outgoing2.push_back( hardOutgoing2[i]);
1027
1028  // Clear positions of intermediate and outgoing particles
1029  PosIntermediate.resize(0);
1030  PosOutgoing1.resize(0);
1031  PosOutgoing2.resize(0);
1032  for(int i =0; i < int(hardIntermediate.size());++i)
1033    PosIntermediate.push_back(0);
1034  for(int i =0; i < int(hardOutgoing1.size());++i)
1035    PosOutgoing1.push_back(0);
1036  for(int i =0; i < int(hardOutgoing2.size());++i)
1037    PosOutgoing2.push_back(0);
1038
1039  // For QCD dijet or e+e- > jets hard process, do not store any candidates,
1040  // as to not discrimintate clusterings
1041  if (  process.compare("pp>jj") == 0
1042    || process.compare("e+e->jj") == 0
1043    || process.compare("e+e->(z>jj)") == 0 ){
1044    for(int i =0; i < int(hardOutgoing1.size());++i)
1045      PosOutgoing1[i] = 0;
1046    for(int i =0; i < int(hardOutgoing2.size());++i)
1047      PosOutgoing2[i] = 0;
1048    // Done
1049    return;
1050  }
1051
1052  // Initialise vector of particles that were already identified as
1053  // hard process particles
1054  vector<int> iPosChecked;
1055
1056  // If the hard process is specified only by containers, then add all
1057  // particles matching with the containers to the hard process.
1058  bool hasOnlyContainers = true;
1059  for(int i =0; i < int(hardOutgoing1.size());++i)
1060    if (  hardOutgoing1[i] != 1100
1061      && hardOutgoing1[i] != 1200
1062      && hardOutgoing1[i] != 5000)
1063      hasOnlyContainers = false;
1064  for(int i =0; i < int(hardOutgoing2.size());++i)
1065    if (  hardOutgoing2[i] != 1100
1066      && hardOutgoing2[i] != 1200
1067      && hardOutgoing2[i] != 5000)
1068      hasOnlyContainers = false;
1069
1070  if (hasOnlyContainers){
1071
1072    PosOutgoing1.resize(0);
1073    PosOutgoing2.resize(0);
1074
1075    // Try to find all unmatched hard process leptons.
1076    // Loop through event to find outgoing lepton
1077    for(int i=0; i < int(event.size()); ++i){
1078
1079      // Skip non-final particles
1080      if ( !event[i].isFinal() ) continue;
1081
1082      // Skip all particles that have already been identified
1083      bool skip = false;
1084      for(int k=0; k < int(iPosChecked.size()); ++k){
1085        if (i == iPosChecked[k])
1086          skip = true;
1087      }
1088      if (skip) continue;
1089
1090      for(int j=0; j < int(outgoing2.size()); ++j){
1091
1092        // If the particle matches an outgoing neutrino, save it
1093        if ( outgoing2[j] == 1100
1094          && ( event[i].idAbs() == 11
1095            || event[i].idAbs() == 13
1096            || event[i].idAbs() == 15) ){
1097          PosOutgoing2.push_back(i);
1098          iPosChecked.push_back(i);
1099        }
1100
1101        // If the particle matches an outgoing lepton, save it
1102        if ( outgoing2[j] == 1200
1103          && ( event[i].idAbs() == 12
1104            || event[i].idAbs() == 14
1105            || event[i].idAbs() == 16) ){
1106          PosOutgoing2.push_back(i);
1107          iPosChecked.push_back(i);
1108        }
1109
1110        // If the particle matches an outgoing b-quark, save it
1111        if ( outgoing2[j] == 5000 && event[i].idAbs() == 5 ){
1112          PosOutgoing2.push_back(i);
1113          iPosChecked.push_back(i);
1114        }
1115
1116      }
1117
1118      // Skip all particles that have already been identified
1119      skip = false;
1120      for(int k=0; k < int(iPosChecked.size()); ++k){
1121        if (i == iPosChecked[k])
1122          skip = true;
1123      }
1124      if (skip) continue;
1125
1126      for(int j=0; j < int(outgoing1.size()); ++j){
1127
1128        // If the particle matches an outgoing neutrino, save it
1129        if ( outgoing1[j] == 1100
1130          && ( event[i].idAbs() == 11
1131            || event[i].idAbs() == 13
1132            || event[i].idAbs() == 15) ){
1133          PosOutgoing1.push_back(i);
1134          iPosChecked.push_back(i);
1135        }
1136
1137        // If the particle matches an outgoing lepton, save it
1138        if ( outgoing1[j] == 1200
1139          && ( event[i].idAbs() == 12
1140            || event[i].idAbs() == 14
1141            || event[i].idAbs() == 16) ){
1142          PosOutgoing1.push_back(i);
1143          iPosChecked.push_back(i);
1144        }
1145
1146        // If the particle matches an outgoing b-quark, save it
1147        if ( outgoing1[j] == 5000 && event[i].idAbs() == 5 ){
1148          PosOutgoing1.push_back(i);
1149          iPosChecked.push_back(i);
1150        }
1151
1152      }
1153    }
1154
1155    // Done
1156    return;
1157  }
1158
1159  // Now begin finding candidates when not only containers are used.
1160
1161  // First try to find final state bosons
1162  for(int i=0; i < int(intermediates.size()); ++i){
1163
1164    // Do nothing if the intermediate boson is absent
1165    if (intermediates[i] == 0) continue;
1166
1167    // Do nothing if this boson does not match any final state boson
1168    bool matchesFinalBoson = false;
1169    for(int j =0; j< int(outgoing1.size()); ++j){
1170      if ( intermediates[i] == outgoing1[j] )
1171        matchesFinalBoson = true;
1172    }
1173    for(int j =0; j< int(outgoing2.size()); ++j){
1174      if ( intermediates[i] == outgoing2[j] )
1175        matchesFinalBoson = true;
1176    }
1177    if (!matchesFinalBoson) continue;
1178
1179    // Loop through event
1180    for(int j=0; j < int(event.size()); ++j) {
1181      // If the particle has a requested intermediate id, check if
1182      // if is a final state boson
1183      if ( (event[j].id() == intermediates[i])
1184        ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
1185        PosIntermediate[i] = j;
1186        intermediates[i] = 0;
1187
1188        for(int k=0; k < int(outgoing1.size()); ++k) {
1189          if (event[j].id() == outgoing1[k]){
1190            PosOutgoing1[k] = j;
1191            outgoing1[k] = 99;
1192          }
1193        }
1194
1195        for(int k=0; k < int(outgoing2.size()); ++k) {
1196          if (event[j].id() == outgoing2[k]){
1197            PosOutgoing2[k] = j;
1198            outgoing2[k] = 99;
1199          }
1200        }
1201
1202        // Check for W-boson container
1203        for(int k=0; k < int(outgoing2.size()); ++k) {
1204          if (event[j].idAbs() == 24 && outgoing2[k] == 2400){
1205            PosOutgoing2[k] = j;
1206            outgoing2[k] = 99;
1207          }
1208        }
1209
1210        iPosChecked.push_back(j);
1211
1212      }
1213    }
1214  }
1215
1216  // Second try to find particles coupled to intermediate bosons
1217  for(int i=0; i < int(intermediates.size()); ++i){
1218
1219    // Do nothing if the intermediate boson is absent
1220    if (intermediates[i] == 0) continue;
1221
1222    // Loop through event
1223    for(int j=0; j < int(event.size()); ++j) {
1224      // If the particle has a requested intermediate id, check if
1225      // daughters are hard process particles
1226      if ( (event[j].id() == intermediates[i])
1227        ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
1228        // If this particle is a potential intermediate
1229        PosIntermediate[i] = j;
1230        intermediates[i] = 0;
1231        // If id's of daughters are good, store position
1232        int iPos1 = event[j].daughter1();
1233        int iPos2 = event[j].daughter2();
1234
1235        // Loop through daughters to check if these contain some hard
1236        // outgoing particles
1237        for( int k=iPos1; k <= iPos2; ++k){
1238          int id = event[k].id();
1239
1240          // Check if daughter is hard outgoing particle
1241          for(int l=0; l < int(outgoing2.size()); ++l)
1242            if ( outgoing2[l] != 99 ){
1243                // Found particle id
1244              if (id == outgoing2[l]
1245                // Found jet
1246                || (id > 0 && abs(id) < 10 && outgoing2[l] == 2212) ){
1247                // Store position
1248                PosOutgoing2[l] = k;
1249                // Remove the matched particle from the list
1250                outgoing2[l] = 99;
1251                iPosChecked.push_back(k);
1252                break;
1253              }
1254
1255            }
1256 
1257          // Check if daughter is hard outgoing antiparticle
1258          for(int l=0; l < int(outgoing1.size()); ++l)
1259            if ( outgoing1[l] != 99 ){
1260                // Found particle id
1261              if (id == outgoing1[l]
1262                // Found jet
1263                || (id < 0 && abs(id) < 10 && outgoing1[l] == 2212) ){
1264                // Store position
1265                PosOutgoing1[l] = k;
1266                // Remove the matched particle from the list
1267                outgoing1[l] = 99;
1268                iPosChecked.push_back(k);
1269                break;
1270            }
1271
1272          }
1273
1274        } // End loop through daughters
1275      } // End if ids match
1276    } // End loop through event
1277  } // End loop though requested intermediates
1278
1279  // If all outgoing particles were found, done
1280  bool done = true;
1281  for(int i=0; i < int(outgoing1.size()); ++i)
1282    if (outgoing1[i] != 99)
1283      done = false;
1284  for(int i=0; i < int(outgoing2.size()); ++i)
1285    if (outgoing2[i] != 99)
1286      done = false;
1287  // Return
1288  if (done) return;
1289
1290  // Leptons not associated with resonance are allowed.
1291  // Try to find all unmatched hard process leptons.
1292  // Loop through event to find outgoing lepton
1293  for(int i=0; i < int(event.size()); ++i){
1294    // Skip non-final particles and final partons
1295    if ( !event[i].isFinal() || event[i].colType() != 0)
1296      continue;
1297    // Skip all particles that have already been identified
1298    bool skip = false;
1299    for(int k=0; k < int(iPosChecked.size()); ++k){
1300      if (i == iPosChecked[k])
1301        skip = true;
1302    }
1303    if (skip) continue;
1304
1305    // Check if any hard outgoing leptons remain
1306    for(int j=0; j < int(outgoing2.size()); ++j){
1307      // Do nothing if this particle has already be found,
1308      // or if this particle is a jet or quark
1309      if (  outgoing2[j] == 99
1310        || outgoing2[j] == 2212
1311        || abs(outgoing2[j]) < 10)
1312        continue;
1313
1314      // If the particle matches an outgoing lepton, save it
1315      if (  event[i].id() == outgoing2[j] ){
1316        PosOutgoing2[j] = i;
1317        outgoing2[j] = 99;
1318        iPosChecked.push_back(i);
1319      }
1320    }
1321
1322    // Check if any hard outgoing antileptons remain
1323    for(int j=0; j < int(outgoing1.size()); ++j){
1324      // Do nothing if this particle has already be found,
1325      // or if this particle is a jet or quark
1326      if (  outgoing1[j] == 99
1327        || outgoing1[j] == 2212
1328        || abs(outgoing1[j]) < 10)
1329        continue;
1330
1331      // If the particle matches an outgoing lepton, save it
1332      if (event[i].id() == outgoing1[j] ){
1333        PosOutgoing1[j] = i;
1334        outgoing1[j] = 99;
1335        iPosChecked.push_back(i);
1336      }
1337    }
1338  }
1339
1340  multimap<int,int> out2copy;
1341  for(int i=0; i < int(event.size()); ++i)
1342    for(int j=0; j < int(outgoing2.size()); ++j)
1343      // Do nothing if this particle has already be found,
1344      // or if this particle is a jet, lepton container or lepton
1345      if (  outgoing2[j] != 99
1346        && outgoing2[j] != 2212
1347        && abs(outgoing2[j]) < 10
1348        && event[i].isFinal()
1349        && event[i].id() == outgoing2[j] ){
1350        out2copy.insert(make_pair(j, i));
1351      }
1352
1353  multimap<int,int> out1copy;
1354  for(int i=0; i < int(event.size()); ++i)
1355    for(int j=0; j < int(outgoing1.size()); ++j)
1356      // Do nothing if this particle has already be found,
1357      // or if this particle is a jet, lepton container or lepton
1358      if (  outgoing1[j] != 99
1359        && outgoing1[j] != 2212
1360        && abs(outgoing1[j]) < 10
1361        && event[i].isFinal()
1362        && event[i].id() == outgoing1[j] ){
1363        out1copy.insert(make_pair(j, i));
1364      }
1365
1366  if ( out1copy.size() >  out2copy.size()){
1367
1368    for ( multimap<int, int>::iterator it = out2copy.begin(); 
1369      it != out2copy.end(); ++it ) {
1370      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1371
1372        // Save parton
1373        PosOutgoing2[it->first] = it->second;
1374        // remove entry form lists
1375        outgoing2[it->first] = 99;
1376        iPosChecked.push_back(it->second);
1377
1378      }
1379    }
1380
1381    for ( multimap<int, int>::iterator it = out1copy.begin(); 
1382      it != out1copy.end(); ++it ) {
1383      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1384
1385        // Save parton
1386        PosOutgoing1[it->first] = it->second;
1387        // remove entry form lists
1388        outgoing1[it->first] = 99;
1389        iPosChecked.push_back(it->second);
1390
1391      }
1392    }
1393
1394  } else {
1395
1396    for ( multimap<int, int>::iterator it = out1copy.begin(); 
1397      it != out1copy.end(); ++it ) {
1398      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1399
1400        // Save parton
1401        PosOutgoing1[it->first] = it->second;
1402        // remove entry form lists
1403        outgoing1[it->first] = 99;
1404        iPosChecked.push_back(it->second);
1405
1406      }
1407    }
1408
1409    for ( multimap<int, int>::iterator it = out2copy.begin(); 
1410      it != out2copy.end(); ++it ) {
1411      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1412
1413        // Save parton
1414        PosOutgoing2[it->first] = it->second;
1415        // remove entry form lists
1416        outgoing2[it->first] = 99;
1417        iPosChecked.push_back(it->second);
1418
1419      }
1420    }
1421  }
1422
1423  // It sometimes happens that MadEvent does not put a
1424  // heavy coloured resonance into the LHE file, even if requested.
1425  // This means that the decay products of this resonance need to be
1426  // found separately.
1427  // Loop through event to find hard process (anti)quarks
1428  for(int i=0; i < int(event.size()); ++i){
1429
1430    // Skip non-final particles and final partons
1431    if ( !event[i].isFinal() || event[i].colType() == 0)
1432      continue;
1433
1434    // Skip all particles that have already been identified
1435    bool skip = false;
1436    for(int k=0; k < int(iPosChecked.size()); ++k){
1437      if (i == iPosChecked[k])
1438        skip = true;
1439    }
1440    if (skip) continue;
1441
1442    // Check if any hard outgoing quarks remain
1443    for(int j=0; j < int(outgoing2.size()); ++j){
1444      // Do nothing if this particle has already be found,
1445      // or if this particle is a jet, lepton container or lepton
1446      if (  outgoing2[j] == 99
1447        || outgoing2[j] == 2212
1448        || abs(outgoing2[j]) > 10)
1449        continue;
1450      // If the particle matches an outgoing quark, save it
1451      if (event[i].id() == outgoing2[j]){
1452        // Save parton
1453        PosOutgoing2[j] = i;
1454        // remove entry form lists
1455        outgoing2[j] = 99;
1456        iPosChecked.push_back(i);
1457      }
1458    }
1459
1460    // Check if any hard outgoing antiquarks remain
1461    for(int j=0; j < int(outgoing1.size()); ++j){
1462      // Do nothing if this particle has already be found,
1463      // or if this particle is a jet, lepton container or lepton
1464      if (  outgoing1[j] == 99
1465        || outgoing1[j] == 2212
1466        || abs(outgoing1[j]) > 10)
1467        continue;
1468      // If the particle matches an outgoing antiquark, save it
1469      if (event[i].id() == outgoing1[j]){
1470        // Save parton
1471        PosOutgoing1[j] = i;
1472        // Remove parton from list
1473        outgoing1[j] = 99;
1474        iPosChecked.push_back(i);
1475      }
1476    }
1477  }
1478
1479  // Done
1480}
1481
1482//--------------------------------------------------------------------------
1483
1484// Function to check if the particle event[iPos] matches any of
1485// the stored outgoing particles of the hard subprocess
1486
1487bool HardProcess::matchesAnyOutgoing(int iPos, const Event& event){
1488
1489  // Match quantum numbers of any first outgoing candidate
1490  bool matchQN1 = false;
1491  // Match quantum numbers of any second outgoing candidate
1492  bool matchQN2 = false;
1493  // Match parton in the hard process,
1494  // or parton from decay of electroweak boson in hard process,
1495  // or parton from decay of electroweak boson from decay of top
1496  bool matchHP = false;
1497
1498  // Check outgoing candidates
1499  for(int i=0; i < int(PosOutgoing1.size()); ++i)
1500    // Compare particle properties
1501    if ( event[iPos].id()         == state[PosOutgoing1[i]].id()
1502     && event[iPos].colType()    == state[PosOutgoing1[i]].colType() 
1503     && event[iPos].chargeType() == state[PosOutgoing1[i]].chargeType() 
1504     && event[iPos].col()        == state[PosOutgoing1[i]].col() 
1505     && event[iPos].acol()       == state[PosOutgoing1[i]].acol()
1506     && event[iPos].charge()     == state[PosOutgoing1[i]].charge() )
1507      matchQN1 = true;
1508
1509  // Check outgoing candidates
1510  for(int i=0; i < int(PosOutgoing2.size()); ++i)
1511    // Compare particle properties
1512    if ( event[iPos].id()         == state[PosOutgoing2[i]].id()
1513     && event[iPos].colType()    == state[PosOutgoing2[i]].colType() 
1514     && event[iPos].chargeType() == state[PosOutgoing2[i]].chargeType() 
1515     && event[iPos].col()        == state[PosOutgoing2[i]].col() 
1516     && event[iPos].acol()       == state[PosOutgoing2[i]].acol()
1517     && event[iPos].charge()     == state[PosOutgoing2[i]].charge() )
1518      matchQN2 = true;
1519
1520  // Check if maps to hard process:
1521  // Check that particle is in hard process
1522  if ( event[iPos].mother1()*event[iPos].mother2() == 12
1523      // Or particle has taken recoil from first splitting
1524      || (  event[iPos].status() == 44
1525         && event[event[iPos].mother1()].mother1()
1526           *event[event[iPos].mother1()].mother2() == 12 )
1527      // Or particle has on-shell resonace as mother
1528      || (  event[iPos].status() == 23
1529         && event[event[iPos].mother1()].mother1()
1530           *event[event[iPos].mother1()].mother2() == 12 )
1531      // Or particle has on-shell resonace as mother,
1532      // which again has and on-shell resonance as mother
1533      || (  event[iPos].status() == 23
1534         && event[event[iPos].mother1()].status() == -22
1535         && event[event[event[iPos].mother1()].mother1()].status() == -22
1536         && event[event[event[iPos].mother1()].mother1()].mother1()
1537           *event[event[event[iPos].mother1()].mother1()].mother2() == 12 ) )
1538      matchHP = true;
1539
1540  // Done
1541  return ( matchHP && (matchQN1 || matchQN2) );
1542
1543}
1544
1545
1546//--------------------------------------------------------------------------
1547
1548// Function to return the type of the ME generator
1549
1550int HardProcess::MEgenType(){ return meGenType;}
1551
1552//--------------------------------------------------------------------------
1553
1554// Function to get the number of coloured final state partons in the
1555// hard process
1556
1557int HardProcess::nQuarksOut(){
1558  int nFin =0;
1559  for(int i =0; i< int(hardOutgoing1.size()); ++i){
1560    if (hardOutgoing1[i] == 2212 || abs(hardOutgoing1[i]) < 10) nFin++;
1561  }
1562  for(int i =0; i< int(hardOutgoing2.size()); ++i){
1563    if (hardOutgoing2[i] == 2212 || abs(hardOutgoing2[i]) < 10) nFin++;
1564  }
1565  // For very loose hard process definition, check number of hard process
1566  // b-quarks explicitly.
1567  for(int i =0; i< int(hardOutgoing1.size()); ++i)
1568    if (hardOutgoing1[i] == 5000)
1569      for(int j =0; j< int(PosOutgoing1.size()); ++j)
1570        if (state[PosOutgoing1[j]].idAbs() == 5)
1571          nFin++;
1572  for(int i =0; i< int(hardOutgoing2.size()); ++i)
1573    if (hardOutgoing2[i] == 5000)
1574      for(int j =0; j< int(PosOutgoing2.size()); ++j)
1575        if (state[PosOutgoing2[j]].idAbs() == 5)
1576          nFin++;
1577  return nFin;
1578}
1579
1580//--------------------------------------------------------------------------
1581
1582// Function to get the number of uncoloured final state particles in the
1583// hard process
1584
1585int HardProcess::nLeptonOut(){
1586  int nFin =0;
1587  for(int i =0; i< int(hardOutgoing1.size()); ++i){
1588    if (abs(hardOutgoing1[i]) > 10 && abs(hardOutgoing1[i]) < 20) nFin++;
1589    // Bookkeep MSSM neutralinos as leptons
1590    if (abs(hardOutgoing1[i]) == 1000022) nFin++;
1591  }
1592  for(int i =0; i< int(hardOutgoing2.size()); ++i){
1593    if (abs(hardOutgoing2[i]) > 10 && abs(hardOutgoing2[i]) < 20) nFin++;
1594    // Bookkeep MSSM neutralinos as leptons
1595    if (abs(hardOutgoing2[i]) == 1000022) nFin++;
1596  }
1597  // For very loose hard process definition, check number of hard process
1598  // lepton explicitly.
1599  // Check lepton / neutrino containers as leptons
1600  for(int i =0; i< int(hardOutgoing1.size()); ++i)
1601    if (hardOutgoing1[i] == 1100)
1602      for(int j =0; j< int(PosOutgoing1.size()); ++j)
1603        if (  abs(state[PosOutgoing1[j]].id()) == 11
1604          || abs(state[PosOutgoing1[j]].id()) == 13
1605          || abs(state[PosOutgoing1[j]].id()) == 15 )
1606          nFin++;
1607  for(int i =0; i< int(hardOutgoing2.size()); ++i)
1608    if (hardOutgoing2[i] == 1200)
1609      for(int j =0; j< int(PosOutgoing2.size()); ++j)
1610        if (  abs(state[PosOutgoing2[j]].id()) == 12
1611          || abs(state[PosOutgoing2[j]].id()) == 14
1612          || abs(state[PosOutgoing2[j]].id()) == 16 )
1613          nFin++;
1614  return nFin;
1615}
1616
1617//--------------------------------------------------------------------------
1618
1619// Function to get the number of uncoloured final state particles in the
1620// hard process
1621
1622int HardProcess::nBosonsOut(){
1623  int nFin =0;
1624  for(int i =0; i< int(hardOutgoing1.size()); ++i){
1625    if (abs(hardOutgoing1[i]) > 20 && abs(hardOutgoing1[i]) <= 25) nFin++;
1626  }
1627  for(int i =0; i< int(hardOutgoing2.size()); ++i){
1628    if (abs(hardOutgoing2[i]) > 20 && abs(hardOutgoing2[i]) <= 25) nFin++;
1629    if ( hardOutgoing2[i] == 2400) nFin++;
1630  }
1631  return nFin;
1632}
1633
1634//--------------------------------------------------------------------------
1635
1636// Function to get the number of coloured initial state partons in the
1637// hard process
1638
1639int HardProcess::nQuarksIn(){
1640  int nIn =0;
1641  if (hardIncoming1 == 2212 || abs(hardIncoming1) < 10) nIn++;
1642  if (hardIncoming2 == 2212 || abs(hardIncoming2) < 10) nIn++;
1643  return nIn;
1644}
1645
1646//--------------------------------------------------------------------------
1647
1648// Function to get the number of uncoloured initial state particles in the
1649// hard process
1650
1651int HardProcess::nLeptonIn(){
1652  int nIn =0;
1653  if (abs(hardIncoming1) > 10 && abs(hardIncoming1) < 20) nIn++;
1654  if (abs(hardIncoming2) > 10 && abs(hardIncoming2) < 20) nIn++;
1655  return nIn;
1656}
1657
1658//--------------------------------------------------------------------------
1659
1660// Function to report if a resonace decay was found in the
1661// 2->2 hard sub-process in the current state
1662
1663int HardProcess::hasResInCurrent(){
1664  for(int i =0; i< int(PosIntermediate.size()); ++i)
1665    if (PosIntermediate[i] == 0) return 0;
1666  // Do not count final state bosons as resonaces
1667  for(int i =0; i< int(PosIntermediate.size()); ++i){
1668    for(int j =0; j< int(PosOutgoing1.size()); ++j){
1669      if ( PosIntermediate[i] == PosOutgoing1[j] )
1670        return 0;
1671    }
1672    for(int j =0; j< int(PosOutgoing2.size()); ++j){
1673      if ( PosIntermediate[i] == PosOutgoing2[j] )
1674        return 0;
1675    }
1676  }
1677  return 1;
1678}
1679
1680//--------------------------------------------------------------------------
1681
1682// Function to report the number of resonace decays in the 2->2 sub-process
1683// of the  current state
1684
1685int HardProcess::nResInCurrent(){
1686  int nRes = 0;
1687  for(int i =0; i< int(PosIntermediate.size()); ++i){
1688    if (PosIntermediate[i] != 0) {
1689      bool matchesFinalBoson = false;
1690      for(int j =0; j< int(PosOutgoing1.size()); ++j){
1691        if ( PosIntermediate[i] == PosOutgoing1[j] )
1692          matchesFinalBoson = true;
1693      }
1694      for(int j =0; j< int(PosOutgoing2.size()); ++j){
1695        if ( PosIntermediate[i] == PosOutgoing2[j] )
1696          matchesFinalBoson = true;
1697      }
1698      if (!matchesFinalBoson) nRes++;
1699    }
1700  }
1701  return nRes;
1702}
1703
1704//--------------------------------------------------------------------------
1705
1706// Function to report if a resonace decay was found in the
1707// 2->2 hard core process
1708
1709bool HardProcess::hasResInProc(){
1710
1711  for(int i =0; i< int(hardIntermediate.size()); ++i)
1712    if (hardIntermediate[i] == 0) return false;
1713  // Do not count final state bosons as resonaces
1714  for(int i =0; i< int(hardIntermediate.size()); ++i){
1715    for(int j =0; j< int(hardOutgoing1.size()); ++j){
1716      if ( hardIntermediate[i] == hardOutgoing1[j] )
1717        return false;
1718    }
1719    for(int j =0; j< int(hardOutgoing2.size()); ++j){
1720      if ( hardIntermediate[i] == hardOutgoing2[j] )
1721        return false;
1722    }
1723  }
1724  return true;
1725}
1726
1727//--------------------------------------------------------------------------
1728
1729// Function to print the hard process (for debug)
1730
1731void HardProcess::list() const {
1732  cout << "   Hard Process: ";
1733  cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
1734  cout << " \t -----> \t ";
1735  for(int i =0; i < int(hardIntermediate.size());++i)
1736    cout << hardIntermediate[i] << " "; 
1737  cout << " \t -----> \t ";
1738  for(int i =0; i < int(hardOutgoing1.size());++i)
1739    cout << hardOutgoing1[i] << " ";
1740  for(int i =0; i < int(hardOutgoing2.size());++i)
1741    cout << hardOutgoing2[i] << " ";
1742  cout << endl;
1743}
1744
1745//--------------------------------------------------------------------------
1746
1747// Function to list the hard process candiates in the matrix element state
1748// (for debug)
1749
1750void HardProcess::listCandidates() const {
1751  cout << "   Hard Process candidates: ";
1752  cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
1753  cout << " \t -----> \t ";
1754  for(int i =0; i < int(PosIntermediate.size());++i)
1755    cout << PosIntermediate[i] << " "; 
1756  cout << " \t -----> \t ";
1757  for(int i =0; i < int(PosOutgoing1.size());++i)
1758    cout << PosOutgoing1[i] << " ";
1759  for(int i =0; i < int(PosOutgoing2.size());++i)
1760    cout << PosOutgoing2[i] << " ";
1761  cout << endl;
1762}
1763
1764//--------------------------------------------------------------------------
1765
1766// Function to clear hard process information
1767
1768void HardProcess::clear() {
1769
1770  // Clear flavour of the first incoming particle
1771  hardIncoming1 = hardIncoming2 = 0;
1772  // Clear outgoing particles
1773  hardOutgoing1.resize(0);
1774  hardOutgoing2.resize(0);
1775  // Clear intermediate bosons in the hard 2->2
1776  hardIntermediate.resize(0);
1777
1778  // Clear reference event
1779  state.clear();
1780
1781  // Clear potential positions of outgoing particles in reference event
1782  PosOutgoing1.resize(0);
1783  PosOutgoing2.resize(0);
1784  // Clear potential positions of intermediate bosons in reference event
1785  PosIntermediate.resize(0);
1786
1787  // Clear merging scale read from LHE file
1788  tms = 0.;
1789  // Clear type of ME generator
1790  meGenType = 0;
1791
1792
1793}
1794
1795//==========================================================================
1796
1797// The MergingHooks class.
1798
1799//--------------------------------------------------------------------------
1800
1801// Initialise MergingHooks class
1802
1803void MergingHooks::init( Settings& settings, Info* infoPtrIn, 
1804  ParticleData* particleDataPtrIn, ostream& os){
1805
1806  // Save pointers
1807  infoPtr          = infoPtrIn;
1808  particleDataPtr  = particleDataPtrIn;
1809
1810  // Initialise AlphaS objects for reweighting
1811  double alphaSvalueFSR = settings.parm("TimeShower:alphaSvalue");
1812  int    alphaSorderFSR = settings.mode("TimeShower:alphaSorder");
1813  AlphaS_FSRSave.init(alphaSvalueFSR,alphaSorderFSR);
1814  double alphaSvalueISR = settings.parm("SpaceShower:alphaSvalue");
1815  int    alphaSorderISR = settings.mode("SpaceShower:alphaSorder");
1816  AlphaS_ISRSave.init(alphaSvalueISR,alphaSorderISR);
1817
1818  // Initialise AlphaS objects for reweighting
1819  int    alphaEMFSRorder = settings.mode("TimeShower:alphaEMorder");
1820  AlphaEM_FSRSave.init(alphaEMFSRorder, &settings);
1821
1822  // Initialise merging switches
1823  doUserMergingSave     =  settings.flag("Merging:doUserMerging");
1824
1825  // Initialise automated MadGraph kT merging
1826  doMGMergingSave       =  settings.flag("Merging:doMGMerging");
1827
1828  // Initialise kT merging
1829  doKTMergingSave       =  settings.flag("Merging:doKTMerging");
1830
1831  // Initialise evolution-pT merging
1832  doPTLundMergingSave   =  settings.flag("Merging:doPTLundMerging");
1833
1834  // Initialise \Delta_R_{ij}, pT_i Q_{ij} merging
1835  doCutBasedMergingSave =  settings.flag("Merging:doCutBasedMerging");
1836
1837  // Initialise exact definition of kT
1838  ktTypeSave            =  settings.mode("Merging:ktType");
1839
1840  // Get core process from user input
1841  processSave = settings.word("Merging:Process");
1842
1843  // Clear hard process
1844  hardProcess.clear();
1845
1846  bool doStandardMerging = doUserMergingSave || doKTMergingSave
1847         || doPTLundMergingSave || doCutBasedMergingSave;
1848
1849  // Initialise the hard process
1850  if ( doStandardMerging )
1851    hardProcess.initOnProcess(processSave, particleDataPtr);
1852  else
1853    hardProcess.initOnLHEF(lheInputFile, particleDataPtr);
1854
1855  // Parameters for reconstruction of evolution scales
1856  includeMassiveSave        = settings.flag("Merging:includeMassive");
1857  enforceStrongOrderingSave = settings.flag("Merging:enforceStrongOrdering");
1858  scaleSeparationFactorSave = settings.parm("Merging:scaleSeparationFactor");
1859  orderInRapiditySave       = settings.flag("Merging:orderInRapidity");
1860
1861  // Parameters for choosing history probabilistically
1862  nonJoinedNormSave     = settings.parm("Merging:nonJoinedNorm");
1863  fsrInRecNormSave      = settings.parm("Merging:fsrInRecNorm");
1864  pickByFullPSave       = settings.flag("Merging:pickByFullP");
1865  pickByPoPT2Save       = settings.flag("Merging:pickByPoPT2");
1866  includeRedundantSave  = settings.flag("Merging:includeRedundant");
1867
1868  // Parameters for scale choices
1869  unorderedScalePrescipSave   =
1870    settings.mode("Merging:unorderedScalePrescrip");
1871  unorderedASscalePrescipSave =
1872    settings.mode("Merging:unorderedASscalePrescrip");
1873  unorderedPDFscalePrescipSave =
1874    settings.mode("Merging:unorderedPDFscalePrescrip");
1875  incompleteScalePrescipSave  =
1876    settings.mode("Merging:incompleteScalePrescrip");
1877
1878  // Parameter for allowing swapping of one colour index while reclustering
1879  allowColourShufflingSave  =
1880    settings.flag("Merging:allowColourShuffling");
1881
1882  // Parameters to allow setting hard process scales to default (dynamical)
1883  // Pythia values.
1884  resetHardQRenSave     =  settings.flag("Merging:usePythiaQRenHard");
1885  resetHardQFacSave     =  settings.flag("Merging:usePythiaQFacHard");
1886
1887  // Parameters for choosing history by sum(|pT|)
1888  pickBySumPTSave       = settings.flag("Merging:pickBySumPT");
1889  herwigAcollFSRSave    = settings.parm("Merging:aCollFSR");
1890  herwigAcollISRSave    = settings.parm("Merging:aCollISR");
1891
1892  // Information on the shower cut-off scale
1893  pT0ISRSave            = settings.parm("SpaceShower:pT0Ref");
1894  pTcutSave             = settings.parm("SpaceShower:pTmin");
1895  pTcutSave             = max(pTcutSave,pT0ISRSave);
1896
1897  // Initialise CKKWL weight
1898  weightSave = 1.;
1899  weightCKKWLSave = 1.;
1900
1901  // Initialise merging scale
1902  tmsValueSave = 0.;
1903  tmsListSave.resize(0);
1904
1905  // Save merging scale on maximal number of jets
1906  if (  doKTMergingSave || doUserMergingSave || doPTLundMergingSave ) {
1907    // Read merging scale (defined in kT) from input parameter.
1908    tmsValueSave    = settings.parm("Merging:TMS");
1909    nJetMaxSave     = settings.mode("Merging:nJetMax");
1910  } else if (doMGMergingSave) {
1911    // Read merging scale (defined in kT) from LHE file.
1912    tmsValueSave    = hardProcess.tms;
1913    nJetMaxSave     = settings.mode("Merging:nJetMax");
1914  } else if (doCutBasedMergingSave) { 
1915    // Save list of cuts defining the merging scale.
1916    nJetMaxSave     = settings.mode("Merging:nJetMax");
1917    // Write tms cut values to list of cut values,
1918    // ordered by DeltaR_{ij}, pT_{i}, Q_{ij}.
1919    tmsListSave.resize(0);
1920    double drms     = settings.parm("Merging:dRijMS");
1921    double ptms     = settings.parm("Merging:pTiMS");
1922    double qms      = settings.parm("Merging:QijMS");
1923    tmsListSave.push_back(drms);
1924    tmsListSave.push_back(ptms);
1925    tmsListSave.push_back(qms);
1926  }
1927
1928  bool writeBanner =  doKTMergingSave || doMGMergingSave || doUserMergingSave
1929                   || doPTLundMergingSave || doCutBasedMergingSave;
1930
1931  if (!writeBanner) return;
1932
1933  // Write banner
1934  os << "\n"
1935     << " *---------- MEPS Merging Initialization  -----------------------*"
1936     << "\n";
1937  if ( doKTMergingSave || doMGMergingSave || doUserMergingSave
1938    || doPTLundMergingSave || doCutBasedMergingSave )
1939     os << " | CKKW-L tree-level merging:\n"
1940        << " | We merge  " << processSave << "  with up to  " << nJetMaxSave
1941        << " additional jets \n";
1942  if ( doKTMergingSave )
1943    os << " | Merging scale is defined in kT with the value ktMS = "
1944       << tmsValueSave << " GeV";
1945  else if ( doMGMergingSave )
1946    os << " | Perform automanted MG/ME merging \n"
1947       << " | Merging scale is defined in kT with the value ktMS = "
1948       << tmsValueSave << " GeV";
1949  else if ( doUserMergingSave )
1950    os << " | Merging scale is defined by the user, with the value tMS = "
1951       << tmsValueSave;
1952  else if ( doPTLundMergingSave )
1953    os << " | Merging scale is defined by Lund pT, with the value tMS = "
1954       << tmsValueSave;
1955  else if ( doCutBasedMergingSave )
1956    os << " | Merging scale is defined by combination of Delta R_{ij}, pT_i\n"
1957       << " | and Q_{ij} cut, with values \n"
1958       << " | Delta R_{ij,min} = " << tmsListSave[0] << "\n"
1959       << " | pT_{i,min} = " << tmsListSave[1] << "\n"
1960       << " | Q_{ij,min} = " << tmsListSave[2];
1961  os << "\n *---------- END MEPS Merging Initialization  -------------------*"
1962     << "\n\n";
1963
1964}
1965
1966//--------------------------------------------------------------------------
1967
1968// Function to return the number of clustering steps for the current event
1969
1970int MergingHooks::getNumberOfClusteringSteps(const Event& event){
1971
1972  // Count the number of final state partons
1973  int nFinalPartons = 0;
1974  for( int i=0; i < event.size(); ++i)
1975    if ( event[i].isFinal() && (event[i].isQuark() || event[i].isGluon()) )
1976      nFinalPartons++;
1977
1978  // Count the number of final state leptons
1979  int nFinalLeptons = 0;
1980  for( int i=0; i < event.size(); ++i)
1981    if ( event[i].isFinal() && event[i].isLepton())
1982      nFinalLeptons++;
1983
1984  // Add neutralinos to number of leptons
1985  for( int i=0; i < event.size(); ++i)
1986    if ( event[i].isFinal()
1987       && event[i].idAbs() == 1000022)
1988      nFinalLeptons++;
1989
1990  // Count the number of final state electroweak bosons
1991  int nFinalBosons = 0;
1992  for( int i=0; i < event.size(); ++i)
1993    if ( event[i].isFinal()
1994      && ( event[i].idAbs() == 22
1995        || event[i].idAbs() == 23
1996        || event[i].idAbs() == 24
1997        || event[i].idAbs() == 25 ) )
1998      nFinalBosons++;
1999
2000  // Save sum of all final state particles
2001  int nFinal = nFinalPartons + nFinalLeptons
2002             + 2*(nFinalBosons - nHardOutBosons() );
2003
2004  // Return the difference to the core process outgoing particles
2005  return (nFinal - nHardOutPartons() - nHardOutLeptons() );
2006}
2007
2008//--------------------------------------------------------------------------
2009
2010// Function to check if event contains an emission not present in the hard
2011// process.
2012
2013bool MergingHooks::isFirstEmission(const Event& event ) {
2014
2015  // Check that only one additional parton has been produced.
2016  // If not, we're already in the PS region (e.g. in MI).
2017  // Then, do not veto.
2018  int nMPI = infoPtr->nMPI();
2019  if (nMPI > 1) return false;
2020
2021  // Count particle types
2022  int nFinalQuarks   = 0;
2023  int nFinalGluons   = 0;
2024  int nFinalLeptons  = 0;
2025  int nFinalBosons   = 0;
2026  int nFinalPhotons  = 0;
2027  int nFinal         = 0;
2028  for( int i=0; i < event.size(); ++i) {
2029    if (event[i].isFinal()){
2030      if (  event[i].id() != 21
2031        && event[i].id() != 22
2032        && event[i].id() != 23
2033        && event[i].idAbs() != 24
2034        && event[i].id() != 25
2035        && event[i].colType() == 0)
2036        nFinalLeptons++;
2037      if (  event[i].id() == 23
2038        || event[i].idAbs() == 24
2039        || event[i].id() == 25)
2040        nFinalBosons++;
2041      if (  event[i].id() == 22)
2042        nFinalPhotons++;
2043      if ( event[i].isQuark())
2044        nFinalQuarks++;
2045      if ( event[i].isGluon())
2046        nFinalGluons++;
2047      if ( !event[i].isDiquark() )
2048        nFinal++;
2049    }
2050  }
2051
2052  // Return highest value if the event does not contain any final state
2053  // coloured particles.
2054  if (nFinalQuarks + nFinalGluons == 0) return false;
2055
2056  // Use MergingHooks functions to get information on the hard process.
2057  int nLeptons      = nHardOutLeptons();
2058
2059  // The state is already in the PS region if the number of leptons had been
2060  // increased bt QED splittings.
2061  if (nFinalLeptons > nLeptons) return false;
2062
2063  // If the mumber of photons if larger than in the hard process, QED
2064  // radiation has pushed the state into the PS region.
2065  int nPhotons = 0;
2066  for(int i =0; i< int(hardProcess.hardOutgoing1.size()); ++i)
2067    if (hardProcess.hardOutgoing1[i] == 22)
2068      nPhotons++;
2069  for(int i =0; i< int(hardProcess.hardOutgoing2.size()); ++i)
2070    if (hardProcess.hardOutgoing2[i] == 22)
2071      nPhotons++;
2072  if (nFinalPhotons > nPhotons) return false;
2073
2074  // Done
2075  return true;
2076}
2077
2078//--------------------------------------------------------------------------
2079
2080// Function to return the minimal kT in the event. If doKTMerging = true, this
2081// function will be used as a merging scale definition.
2082
2083double MergingHooks::kTms(const Event& event) {
2084
2085  // Only check first emission.
2086  if (!isFirstEmission(event)) return 0.;
2087
2088  // Find all electroweak decayed bosons in the state.
2089  vector<int> ewResonancePos;
2090  ewResonancePos.clear();
2091  for (int i=0; i < event.size(); ++i)
2092    if ( abs(event[i].status()) == 22 )
2093      ewResonancePos.push_back(i);
2094
2095  // Declare final parton vectors
2096  vector <int> FinalPartPos;
2097  FinalPartPos.clear();
2098  // Search inEvent record for final state partons.
2099  // Exclude decay products of ew resonance.
2100  for (int i=0; i < event.size(); ++i){
2101    if ( event[i].isFinal()
2102      && event[i].colType() != 0
2103      && event[i].idAbs()   != 6 ){
2104      bool isDecayProduct = false;
2105      for(int j=0; j < int(ewResonancePos.size()); ++j)
2106        if ( event.isAncestor(i, ewResonancePos[j])
2107          && event[i].mother2() == 0 )
2108          isDecayProduct = true;
2109      // Except for e+e- -> jets, do not check radiation in resonance decays.
2110      if ( !isDecayProduct
2111        || getProcessString().compare("e+e->jj") == 0 
2112        || getProcessString().compare("e+e->(z>jj)") == 0 )
2113        FinalPartPos.push_back(i);
2114    }
2115  }
2116
2117  // Declare kT algorithm parameters
2118  double Dparam = 0.4;
2119
2120  // Find minimal Durham kT in event, using own function: Check
2121  // definition of separation
2122  int type = (event[3].colType() == 0
2123           && event[4].colType() == 0) ? -1 : ktTypeSave;
2124  // Find minimal kT
2125  double ktmin = event[0].e();
2126  for(int i=0; i < int(FinalPartPos.size()); ++i){
2127    double kt12  = ktmin;
2128    // Compute separation to the beam axis for hadronic collisions
2129    if (type == 1 || type == 2) {
2130      double temp = event[FinalPartPos[i]].pT();
2131      kt12 = min(kt12, temp);
2132    }
2133    // Compute separation to other final state jets
2134    for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
2135      double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
2136                      type, Dparam);
2137      kt12 = min(kt12, temp);
2138    }
2139    // Keep the minimal Durham separation
2140    ktmin = min(ktmin,kt12);
2141  }
2142
2143  // Return minimal Durham kT
2144  return ktmin;
2145
2146}
2147
2148//--------------------------------------------------------------------------
2149
2150// Function to compute durham y separation from Particle input.
2151
2152double MergingHooks::kTdurham(const Particle& RadAfterBranch,
2153  const Particle& EmtAfterBranch, int Type, double D ){
2154
2155  // Declare return variable
2156  double ktdur;
2157  // Save 4-momenta of final state particles
2158  Vec4 jet1 = RadAfterBranch.p();
2159  Vec4 jet2 = EmtAfterBranch.p();
2160
2161  if ( Type == -1) {
2162    // Get angle between jets for e+e- collisions, make sure that
2163    // -1 <= cos(theta) <= 1
2164    double costh;
2165    if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
2166    else {
2167      costh = costheta(jet1,jet2);
2168    }
2169    // Calculate kt durham separation between jets for e+e- collisions
2170    ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
2171  } else if ( Type == 1 ){
2172    // Get delta_y for hadronic collisions:
2173    // Get mT of first jet
2174    double mT1sq = jet1.m2Calc() + jet1.pT2();
2175    double mT1 = 0.;
2176    if (mT1sq < 0) mT1 = - sqrt(-mT1sq);
2177    else mT1 = sqrt(mT1sq);
2178    // Get mT of second jet
2179    double mT2sq = jet2.m2Calc() + jet2.pT2();
2180    double mT2 = 0.;
2181    if (mT2sq < 0) mT2 = - sqrt(-mT2sq);
2182    else mT2 = sqrt(mT2sq);
2183    // Get rapidity of first jet
2184    double y1 = log( ( jet1.e() + abs(jet1.pz()) ) / mT1 );
2185    if (jet1.pz() < 0) y1 *= -1.;
2186    // Get rapidity of second jet
2187    double y2 = log( ( jet2.e() + abs(jet2.pz()) ) / mT2 );
2188    if (jet2.pz() < 0) y2 *= -1.;
2189    // Get delta_phi for hadronic collisions
2190    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
2191    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); 
2192    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
2193    double dPhi = acos( cosdPhi );
2194    // Calculate kT durham like fastjet,
2195    // but with rapidity instead of pseudo-rapidity
2196    ktdur = min( pow(pt1,2),pow(pt2,2) )
2197          * ( pow(y1-y2,2) + pow(dPhi,2) ) / pow(D,2);
2198  } else if ( Type == 2 ){
2199    // Get delta_eta for hadronic collisions
2200    double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
2201    double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
2202    // Get delta_phi and cos(Delta_phi) for hadronic collisions
2203    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
2204    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); 
2205    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
2206    double dPhi = acos( cosdPhi );
2207    // Calculate kT durham like fastjet
2208    ktdur = min( pow(pt1,2),pow(pt2,2) )
2209          * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
2210  } else if ( Type == 3 ){
2211    // Get cosh(Delta_eta) for hadronic collisions
2212    double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
2213    double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
2214    double coshdEta = cosh( eta1 - eta2 );
2215    // Get delta_phi and cos(Delta_phi) for hadronic collisions
2216    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
2217    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); 
2218    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
2219    // Calculate kT durham separation "SHERPA-like"
2220    ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
2221          * ( coshdEta - cosdPhi ) / pow(D,2);
2222  } else {
2223    ktdur = 0.0;
2224  }
2225  // Return kT
2226  return sqrt(ktdur);
2227}
2228
2229//--------------------------------------------------------------------------
2230
2231// Find the minimal Lund pT between coloured partons in the input
2232// event. If doPTLundMerging = true, this function will be used as a merging
2233// scale definition.
2234
2235double MergingHooks::rhoms( const Event& event, bool withColour){
2236
2237  // Only check first emission.
2238  if (!isFirstEmission(event)) return 0.;
2239
2240  // Find all electroweak decayed bosons in the state.
2241  vector<int> ewResonancePos;
2242  ewResonancePos.clear();
2243  for (int i=0; i < event.size(); ++i)
2244    if ( abs(event[i].status()) == 22 )
2245      ewResonancePos.push_back(i);
2246
2247  // Declare final parton vectors
2248  vector <int> FinalPartPos;
2249  FinalPartPos.clear();
2250  // Search inEvent record for final state partons.
2251  // Exclude decay products of ew resonance.
2252  for (int i=0; i < event.size(); ++i){
2253    if ( event[i].isFinal()
2254      && event[i].colType() != 0
2255      && event[i].idAbs()   != 6 ){
2256      bool isDecayProduct = false;
2257      for(int j=0; j < int(ewResonancePos.size()); ++j)
2258        if ( event.isAncestor(i, ewResonancePos[j])
2259          && event[i].mother2() == 0 )
2260          isDecayProduct = true;
2261      // Except for e+e- -> jets, do not check radiation in resonance decays.
2262      if ( !isDecayProduct
2263        || getProcessString().compare("e+e->jj") == 0 
2264        || getProcessString().compare("e+e->(z>jj)") == 0 )
2265        FinalPartPos.push_back(i);
2266    }
2267  }
2268
2269  // Get index of first incoming
2270  int in1 = 0;
2271  for (int i=0; i < event.size(); ++i)
2272    if (abs(event[i].status()) == 41 ){
2273      in1 = i;
2274      break;
2275    }
2276
2277  // Get index of second incoming
2278  int in2 = 0;
2279  for (int i=0; i < event.size(); ++i)
2280    if (abs(event[i].status()) == 42 ){
2281      in2 = i;
2282      break;
2283    }
2284
2285  // If no incoming of the cascade are found, try incoming
2286  if (in1 == 0 || in2 == 0){
2287    // Find current incoming partons
2288    for(int i=0; i < int(event.size()); ++i){
2289      if (event[i].mother1() == 1) in1 = i;
2290      if (event[i].mother1() == 2) in2 = i;
2291    }
2292  }
2293
2294  // Find minimal pythia pt in event
2295  double ptmin = event[0].e();
2296  for(int i=0; i < int(FinalPartPos.size()); ++i){
2297    double pt12  = ptmin;
2298    // Compute pythia ISR separation i-jet and first incoming
2299    if (event[in1].colType() != 0) {
2300      double temp = rhoPythia( event[in1],
2301                      event[FinalPartPos[i]], event[in2], -1 );
2302      pt12 = min(pt12, temp);
2303    }
2304    // Compute pythia ISR separation i-jet and second incoming
2305    if ( event[in2].colType() != 0) {
2306      double temp = rhoPythia( event[in2],
2307                      event[FinalPartPos[i]], event[in1], -1 );
2308      pt12 = min(pt12, temp);
2309    }
2310
2311    if (withColour) {
2312      // Compute pythia FSR separation between two jets,
2313      // with knowledge of colour connections
2314      for(int j=0; j < int(FinalPartPos.size()); ++j) {
2315
2316        // Find recoiler in event
2317        if ( i!=j ){
2318          bool isHard = false;
2319          int radAcl = event[FinalPartPos[i]].acol();
2320          int radCol = event[FinalPartPos[i]].col();
2321          int emtAcl = event[FinalPartPos[j]].acol();
2322          int emtCol = event[FinalPartPos[j]].col();
2323          int iRec = -1;
2324          // Check in final state
2325          if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
2326            iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
2327                     event,1, isHard);
2328          if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
2329            iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
2330                     event,1, isHard);
2331          if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
2332            iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
2333                     event,1, isHard);
2334          if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
2335            iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
2336                     event,1, isHard);
2337
2338          // Check in initial state
2339          if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
2340            iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
2341                     event,2, isHard);
2342          if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
2343            iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
2344                     event,2, isHard);
2345          if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
2346            iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
2347                     event,2, isHard);
2348          if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
2349            iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
2350                     event,2, isHard);
2351
2352          if (iRec > 0) {
2353            double temp = rhoPythia( event[FinalPartPos[i]],
2354                                     event[FinalPartPos[j]],
2355                                     event[iRec], 1 );
2356            pt12 = min(pt12, temp);
2357          }
2358        }
2359
2360        // If minimal pT below shower cut-off, return
2361        if (pt12 < 0.4) return pt12;
2362
2363      }
2364
2365    } else {
2366      // Compute pythia FSR separation between two jets,
2367      // without any knowledge of colour connections
2368      for(int j=0; j < int(FinalPartPos.size()); ++j) {
2369        for(int k=0; k < int(FinalPartPos.size()); ++k) {
2370          // Allow any parton as recoiler
2371          if ( (i != j) && (i != k) && (j != k) ){
2372
2373            double temp = 0.;
2374            // Only check splittings allowed by flavour, e.g.
2375            // only q -> qg and g -> qqbar
2376            temp = rhoPythia( event[FinalPartPos[i]],
2377                              event[FinalPartPos[j]],
2378                              event[FinalPartPos[k]], 1 );
2379            pt12 = min(pt12, temp);
2380          }
2381        }
2382        // If minimal pT below shower cut-off, return
2383        if (pt12 < 0.4) return pt12;
2384      }
2385
2386      // Compute pythia FSR separation between two jets, with initial recoiler
2387      // without any knowledge of colour connections
2388      if ( event[in1].colType() != 0 && event[in2].colType() != 0) {
2389        for(int j=0; j < int(FinalPartPos.size()); ++j) {
2390          // Allow both initial partons as recoiler
2391          if ( i != j ){
2392            // Check with first initial as recoiler
2393            double temp = rhoPythia( event[FinalPartPos[i]],
2394                                     event[FinalPartPos[j]],
2395                                     event[in1], 1 );
2396            pt12 = min(pt12, temp);
2397            // Check with second initial as recoiler
2398            temp        = rhoPythia( event[FinalPartPos[i]],
2399                                     event[FinalPartPos[j]],
2400                                     event[in2], 1 );
2401            pt12 = min(pt12, temp);
2402          }
2403
2404          // If minimal pT below shower cut-off, return
2405          if (pt12 < 0.4) return pt12;
2406        }
2407      }
2408
2409    }
2410    // Reset minimal y separation
2411    ptmin = min(ptmin,pt12);
2412  }
2413
2414  return ptmin;
2415
2416}
2417
2418//--------------------------------------------------------------------------
2419
2420// Function to compute "pythia pT separation" from Particle input, as a helper
2421// for rhoms(...).
2422
2423double MergingHooks::rhoPythia(const Particle& RadAfterBranch,
2424              const Particle& EmtAfterBranch,
2425              const Particle& RecAfterBranch, int ShowerType){
2426
2427  // Save type: 1 = FSR pT definition, else ISR definition
2428  int Type   = ShowerType;
2429  // Calculate virtuality of splitting
2430  int sign = (Type==1) ? 1 : -1;
2431  Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p());
2432  double Qsq = sign * Q.m2Calc();
2433  // Mass term of radiator
2434  double m2Rad = ( includeMassive()
2435               && abs(RadAfterBranch.id()) >= 4
2436               && abs(RadAfterBranch.id()) < 7)
2437               ? pow(particleDataPtr->m0(RadAfterBranch.id()), 2)
2438               : 0.;
2439  // Construct 2->3 variables for FSR
2440  Vec4   sum     = RadAfterBranch.p() + RecAfterBranch.p()
2441                 + EmtAfterBranch.p();
2442  double m2Dip = sum.m2Calc();
2443  double x1 = 2. * (sum * RadAfterBranch.p()) / m2Dip;
2444  double x3 = 2. * (sum * EmtAfterBranch.p()) / m2Dip;
2445  // Construct momenta of dipole before/after splitting for ISR
2446  Vec4 qBR(RadAfterBranch.p() - EmtAfterBranch.p() + RecAfterBranch.p());
2447  Vec4 qAR(RadAfterBranch.p() + RecAfterBranch.p());
2448  // Calculate z of splitting, different for FSR and ISR
2449  double z = (Type==1) ? x1/(x1+x3)
2450                     : (qBR.m2Calc())/( qAR.m2Calc());
2451  // Separation of splitting, different for FSR and ISR
2452  double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
2453  // pT^2 = separation*virtuality
2454  pTpyth *= (Qsq - sign*m2Rad);
2455  if (pTpyth < 0.) pTpyth = 0.;
2456  // Return pT
2457  return sqrt(pTpyth);
2458}
2459
2460//--------------------------------------------------------------------------
2461
2462// Function to find a colour (anticolour) index in the input event.
2463// Helper for rhoms
2464// IN  int col       : Colour tag to be investigated
2465//     int iExclude1 : Identifier of first particle to be excluded
2466//                     from search
2467//     int iExclude2 : Identifier of second particle to be excluded
2468//                     from  search
2469//     Event event   : event to be searched for colour tag
2470//     int type      : Tag to define if col should be counted as
2471//                      colour (type = 1) [->find anti-colour index
2472//                                         contracted with col]
2473//                      anticolour (type = 2) [->find colour index
2474//                                         contracted with col]
2475// OUT int           : Position of particle in event record
2476//                     contraced with col [0 if col is free tag]
2477
2478int MergingHooks::findColour(int col, int iExclude1, int iExclude2,
2479      const Event& event, int type, bool isHardIn){
2480
2481  bool isHard = isHardIn;
2482  int index = 0;
2483
2484  if (isHard){
2485    // Search event record for matching colour & anticolour
2486    for(int n = 0; n < event.size(); ++n) {
2487      if ( n != iExclude1 && n != iExclude2
2488        && event[n].colType() != 0
2489        &&(   event[n].status() > 0          // Check outgoing
2490           || event[n].status() == -21) ) {  // Check incoming
2491         if ( event[n].acol() == col ) {
2492          index = -n;
2493          break;
2494        }
2495        if ( event[n].col()  == col ){
2496          index =  n;
2497          break;
2498        }
2499      }
2500    }
2501  } else {
2502
2503    // Search event record for matching colour & anticolour
2504    for(int n = 0; n < event.size(); ++n) {
2505      if (  n != iExclude1 && n != iExclude2
2506        && event[n].colType() != 0
2507        &&(   event[n].status() == 43        // Check outgoing from ISR
2508           || event[n].status() == 51        // Check outgoing from FSR
2509           || event[n].status() == 52        // Check outgoing from FSR
2510           || event[n].status() == -41       // first initial
2511           || event[n].status() == -42) ) {  // second initial
2512        if ( event[n].acol() == col ) {
2513          index = -n;
2514          break;
2515        }
2516        if ( event[n].col()  == col ){
2517          index =  n;
2518          break;
2519        }
2520      }
2521    }
2522  }
2523  // if no matching colour / anticolour has been found, return false
2524  if ( type == 1 && index < 0) return abs(index);
2525  if ( type == 2 && index > 0) return abs(index);
2526
2527  return 0;
2528}
2529
2530//--------------------------------------------------------------------------
2531
2532// Function to check if the properties of the input particle should be
2533// checked against the cut-based merging scale defintion.
2534
2535bool MergingHooks::checkAgainstCut( const Particle& particle){
2536
2537  // Do not check uncoloured particles.
2538  if (particle.colType() == 0)
2539    return false;
2540  // Do not check tops and bottoms.
2541  if (particle.idAbs() == 5 || particle.idAbs() == 6)
2542    return false;
2543  // Done
2544  return true;
2545
2546}
2547
2548//--------------------------------------------------------------------------
2549
2550// Find the if the event passes the Delta R_{ij}, pT_{i} and Q_{ij} cuts on
2551// the matrix element, as needed for cut-based merging scale definition.
2552
2553double MergingHooks::cutbasedms( const Event& event ){
2554
2555  // Only check first emission.
2556  if (!isFirstEmission(event)) return -1.;
2557
2558  // Save allowed final state particles
2559  vector<int> partons;
2560  for( int i=0; i < event.size(); ++i)
2561    if ( event[i].isFinal() && checkAgainstCut(event[i]) )
2562      partons.push_back(i);
2563
2564  // Declare overall veto
2565  bool doVeto = false;
2566  // Declare vetoes
2567  bool vetoPT  = false;
2568  bool vetoRjj = false;
2569  bool vetoMjj = false;
2570  // Declare cuts used in matrix element
2571  double pTjmin = pTiMS();
2572  double mjjmin = QijMS();
2573  double rjjmin = dRijMS();
2574
2575  // Declare minimum values
2576  double minPT  = event[0].e();
2577  double minRJJ = 10.;
2578  double minMJJ = event[0].e();
2579
2580  // Check matrix element cuts
2581  for( int i=0; i < int(partons.size()); ++i){
2582    // Save pT value
2583    minPT = min(minPT,event[partons[i]].pT());
2584
2585    // Check two-parton cuts
2586    for( int j=0; j < int(partons.size()); ++j){
2587      if (i!=j){
2588
2589        // Save delta R value
2590        minRJJ = min(minRJJ, deltaRij( event[partons[i]].p(),
2591              event[partons[j]].p()) );
2592        // Save delta R value
2593        minMJJ = min(minMJJ, ( event[partons[i]].p()
2594                              +event[partons[j]].p() ).mCalc() );
2595
2596      }
2597    }
2598  // Done with cut evaluation
2599  }
2600
2601  // Check if all partons are in the matrix element region
2602  if (minPT  > pTjmin) vetoPT  = true;
2603  if (minRJJ > rjjmin) vetoRjj = true;
2604  if (minMJJ > mjjmin) vetoMjj = true;
2605
2606  // In the matrix element calculation, we have rejected the events if any of
2607  // the cuts had not been fulfilled,
2608  // i.e. we are in the matrix element domain if all cuts are fulfilled.
2609  // We veto any emission in the ME region.
2610  // Disregard the two-parton cuts if only one parton in the final state
2611  if (int(partons.size() == 1))
2612    doVeto = vetoPT;
2613  else
2614  // Veto if the combination of cuts would be in the ME region
2615    doVeto = vetoPT && vetoRjj && vetoMjj;
2616
2617  // If event is above merging scale, veto
2618  if (doVeto) return 1.;
2619
2620  // Else, do nothing
2621  return -1.;
2622
2623}
2624
2625//--------------------------------------------------------------------------
2626
2627// Function to compute Delta R separation from 4-vector input.
2628
2629double MergingHooks::deltaRij(Vec4 jet1, Vec4 jet2){
2630
2631  // Declare return variable
2632  double deltaR = 0.;
2633  // Get delta_eta and cosh(Delta_eta) for hadronic collisions
2634  double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
2635  double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
2636  // Get delta_phi and cos(Delta_phi) for hadronic collisions
2637  double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
2638  double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) ); 
2639  double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
2640  double dPhi = acos( cosdPhi );
2641  // Calculate kT durham like fastjet
2642  deltaR = sqrt(pow(eta1-eta2,2) + pow(dPhi,2));
2643  // Return kT
2644  return deltaR;
2645}
2646
2647//==========================================================================
2648
2649} // end namespace Pythia8
Note: See TracBrowser for help on using the repository browser.