Line data Source code
1 : // MergingHooks.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "Pythia8/MergingHooks.h"
11 :
12 : namespace 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 :
26 : void HardProcess::initOnProcess( string process, ParticleData* particleData) {
27 0 : state.init("(hard process)", particleData);
28 0 : translateProcessString(process);
29 0 : }
30 :
31 : //--------------------------------------------------------------------------
32 :
33 : // Initialisation on the path to LHE file
34 :
35 : void HardProcess::initOnLHEF( string LHEfile, ParticleData* particleData) {
36 0 : state.init("(hard process)", particleData);
37 0 : translateLHEFString(LHEfile);
38 0 : }
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 :
51 : void HardProcess::translateLHEFString( string LHEpath){
52 :
53 : // Open path to LHEF and extract merging scale
54 0 : ifstream infile;
55 0 : 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 0 : string lineGenerator;
61 0 : while( iLine < nLinesMax
62 0 : && lineGenerator.find("SHERPA", 0) == string::npos
63 0 : && lineGenerator.find("POWHEG-BOX", 0) == string::npos
64 0 : && lineGenerator.find("Pythia8", 0) == string::npos
65 0 : && lineGenerator.find("MadGraph", 0) == string::npos){
66 0 : iLine++;
67 0 : lineGenerator = " ";
68 0 : getline(infile,lineGenerator);
69 : }
70 0 : infile.close();
71 :
72 0 : vector <int> incom;
73 0 : vector <int> inter;
74 0 : 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 0 : 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 0 : string inParticleNamesSH[] = {
86 : // Leptons
87 0 : "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
88 : // Jet container
89 0 : "-93","93","-90","90","-91","91",
90 : // Quarks
91 0 : "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6"};
92 0 : string inParticleNamesMG[] = {
93 : // Leptons
94 0 : "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
95 : // Jet container
96 0 : "p~","p","l+","l-","vl~","vl",
97 : // Quarks
98 0 : "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};
99 :
100 : // Declare intermediate particle identifiers
101 0 : 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 0 : string interParticleNamesMG[] = {
114 : // Electroweak gauge bosons
115 0 : "a","z","w-","w+","h","W",
116 : // Top quarks
117 0 : "t~","t",
118 : // Dummy index as back-up
119 0 : "xx",
120 : // All squarks
121 0 : "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
122 0 : "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};
123 :
124 : // Declare final state particle identifiers
125 0 : 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 0 : string outParticleNamesMG[] = {
142 : // Leptons
143 0 : "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
144 : // Jet container and lepton containers
145 0 : "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
146 : // Quarks
147 0 : "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
148 : // SM uncoloured bosons
149 0 : "a","z","w-","w+","h","W",
150 : // Neutralino in SUSY
151 0 : "n1",
152 : // All squarks
153 0 : "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
154 0 : "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};
155 :
156 0 : string outParticleNamesSH[] = {
157 : // Leptons
158 0 : "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
159 : // Jet container and lepton containers
160 0 : "-93","93","-90","90","-91","91","0","0","0",
161 : // Quarks
162 0 : "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6",
163 : // SM uncoloured bosons
164 0 : "22","23","-24","24","25","0",
165 : // Neutralino in SUSY
166 0 : "1000022",
167 : // All squarks
168 0 : "-1000001","1000001","-1000002","1000002","-1000003","1000003",
169 0 : "-1000004","1000004",
170 0 : "-1000005","1000005","-1000006","1000006","-2000001","2000001",
171 0 : "-2000002","2000002",
172 0 : "-2000003","2000003","-2000004","2000004","-2000005","2000005",
173 0 : "-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 0 : int meGenType = (lineGenerator.find("MadGraph", 0) != string::npos) ? -1
183 0 : : (lineGenerator.find("SHERPA", 0) != string::npos) ? -2
184 0 : : (lineGenerator.find("POWHEG-BOX", 0) != string::npos) ? -3
185 0 : : (lineGenerator.find("Pythia8", 0) != string::npos) ? -4
186 : : 0;
187 :
188 0 : if (meGenType == -2){
189 : // Now read merging scale
190 : // Open path to LHEF and extract merging scale
191 0 : infile.open( (char*)( LHEpath +"_1.lhe").c_str());
192 0 : string lineTMS;
193 0 : while(lineTMS.find("NJetFinder ", 0) == string::npos){
194 0 : lineTMS = " ";
195 0 : getline(infile,lineTMS);
196 : }
197 0 : infile.close();
198 0 : lineTMS = lineTMS.substr(0,lineTMS.find(" 0.0 ",0));
199 0 : lineTMS = lineTMS.substr(lineTMS.find(" ", 0)+3,lineTMS.size());
200 : // Remove whitespaces
201 0 : while(lineTMS.find(" ", 0) != string::npos)
202 0 : lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
203 : // Replace d with e
204 0 : if ( lineTMS.find("d", 0) != string::npos)
205 0 : lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
206 0 : tms = atof((char*)lineTMS.c_str());
207 :
208 : // Now read hard process
209 : // Open path to LHEF and extract hard process
210 0 : infile.open( (char*)( LHEpath +"_0.lhe").c_str());
211 0 : string line;
212 0 : while(line.find("Process", 0) == string::npos){
213 0 : line = " ";
214 0 : getline(infile,line);
215 : }
216 0 : infile.close();
217 0 : line = line.substr(line.find(" ",0),line.size());
218 :
219 : // Cut string into incoming and outgoing pieces
220 0 : vector <string> pieces;
221 0 : pieces.push_back( line.substr(0,line.find("->", 0)) );
222 : // Do not count additional final jets
223 0 : int end = (line.find("{", 0) != string::npos) ? line.find("{", 0)-2
224 0 : : line.size();
225 0 : pieces.push_back( line.substr(line.find("->", 0)+2, end) );
226 :
227 : // Get incoming particles
228 0 : for(int i=0; i < nIn; ++i) {
229 0 : for(int n = pieces[0].find(inParticleNamesSH[i], 0);
230 0 : n != int(string::npos);
231 0 : n = pieces[0].find(inParticleNamesSH[i], n)) {
232 0 : incom.push_back(inParticleNumbers[i]);
233 0 : pieces[0].erase(pieces[0].begin()+n,
234 0 : 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 0 : inter.push_back(0);
241 : // Get final particles
242 0 : for(int i=0; i < nOut; ++i) {
243 0 : for(int n = pieces[1].find(outParticleNamesSH[i], 0);
244 0 : n != int(string::npos);
245 0 : n = pieces[1].find(outParticleNamesSH[i], n)) {
246 0 : outgo.push_back(outParticleNumbers[i]);
247 0 : pieces[1].erase(pieces[1].begin()+n,
248 0 : pieces[1].begin()+n+outParticleNamesSH[i].size());
249 : n=0;
250 : }
251 : }
252 :
253 0 : } else if (meGenType == -1 || meGenType == -3 || meGenType == -4){
254 :
255 : // Now read merging scale
256 0 : string lineTMS;
257 :
258 0 : if (meGenType == -1) {
259 : // Open path to LHEF and extract merging scale
260 0 : infile.open( (char*)( LHEpath +"_1.lhe").c_str());
261 0 : while(lineTMS.find("ktdurham", 0) == string::npos){
262 0 : lineTMS = " ";
263 0 : getline(infile,lineTMS);
264 : }
265 0 : lineTMS = lineTMS.substr(0,lineTMS.find("=",0));
266 0 : infile.close();
267 : } else {
268 0 : lineTMS = "30.";
269 : }
270 :
271 : // Remove whitespaces
272 0 : while(lineTMS.find(" ", 0) != string::npos)
273 0 : lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
274 : // Replace d with e
275 0 : if ( lineTMS.find("d", 0) != string::npos)
276 0 : lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
277 0 : tms = atof((char*)lineTMS.c_str());
278 :
279 : // Now read hard process
280 : // Open path to LHEF and extract hard process
281 0 : infile.open( (char*)( LHEpath +"_0.lhe").c_str());
282 0 : string line;
283 0 : while(line.find("@1", 0) == string::npos){
284 0 : line = " ";
285 0 : getline(infile,line);
286 : }
287 0 : infile.close();
288 0 : line = line.substr(0,line.find("@",0));
289 :
290 : // Count number of resonances
291 : int appearances = 0;
292 0 : for(int n = line.find("(", 0); n != int(string::npos);
293 0 : n = line.find("(", n)) {
294 0 : appearances++;
295 0 : n++;
296 : }
297 :
298 : // Cut string in incoming, resonance+decay and outgoing pieces
299 0 : vector <string> pieces;
300 0 : for(int i =0; i < appearances;++i) {
301 0 : int n = line.find("(", 0);
302 0 : pieces.push_back(line.substr(0,n));
303 0 : line = line.substr(n+1,line.size());
304 : }
305 : // Cut last resonance from rest
306 0 : if (appearances > 0) {
307 0 : pieces.push_back( line.substr(0,line.find(")",0)) );
308 0 : pieces.push_back( line.substr(line.find(")",0)+1,line.size()) );
309 0 : }
310 :
311 : // If the string was not cut into pieces, i.e. no resonance was
312 : // required, cut string using '>' as delimiter
313 0 : if (pieces.empty() ){
314 : appearances = 0;
315 0 : for(int n = line.find(">", 0); n != int(string::npos);
316 0 : n = line.find(">", n)) {
317 0 : appearances++;
318 0 : n++;
319 : }
320 :
321 : // Cut string in incoming and outgoing pieces
322 0 : for(int i =0; i < appearances;++i) {
323 0 : int n = line.find(">", 0);
324 0 : pieces.push_back(line.substr(0,n));
325 0 : line = line.substr(n+1,line.size());
326 : }
327 :
328 0 : if (appearances == 1) pieces.push_back(line);
329 0 : if (appearances > 1) {
330 0 : pieces.push_back( line.substr(0,line.find(">",0)) );
331 0 : pieces.push_back( line.substr(line.find(">",0)+1,line.size()) );
332 0 : }
333 : }
334 :
335 : // Get incoming particles
336 0 : for(int i=0; i < nIn; ++i) {
337 0 : for(int n = pieces[0].find(inParticleNamesMG[i], 0);
338 0 : n != int(string::npos);
339 0 : n = pieces[0].find(inParticleNamesMG[i], n)) {
340 0 : incom.push_back(inParticleNumbers[i]);
341 0 : pieces[0].erase(pieces[0].begin()+n,
342 0 : pieces[0].begin()+n+inParticleNamesMG[i].size());
343 : n=0;
344 : }
345 : }
346 :
347 : // Check intermediate resonances and decay products
348 0 : for(int i =1; i < int(pieces.size()); ++i){
349 : // Seperate strings into intermediate and outgoing, if not already done
350 0 : int k = pieces[i].find(">", 0);
351 :
352 0 : string intermediate = (pieces[i].find(">", 0) != string::npos) ?
353 0 : pieces[i].substr(0,k) : "";
354 0 : string outgoing = (pieces[i].find(">", 0) != string::npos) ?
355 0 : pieces[i].substr(k+1,pieces[i].size()) : pieces[i];
356 :
357 : // Get intermediate particles
358 0 : for(int j=0; j < nInt; ++j) {
359 0 : for(int n = intermediate.find(interParticleNamesMG[j], 0);
360 0 : n != int(string::npos);
361 0 : n = intermediate.find(interParticleNamesMG[j], n)) {
362 0 : inter.push_back(interParticleNumbers[j]);
363 0 : intermediate.erase(intermediate.begin()+n,
364 0 : intermediate.begin()+n+interParticleNamesMG[j].size());
365 : n=0;
366 : }
367 : }
368 :
369 : // Get outgoing particles
370 0 : for(int j=0; j < nOut; ++j) {
371 0 : for(int n = outgoing.find(outParticleNamesMG[j], 0);
372 0 : n != int(string::npos);
373 0 : n = outgoing.find(outParticleNamesMG[j], n)) {
374 0 : outgo.push_back(outParticleNumbers[j]);
375 0 : outgoing.erase(outgoing.begin()+n,
376 0 : 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 0 : if (inter.empty()) {
384 :
385 : // For final state bosons, bookkeep the final state boson as
386 : // intermediate as well.
387 : int nBosons = 0;
388 0 : for(int l=0; l < int(outgo.size()); ++l)
389 0 : if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
390 0 : nBosons++;
391 :
392 0 : int nZeros = (outgo.size() - nBosons)/2;
393 0 : for(int l=0; l < nZeros; ++l)
394 0 : inter.push_back(0);
395 0 : }
396 :
397 : // For final state bosons, bookkeep the final state boson as
398 : // intermediate as well.
399 0 : for(int l=0; l < int(outgo.size()); ++l)
400 0 : if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
401 0 : inter.push_back(outgo[l]);
402 :
403 0 : }
404 :
405 0 : } else {
406 :
407 0 : cout << "Reading of tms and hard process information from LHEF currently"
408 0 : << " only automated for MadEvent- or SHERPA-produced LHEF" << endl;
409 0 : int tempInt = 0;
410 0 : cout << "Use default process pp -> e+ve + jets? (0:no / 1:yes): ";
411 0 : cin >> tempInt;
412 0 : cout << endl;
413 :
414 0 : if (tempInt == 0){
415 0 : tempInt = 0;
416 0 : double tempDouble = 0.0;
417 0 : cout << "Please specify merging scale (kT Durham, in GeV): ";
418 0 : cin >> tempDouble;
419 0 : tms = tempDouble;
420 0 : cout << endl;
421 0 : cout << "Please specify first incoming particle ";
422 0 : cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
423 0 : cin >> tempInt;
424 0 : incom.push_back(tempInt);
425 0 : tempInt = 0;
426 0 : cout << endl;
427 0 : cout << "Please specify second incoming particle ";
428 0 : cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
429 0 : cin >> tempInt;
430 0 : incom.push_back(tempInt);
431 0 : tempInt = 0;
432 0 : cout << endl;
433 0 : cout << "Please specify intermediate particle, if any ";
434 0 : cout << "(0 for none, else PDG code): ";
435 0 : cin >> tempInt;
436 0 : inter.push_back(tempInt);
437 0 : cout << endl;
438 : do {
439 0 : tempInt = 0;
440 0 : cout << "Please specify outgoing particle ";
441 0 : cout << "(jet=2212, else PDG code, exit with 99): ";
442 0 : cin >> tempInt;
443 0 : if (tempInt != 99) outgo.push_back(tempInt);
444 0 : } while(tempInt != 99);
445 0 : cout << endl;
446 0 : } else {
447 0 : cout << "LHE file not produced by SHERPA or MG/ME - ";
448 0 : cout << "Using default process and tms" << endl;
449 0 : incom.push_back(2212);
450 0 : incom.push_back(2212);
451 0 : inter.push_back(24);
452 0 : outgo.push_back(-11);
453 0 : outgo.push_back(12);
454 0 : tms = 10.;
455 : }
456 0 : }
457 :
458 : // Now store incoming, intermediate and outgoing
459 : // Set intermediate tags
460 0 : for(int i=0; i < int(inter.size()); ++i)
461 0 : hardIntermediate.push_back(inter[i]);
462 :
463 : // Set the incoming particle tags
464 0 : if (incom.size() != 2)
465 0 : cout << "Only two incoming particles allowed" << endl;
466 : else {
467 0 : hardIncoming1 = incom[0];
468 0 : hardIncoming2 = incom[1];
469 : }
470 :
471 : // Remember final state bosons
472 : int nBosons = 0;
473 0 : for(int i=0; i < int(outgo.size()); ++i)
474 0 : if ( (abs(outgo[i]) > 20 && abs(outgo[i]) <= 25) || outgo[i] == 2400)
475 0 : nBosons++;
476 : // Remember b-quark container
477 : int nBQuarks = 0;
478 0 : for(int i=0; i < int(outgo.size()); ++i)
479 0 : if ( outgo[i] == 5000)
480 0 : nBQuarks++;
481 : // Remember jet container
482 : int nJets = 0;
483 0 : for(int i=0; i < int(outgo.size()); ++i)
484 0 : if ( outgo[i] == 2212)
485 0 : nJets++;
486 : // Remember lepton container
487 : int nLeptons = 0;
488 0 : for(int i=0; i < int(outgo.size()); ++i)
489 0 : if ( outgo[i] == 1100)
490 0 : nLeptons++;
491 : // Remember lepton container
492 : int nNeutrinos = 0;
493 0 : for(int i=0; i < int(outgo.size()); ++i)
494 0 : if ( outgo[i] == 1200)
495 0 : nNeutrinos++;
496 0 : int nContainers = nLeptons + nNeutrinos + nJets + nBQuarks;
497 :
498 : // Set final particle identifiers
499 0 : if ( (outgo.size() - nBosons - nContainers)%2 == 1) {
500 0 : cout << "Only even number of outgoing particles allowed" << endl;
501 0 : for(int i=0; i < int(outgo.size()); ++i)
502 0 : cout << outgo[i] << endl;
503 0 : } else {
504 :
505 : // Push back particles / antiparticles
506 0 : for(int i=0; i < int(outgo.size()); ++i)
507 0 : if (outgo[i] > 0
508 0 : && outgo[i] != 2212
509 0 : && outgo[i] != 5000
510 0 : && outgo[i] != 1100
511 0 : && outgo[i] != 1200
512 0 : && outgo[i] != 2400
513 0 : && outgo[i] != 1000022)
514 0 : hardOutgoing2.push_back( outgo[i]);
515 0 : else if (outgo[i] < 0)
516 0 : hardOutgoing1.push_back( outgo[i]);
517 :
518 : // Save final state W-boson container as particle
519 0 : for(int i=0; i < int(outgo.size()); ++i)
520 0 : if ( outgo[i] == 2400)
521 0 : hardOutgoing2.push_back( outgo[i]);
522 :
523 : // Push back jets, distribute evenly amongst particles / antiparticles
524 : // Push back majorana particles, distribute evenly
525 : int iNow = 0;
526 0 : for(int i=0; i < int(outgo.size()); ++i)
527 0 : if ( (outgo[i] == 2212
528 0 : || outgo[i] == 5000
529 0 : || outgo[i] == 1200
530 0 : || outgo[i] == 1000022)
531 0 : && iNow%2 == 0 ){
532 0 : hardOutgoing2.push_back( outgo[i]);
533 0 : iNow++;
534 0 : } else if ( (outgo[i] == 2212
535 0 : || outgo[i] == 5000
536 0 : || outgo[i] == 1100
537 0 : || outgo[i] == 1000022)
538 0 : && iNow%2 == 1 ){
539 0 : hardOutgoing1.push_back( outgo[i]);
540 0 : iNow++;
541 0 : }
542 : }
543 :
544 : // Done
545 0 : }
546 :
547 : //--------------------------------------------------------------------------
548 :
549 : // Function to translate a string specitying the core process into the
550 : // internal notation
551 : // Currently, the input string has to be in MadEvent notation
552 :
553 : void HardProcess::translateProcessString( string process){
554 :
555 0 : vector <int> incom;
556 0 : vector <int> inter;
557 0 : vector <int> outgo;
558 : // Particle identifiers, ordered in such a way that e.g. the "u"
559 : // in a mu is not mistaken for an u quark
560 0 : int inParticleNumbers[] = {
561 : // Leptons
562 : -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
563 : // Jet container
564 : 2212,2212,0,0,0,0,
565 : // Quarks
566 : -1,1,-2,2,-3,3,-4,4,-5,5,-6,6};
567 0 : string inParticleNamesMG[] = {
568 : // Leptons
569 0 : "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
570 : // Jet container
571 0 : "p~","p","l+","l-","vl~","vl",
572 : // Quarks
573 0 : "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};
574 :
575 : // Declare intermediate particle identifiers
576 0 : int interParticleNumbers[] = {
577 : // Electroweak gauge bosons
578 : 22,23,-24,24,25,2400,
579 : // All squarks
580 : -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
581 : -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
582 : -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006,
583 : // Top quarks
584 : -6,6,
585 : // Dummy index as back-up
586 : 0};
587 : // Declare names of intermediate particles
588 0 : string interParticleNamesMG[] = {
589 : // Electroweak gauge bosons
590 0 : "a","z","w-","w+","h","W",
591 : // All squarks
592 0 : "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
593 0 : "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2",
594 : // Top quarks
595 0 : "t~","t",
596 : // Dummy index as back-up
597 0 : "xx"};
598 :
599 : // Declare final state particle identifiers
600 0 : int outParticleNumbers[] = {
601 : // Leptons
602 : -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
603 : // Jet container and lepton containers
604 : 2212,2212,0,0,0,0,1200,1100,5000,
605 : // All squarks
606 : -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
607 : -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
608 : -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006,
609 : // Quarks
610 : -1,1,-2,2,-3,3,-4,4,-5,5,-6,6,
611 : // SM uncoloured bosons
612 : 22,23,-24,24,25,2400,
613 : // Neutralino in SUSY
614 : 1000022};
615 : // Declare names of final state particles
616 0 : string outParticleNamesMG[] = {
617 : // Leptons
618 0 : "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
619 : // Jet container and lepton containers
620 0 : "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
621 : // All squarks
622 0 : "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
623 0 : "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2",
624 : // Quarks
625 0 : "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
626 : // SM uncoloured bosons
627 0 : "a","z","w-","w+","h","W",
628 : // Neutralino in SUSY
629 0 : "n1"};
630 :
631 : // Declare size of particle name arrays
632 : int nIn = 30;
633 : int nInt = 33;
634 : int nOut = 64;
635 :
636 : // Start mapping user-defined particles onto particle ids.
637 : //string fullProc = "pp>{blaa,124}LEPTONS,NEUTRINOS";
638 0 : string fullProc = process;
639 :
640 : // Find user-defined hard process content
641 : // Count number of user particles
642 : int nUserParticles = 0;
643 0 : for(int n = fullProc.find("{", 0); n != int(string::npos);
644 0 : n = fullProc.find("{", n)) {
645 0 : nUserParticles++;
646 0 : n++;
647 : }
648 : // Cut user-defined particles from remaining process
649 0 : vector <string> userParticleStrings;
650 0 : for(int i =0; i < nUserParticles;++i) {
651 0 : int n = fullProc.find("{", 0);
652 0 : userParticleStrings.push_back(fullProc.substr(0,n));
653 0 : fullProc = fullProc.substr(n+1,fullProc.size());
654 : }
655 : // Cut remaining process string from rest
656 0 : if (nUserParticles > 0)
657 0 : userParticleStrings.push_back(
658 0 : fullProc.substr( 0, fullProc.find("}",0) ) );
659 : // Remove curly brackets and whitespace
660 0 : for(int i =0; i < int(userParticleStrings.size());++i) {
661 0 : while(userParticleStrings[i].find("{", 0) != string::npos)
662 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
663 0 : +userParticleStrings[i].find("{", 0));
664 0 : while(userParticleStrings[i].find("}", 0) != string::npos)
665 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
666 0 : +userParticleStrings[i].find("}", 0));
667 0 : while(userParticleStrings[i].find(" ", 0) != string::npos)
668 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
669 0 : +userParticleStrings[i].find(" ", 0));
670 : }
671 : // Convert particle numbers in user particle to integers
672 0 : vector<int>userParticleNumbers;
673 0 : if ( int(userParticleStrings.size()) > 1) {
674 0 : for( int i = 1; i < int(userParticleStrings.size()); ++i) {
675 0 : userParticleNumbers.push_back(
676 0 : atoi((char*)userParticleStrings[i].substr(
677 0 : userParticleStrings[i].find(",",0)+1,
678 0 : userParticleStrings[i].size()).c_str() ) );
679 : }
680 0 : }
681 : // Save remaining process string
682 0 : if (nUserParticles > 0)
683 0 : userParticleStrings.push_back(
684 0 : fullProc.substr(
685 0 : fullProc.find("}",0)+1, fullProc.size() ) );
686 : // Remove curly brackets and whitespace
687 0 : for( int i = 0; i < int(userParticleStrings.size()); ++i) {
688 0 : while(userParticleStrings[i].find("{", 0) != string::npos)
689 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
690 0 : +userParticleStrings[i].find("{", 0));
691 0 : while(userParticleStrings[i].find("}", 0) != string::npos)
692 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
693 0 : +userParticleStrings[i].find("}", 0));
694 0 : while(userParticleStrings[i].find(" ", 0) != string::npos)
695 0 : userParticleStrings[i].erase(userParticleStrings[i].begin()
696 0 : +userParticleStrings[i].find(" ", 0));
697 : }
698 :
699 : // Start mapping residual process string onto particle IDs.
700 : // Declare leftover process after user-defined particles have been converted
701 0 : string residualProc;
702 0 : if ( int(userParticleStrings.size()) > 1 )
703 0 : residualProc = userParticleStrings.front() + userParticleStrings.back();
704 : else
705 0 : residualProc = fullProc;
706 :
707 : // Remove comma separation
708 0 : while(residualProc.find(",", 0) != string::npos)
709 0 : residualProc.erase(residualProc.begin()+residualProc.find(",",0));
710 :
711 : // Count number of resonances
712 : int appearances = 0;
713 0 : for(int n = residualProc.find("(", 0); n != int(string::npos);
714 0 : n = residualProc.find("(", n)) {
715 0 : appearances++;
716 0 : n++;
717 : }
718 :
719 : // Cut string in incoming, resonance+decay and outgoing pieces
720 0 : vector <string> pieces;
721 0 : for(int i =0; i < appearances;++i) {
722 0 : int n = residualProc.find("(", 0);
723 0 : pieces.push_back(residualProc.substr(0,n));
724 0 : residualProc = residualProc.substr(n+1,residualProc.size());
725 : }
726 : // Cut last resonance from rest
727 0 : if (appearances > 0) {
728 0 : pieces.push_back( residualProc.substr(0,residualProc.find(")",0)) );
729 0 : pieces.push_back( residualProc.substr(
730 0 : residualProc.find(")",0)+1, residualProc.size()) );
731 0 : }
732 :
733 : // If the string was not cut into pieces, i.e. no resonance was
734 : // required, cut string using '>' as delimiter
735 0 : if (pieces.empty() ){
736 : appearances = 0;
737 0 : for(int n = residualProc.find(">", 0); n != int(string::npos);
738 0 : n = residualProc.find(">", n)) {
739 0 : appearances++;
740 0 : n++;
741 : }
742 :
743 : // Cut string in incoming and outgoing pieces
744 0 : for(int i =0; i < appearances;++i) {
745 0 : int n = residualProc.find(">", 0);
746 0 : pieces.push_back(residualProc.substr(0,n));
747 0 : residualProc = residualProc.substr(n+1,residualProc.size());
748 : }
749 :
750 0 : if (appearances == 1) pieces.push_back(residualProc);
751 0 : if (appearances > 1) {
752 0 : pieces.push_back( residualProc.substr(0,residualProc.find(">",0)) );
753 0 : pieces.push_back( residualProc.substr(
754 0 : residualProc.find(">",0)+1, residualProc.size()) );
755 0 : }
756 : }
757 :
758 : // Get incoming particles
759 0 : for(int i=0; i < nIn; ++i) {
760 0 : for(int n = pieces[0].find(inParticleNamesMG[i], 0);
761 0 : n != int(string::npos);
762 0 : n = pieces[0].find(inParticleNamesMG[i], n)) {
763 0 : incom.push_back(inParticleNumbers[i]);
764 0 : pieces[0].erase(pieces[0].begin()+n,
765 0 : pieces[0].begin()+n+inParticleNamesMG[i].size());
766 : n=0;
767 : }
768 : }
769 :
770 : // Check intermediate resonances and decay products
771 0 : for(int i =1; i < int(pieces.size()); ++i){
772 : // Seperate strings into intermediate and outgoing, if not already done
773 0 : int k = pieces[i].find(">", 0);
774 :
775 0 : string intermediate = (pieces[i].find(">", 0) != string::npos) ?
776 0 : pieces[i].substr(0,k) : "";
777 0 : string outgoing = (pieces[i].find(">", 0) != string::npos) ?
778 0 : pieces[i].substr(k+1,pieces[i].size()) : pieces[i];
779 :
780 : // Get intermediate particles
781 0 : for(int j=0; j < nInt; ++j) {
782 0 : for(int n = intermediate.find(interParticleNamesMG[j], 0);
783 0 : n != int(string::npos);
784 0 : n = intermediate.find(interParticleNamesMG[j], n)) {
785 0 : inter.push_back(interParticleNumbers[j]);
786 0 : intermediate.erase(intermediate.begin()+n,
787 0 : intermediate.begin()+n+interParticleNamesMG[j].size());
788 : n=0;
789 : }
790 : }
791 :
792 : // Get outgoing particles
793 0 : for(int j=0; j < nOut; ++j) {
794 0 : for(int n = outgoing.find(outParticleNamesMG[j], 0);
795 0 : n != int(string::npos);
796 0 : n = outgoing.find(outParticleNamesMG[j], n)) {
797 0 : outgo.push_back(outParticleNumbers[j]);
798 0 : outgoing.erase(outgoing.begin()+n,
799 0 : outgoing.begin()+n+outParticleNamesMG[j].size());
800 : n=0;
801 : }
802 : }
803 :
804 : // For arbitrary or non-existing intermediate, remember zero for each
805 : // two outgoing particles, without bosons.
806 0 : if (inter.empty()) {
807 :
808 : // For final state bosons, bookkeep the final state boson as
809 : // intermediate as well.
810 : int nBosons = 0;
811 0 : for(int l=0; l < int(outgo.size()); ++l)
812 0 : if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
813 0 : nBosons++;
814 :
815 0 : int nZeros = (outgo.size() - nBosons)/2;
816 0 : for(int l=0; l < nZeros; ++l)
817 0 : inter.push_back(0);
818 0 : }
819 :
820 : // For final state bosons, bookkeep the final state boson as
821 : // intermediate as well.
822 0 : for(int l=0; l < int(outgo.size()); ++l)
823 0 : if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
824 0 : inter.push_back(outgo[l]);
825 :
826 0 : }
827 :
828 : // Now store incoming, intermediate and outgoing
829 : // Set intermediate tags
830 0 : for(int i=0; i < int(inter.size()); ++i)
831 0 : hardIntermediate.push_back(inter[i]);
832 :
833 : // Set the incoming particle tags
834 0 : if (incom.size() != 2)
835 0 : cout << "Only two incoming particles allowed" << endl;
836 : else {
837 0 : hardIncoming1 = incom[0];
838 0 : hardIncoming2 = incom[1];
839 : }
840 :
841 : // Now store final particle identifiers
842 : // Start with user-defined particles.
843 0 : for( int i = 0; i < int(userParticleNumbers.size()); ++i)
844 0 : if (userParticleNumbers[i] > 0) {
845 0 : hardOutgoing2.push_back( userParticleNumbers[i]);
846 0 : hardIntermediate.push_back(0);
847 : // For non-existing intermediate, remember zero.
848 0 : } else if (userParticleNumbers[i] < 0) {
849 0 : hardOutgoing1.push_back( userParticleNumbers[i]);
850 : // For non-existing intermediate, remember zero.
851 0 : hardIntermediate.push_back(0);
852 : }
853 :
854 : // Push back particles / antiparticles
855 0 : for(int i=0; i < int(outgo.size()); ++i)
856 0 : if (outgo[i] > 0
857 0 : && outgo[i] != 2212
858 0 : && outgo[i] != 5000
859 0 : && outgo[i] != 1100
860 0 : && outgo[i] != 1200
861 0 : && outgo[i] != 2400
862 0 : && outgo[i] != 1000022)
863 0 : hardOutgoing2.push_back( outgo[i]);
864 0 : else if (outgo[i] < 0)
865 0 : hardOutgoing1.push_back( outgo[i]);
866 :
867 : // Save final state W-boson container as particle
868 0 : for(int i=0; i < int(outgo.size()); ++i)
869 0 : if ( outgo[i] == 2400)
870 0 : hardOutgoing2.push_back( outgo[i]);
871 :
872 : // Push back jets, distribute evenly among particles / antiparticles
873 : // Push back majorana particles, distribute evenly
874 : int iNow = 0;
875 0 : for(int i=0; i < int(outgo.size()); ++i)
876 0 : if ( (outgo[i] == 2212
877 0 : || outgo[i] == 5000
878 0 : || outgo[i] == 1200
879 0 : || outgo[i] == 1000022)
880 0 : && iNow%2 == 0 ){
881 0 : hardOutgoing2.push_back( outgo[i]);
882 0 : iNow++;
883 0 : } else if ( (outgo[i] == 2212
884 0 : || outgo[i] == 5000
885 0 : || outgo[i] == 1100
886 0 : || outgo[i] == 1000022)
887 0 : && iNow%2 == 1 ){
888 0 : hardOutgoing1.push_back( outgo[i]);
889 0 : iNow++;
890 0 : }
891 :
892 : // Done
893 0 : }
894 :
895 : //--------------------------------------------------------------------------
896 :
897 : // Function to check if the candidates stored in Pos1 and Pos2, together with
898 : // a proposed candidate iPos are allowed.
899 :
900 : bool HardProcess::allowCandidates(int iPos, vector<int> Pos1,
901 : vector<int> Pos2, const Event& event){
902 :
903 : bool allowed = true;
904 :
905 : // Find colour-partner of new candidate
906 0 : int type = (event[iPos].col() > 0) ? 1 : (event[iPos].acol() > 0) ? -1 : 0;
907 :
908 0 : if (type == 0) return true;
909 :
910 0 : if (type == 1){
911 0 : int col = event[iPos].col();
912 : int iPartner = 0;
913 0 : for(int i=0; i < int(event.size()); ++i)
914 0 : if ( i != iPos
915 0 : && (( event[i].isFinal() && event[i].acol() == col)
916 0 : ||( event[i].status() == -21 && event[i].col() == col) ))
917 0 : iPartner = i;
918 :
919 0 : vector<int> partners;
920 0 : for(int i=0; i < int(event.size()); ++i)
921 0 : for(int j=0; j < int(Pos1.size()); ++j)
922 0 : if ( Pos1[j] != 0 && i != Pos1[j] && event[Pos1[j]].colType() != 0
923 0 : && (( event[i].isFinal() && event[i].col() == event[Pos1[j]].acol())
924 0 : ||( event[i].status() == -21
925 0 : && event[i].acol() == event[Pos1[j]].acol()) ))
926 0 : partners.push_back(i);
927 :
928 : // Never allow equal initial partners!
929 0 : if (event[iPartner].status() == -21){
930 0 : for(int i=0; i < int(partners.size()); ++i)
931 0 : if ( partners[i] == iPartner)
932 0 : allowed = false;
933 0 : }
934 :
935 0 : } else {
936 0 : int col = event[iPos].acol();
937 : int iPartner = 0;
938 0 : for(int i=0; i < int(event.size()); ++i)
939 0 : if ( i != iPos
940 0 : && (( event[i].isFinal() && event[i].col() == col)
941 0 : ||(!event[i].isFinal() && event[i].acol() == col) ))
942 0 : iPartner = i;
943 :
944 0 : vector<int> partners;
945 0 : for(int i=0; i < int(event.size()); ++i)
946 0 : for(int j=0; j < int(Pos2.size()); ++j)
947 0 : if ( Pos2[j] != 0 && i != Pos2[j] && event[Pos2[j]].colType() != 0
948 0 : && (( event[i].isFinal() && event[i].acol() == event[Pos2[j]].col())
949 0 : ||( event[i].status() == -21
950 0 : && event[i].col() == event[Pos2[j]].col()) ))
951 0 : partners.push_back(i);
952 :
953 : // Never allow equal initial partners!
954 0 : if (event[iPartner].status() == -21){
955 0 : for(int i=0; i < int(partners.size()); ++i){
956 0 : if ( partners[i] == iPartner)
957 0 : allowed = false;
958 : }
959 0 : }
960 :
961 0 : }
962 :
963 0 : return allowed;
964 :
965 0 : }
966 :
967 : //--------------------------------------------------------------------------
968 :
969 : // Function to identify the hard subprocess in the current event
970 :
971 : void HardProcess::storeCandidates( const Event& event, string process){
972 :
973 : // Store the reference event
974 0 : state.clear();
975 0 : state = event;
976 :
977 : // Local copy of intermediate bosons
978 0 : vector<int> intermediates;
979 0 : for(int i =0; i < int(hardIntermediate.size());++i)
980 0 : intermediates.push_back( hardIntermediate[i]);
981 :
982 : // Local copy of outpoing partons
983 0 : vector<int> outgoing1;
984 0 : for(int i =0; i < int(hardOutgoing1.size());++i)
985 0 : outgoing1.push_back( hardOutgoing1[i]);
986 0 : vector<int> outgoing2;
987 0 : for(int i =0; i < int(hardOutgoing2.size());++i)
988 0 : outgoing2.push_back( hardOutgoing2[i]);
989 :
990 : // Clear positions of intermediate and outgoing particles
991 0 : PosIntermediate.resize(0);
992 0 : PosOutgoing1.resize(0);
993 0 : PosOutgoing2.resize(0);
994 0 : for(int i =0; i < int(hardIntermediate.size());++i)
995 0 : PosIntermediate.push_back(0);
996 0 : for(int i =0; i < int(hardOutgoing1.size());++i)
997 0 : PosOutgoing1.push_back(0);
998 0 : for(int i =0; i < int(hardOutgoing2.size());++i)
999 0 : PosOutgoing2.push_back(0);
1000 :
1001 : // For QCD dijet or e+e- > jets hard process, do not store any candidates,
1002 : // as to not discrimintate clusterings
1003 0 : if ( process.compare("pp>jj") == 0
1004 0 : || process.compare("e+e->jj") == 0
1005 0 : || process.compare("e+e->(z>jj)") == 0 ){
1006 0 : for(int i =0; i < int(hardOutgoing1.size());++i)
1007 0 : PosOutgoing1[i] = 0;
1008 0 : for(int i =0; i < int(hardOutgoing2.size());++i)
1009 0 : PosOutgoing2[i] = 0;
1010 : // Done
1011 0 : return;
1012 : }
1013 :
1014 : // Initialise vector of particles that were already identified as
1015 : // hard process particles
1016 0 : vector<int> iPosChecked;
1017 :
1018 : // If the hard process is specified only by containers, then add all
1019 : // particles matching with the containers to the hard process.
1020 : bool hasOnlyContainers = true;
1021 0 : for(int i =0; i < int(hardOutgoing1.size());++i)
1022 0 : if ( hardOutgoing1[i] != 1100
1023 0 : && hardOutgoing1[i] != 1200
1024 0 : && hardOutgoing1[i] != 5000)
1025 0 : hasOnlyContainers = false;
1026 0 : for(int i =0; i < int(hardOutgoing2.size());++i)
1027 0 : if ( hardOutgoing2[i] != 1100
1028 0 : && hardOutgoing2[i] != 1200
1029 0 : && hardOutgoing2[i] != 5000)
1030 0 : hasOnlyContainers = false;
1031 :
1032 0 : if (hasOnlyContainers){
1033 :
1034 0 : PosOutgoing1.resize(0);
1035 0 : PosOutgoing2.resize(0);
1036 :
1037 : // Try to find all unmatched hard process leptons.
1038 : // Loop through event to find outgoing lepton
1039 0 : for(int i=0; i < int(event.size()); ++i){
1040 :
1041 : // Skip non-final particles
1042 0 : if ( !event[i].isFinal() ) continue;
1043 :
1044 : // Skip all particles that have already been identified
1045 : bool skip = false;
1046 0 : for(int k=0; k < int(iPosChecked.size()); ++k){
1047 0 : if (i == iPosChecked[k])
1048 0 : skip = true;
1049 : }
1050 0 : if (skip) continue;
1051 :
1052 0 : for(int j=0; j < int(outgoing2.size()); ++j){
1053 :
1054 : // If the particle matches an outgoing neutrino, save it
1055 0 : if ( outgoing2[j] == 1100
1056 0 : && ( event[i].idAbs() == 11
1057 0 : || event[i].idAbs() == 13
1058 0 : || event[i].idAbs() == 15) ){
1059 0 : PosOutgoing2.push_back(i);
1060 0 : iPosChecked.push_back(i);
1061 : }
1062 :
1063 : // If the particle matches an outgoing lepton, save it
1064 0 : if ( outgoing2[j] == 1200
1065 0 : && ( event[i].idAbs() == 12
1066 0 : || event[i].idAbs() == 14
1067 0 : || event[i].idAbs() == 16) ){
1068 0 : PosOutgoing2.push_back(i);
1069 0 : iPosChecked.push_back(i);
1070 : }
1071 :
1072 : // If the particle matches an outgoing b-quark, save it
1073 0 : if ( outgoing2[j] == 5000 && event[i].idAbs() == 5 ){
1074 0 : PosOutgoing2.push_back(i);
1075 0 : iPosChecked.push_back(i);
1076 : }
1077 :
1078 : }
1079 :
1080 : // Skip all particles that have already been identified
1081 : skip = false;
1082 0 : for(int k=0; k < int(iPosChecked.size()); ++k){
1083 0 : if (i == iPosChecked[k])
1084 0 : skip = true;
1085 : }
1086 0 : if (skip) continue;
1087 :
1088 0 : for(int j=0; j < int(outgoing1.size()); ++j){
1089 :
1090 : // If the particle matches an outgoing neutrino, save it
1091 0 : if ( outgoing1[j] == 1100
1092 0 : && ( event[i].idAbs() == 11
1093 0 : || event[i].idAbs() == 13
1094 0 : || event[i].idAbs() == 15) ){
1095 0 : PosOutgoing1.push_back(i);
1096 0 : iPosChecked.push_back(i);
1097 : }
1098 :
1099 : // If the particle matches an outgoing lepton, save it
1100 0 : if ( outgoing1[j] == 1200
1101 0 : && ( event[i].idAbs() == 12
1102 0 : || event[i].idAbs() == 14
1103 0 : || event[i].idAbs() == 16) ){
1104 0 : PosOutgoing1.push_back(i);
1105 0 : iPosChecked.push_back(i);
1106 : }
1107 :
1108 : // If the particle matches an outgoing b-quark, save it
1109 0 : if ( outgoing1[j] == 5000 && event[i].idAbs() == 5 ){
1110 0 : PosOutgoing1.push_back(i);
1111 0 : iPosChecked.push_back(i);
1112 : }
1113 :
1114 : }
1115 0 : }
1116 :
1117 : // Done
1118 0 : return;
1119 : }
1120 :
1121 : // Now begin finding candidates when not only containers are used.
1122 :
1123 : // First try to find final state bosons
1124 0 : for(int i=0; i < int(intermediates.size()); ++i){
1125 :
1126 : // Do nothing if the intermediate boson is absent
1127 0 : if (intermediates[i] == 0) continue;
1128 :
1129 : // Do nothing if this boson does not match any final state boson
1130 : bool matchesFinalBoson = false;
1131 0 : for(int j =0; j< int(outgoing1.size()); ++j){
1132 0 : if ( intermediates[i] == outgoing1[j] )
1133 0 : matchesFinalBoson = true;
1134 : }
1135 0 : for(int j =0; j< int(outgoing2.size()); ++j){
1136 0 : if ( intermediates[i] == outgoing2[j] )
1137 0 : matchesFinalBoson = true;
1138 : }
1139 0 : if (!matchesFinalBoson) continue;
1140 :
1141 : // Loop through event
1142 0 : for(int j=0; j < int(event.size()); ++j) {
1143 :
1144 : // Skip all particles that have already been identified
1145 : bool skip = false;
1146 0 : for(int m=0; m < int(iPosChecked.size()); ++m)
1147 0 : if (j == iPosChecked[m]) skip = true;
1148 0 : if (skip) continue;
1149 :
1150 : // If the particle has a requested intermediate id, check if
1151 : // if is a final state boson
1152 0 : if ( (event[j].id() == intermediates[i])
1153 0 : ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
1154 :
1155 0 : PosIntermediate[i] = j;
1156 0 : intermediates[i] = 0;
1157 : // Be careful only to replace one index at a time!
1158 : bool indexSet = false;
1159 :
1160 0 : for(int k=0; k < int(outgoing1.size()); ++k) {
1161 0 : if (event[j].id() == outgoing1[k] && !indexSet){
1162 0 : PosOutgoing1[k] = j;
1163 0 : outgoing1[k] = 99;
1164 : indexSet = true;
1165 0 : }
1166 : }
1167 :
1168 0 : for(int k=0; k < int(outgoing2.size()); ++k) {
1169 0 : if (event[j].id() == outgoing2[k] && !indexSet){
1170 0 : PosOutgoing2[k] = j;
1171 0 : outgoing2[k] = 99;
1172 : indexSet = true;
1173 0 : }
1174 : }
1175 :
1176 : // Check for W-boson container
1177 0 : for(int k=0; k < int(outgoing2.size()); ++k) {
1178 0 : if (event[j].idAbs() == 24 && outgoing2[k] == 2400 && !indexSet ){
1179 0 : PosOutgoing2[k] = j;
1180 0 : outgoing2[k] = 99;
1181 : indexSet = true;
1182 0 : }
1183 : }
1184 :
1185 0 : iPosChecked.push_back(j);
1186 :
1187 0 : }
1188 0 : }
1189 0 : }
1190 :
1191 : // Second try to find particles coupled to intermediate bosons
1192 0 : for(int i=0; i < int(intermediates.size()); ++i){
1193 :
1194 : // Do nothing if the intermediate boson is absent
1195 0 : if (intermediates[i] == 0) continue;
1196 :
1197 : // Loop through event
1198 0 : for(int j=0; j < int(event.size()); ++j) {
1199 : // If the particle has a requested intermediate id, check if
1200 : // daughters are hard process particles
1201 0 : if ( (event[j].id() == intermediates[i])
1202 0 : ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
1203 : // If this particle is a potential intermediate
1204 0 : PosIntermediate[i] = j;
1205 0 : intermediates[i] = 0;
1206 : // If id's of daughters are good, store position
1207 0 : int iPos1 = event[j].daughter1();
1208 0 : int iPos2 = event[j].daughter2();
1209 :
1210 : // Loop through daughters to check if these contain some hard
1211 : // outgoing particles
1212 0 : for( int k=iPos1; k <= iPos2; ++k){
1213 0 : int id = event[k].id();
1214 :
1215 : // Skip all particles that have already been identified
1216 : bool skip = false;
1217 0 : for(int m=0; m < int(iPosChecked.size()); ++m)
1218 0 : if (k == iPosChecked[m]) skip = true;
1219 0 : if (skip) continue;
1220 :
1221 : // Check if daughter is hard outgoing particle
1222 0 : for(int l=0; l < int(outgoing2.size()); ++l)
1223 0 : if ( outgoing2[l] != 99 ){
1224 : // Found particle id
1225 0 : if (id == outgoing2[l]
1226 : // Found jet
1227 0 : || (id > 0 && abs(id) < 10 && outgoing2[l] == 2212) ){
1228 : // Store position
1229 0 : PosOutgoing2[l] = k;
1230 : // Remove the matched particle from the list
1231 0 : outgoing2[l] = 99;
1232 0 : iPosChecked.push_back(k);
1233 0 : break;
1234 : }
1235 :
1236 : }
1237 :
1238 : // Check if daughter is hard outgoing antiparticle
1239 0 : for(int l=0; l < int(outgoing1.size()); ++l)
1240 0 : if ( outgoing1[l] != 99 ){
1241 : // Found particle id
1242 0 : if (id == outgoing1[l]
1243 : // Found jet
1244 0 : || (id < 0 && abs(id) < 10 && outgoing1[l] == 2212) ){
1245 : // Store position
1246 0 : PosOutgoing1[l] = k;
1247 : // Remove the matched particle from the list
1248 0 : outgoing1[l] = 99;
1249 0 : iPosChecked.push_back(k);
1250 0 : break;
1251 : }
1252 :
1253 : }
1254 :
1255 0 : } // End loop through daughters
1256 0 : } // End if ids match
1257 : } // End loop through event
1258 0 : } // End loop though requested intermediates
1259 :
1260 : // If all outgoing particles were found, done
1261 : bool done = true;
1262 0 : for(int i=0; i < int(outgoing1.size()); ++i)
1263 0 : if (outgoing1[i] != 99)
1264 0 : done = false;
1265 0 : for(int i=0; i < int(outgoing2.size()); ++i)
1266 0 : if (outgoing2[i] != 99)
1267 0 : done = false;
1268 : // Return
1269 0 : if (done) return;
1270 :
1271 : // Leptons not associated with resonance are allowed.
1272 : // Try to find all unmatched hard process leptons.
1273 : // Loop through event to find outgoing lepton
1274 0 : for(int i=0; i < int(event.size()); ++i){
1275 : // Skip non-final particles and final partons
1276 0 : if ( !event[i].isFinal() || event[i].colType() != 0)
1277 : continue;
1278 : // Skip all particles that have already been identified
1279 : bool skip = false;
1280 0 : for(int k=0; k < int(iPosChecked.size()); ++k){
1281 0 : if (i == iPosChecked[k])
1282 0 : skip = true;
1283 : }
1284 0 : if (skip) continue;
1285 :
1286 : // Check if any hard outgoing leptons remain
1287 0 : for(int j=0; j < int(outgoing2.size()); ++j){
1288 : // Do nothing if this particle has already be found,
1289 : // or if this particle is a jet or quark
1290 0 : if ( outgoing2[j] == 99
1291 0 : || outgoing2[j] == 2212
1292 0 : || abs(outgoing2[j]) < 10)
1293 : continue;
1294 :
1295 : // If the particle matches an outgoing lepton, save it
1296 0 : if ( event[i].id() == outgoing2[j] ){
1297 0 : PosOutgoing2[j] = i;
1298 0 : outgoing2[j] = 99;
1299 0 : iPosChecked.push_back(i);
1300 : }
1301 : }
1302 :
1303 : // Check if any hard outgoing antileptons remain
1304 0 : for(int j=0; j < int(outgoing1.size()); ++j){
1305 : // Do nothing if this particle has already be found,
1306 : // or if this particle is a jet or quark
1307 0 : if ( outgoing1[j] == 99
1308 0 : || outgoing1[j] == 2212
1309 0 : || abs(outgoing1[j]) < 10)
1310 : continue;
1311 :
1312 : // If the particle matches an outgoing lepton, save it
1313 0 : if (event[i].id() == outgoing1[j] ){
1314 0 : PosOutgoing1[j] = i;
1315 0 : outgoing1[j] = 99;
1316 0 : iPosChecked.push_back(i);
1317 : }
1318 : }
1319 0 : }
1320 :
1321 0 : multimap<int,int> out2copy;
1322 0 : for(int i=0; i < int(event.size()); ++i)
1323 0 : for(int j=0; j < int(outgoing2.size()); ++j)
1324 : // Do nothing if this particle has already be found,
1325 : // or if this particle is a jet.
1326 0 : if ( outgoing2[j] != 99
1327 0 : && outgoing2[j] != 2212
1328 0 : && ( abs(outgoing2[j]) < 10
1329 0 : || (abs(outgoing2[j]) > 1000000 && abs(outgoing2[j]) < 1000010)
1330 0 : || (abs(outgoing2[j]) > 2000000 && abs(outgoing2[j]) < 2000010)
1331 0 : || abs(outgoing2[j]) == 1000021 )
1332 0 : && event[i].isFinal()
1333 0 : && event[i].id() == outgoing2[j] ){
1334 0 : out2copy.insert(make_pair(j, i));
1335 0 : }
1336 :
1337 0 : multimap<int,int> out1copy;
1338 0 : for(int i=0; i < int(event.size()); ++i)
1339 0 : for(int j=0; j < int(outgoing1.size()); ++j)
1340 : // Do nothing if this particle has already be found,
1341 : // or if this particle is a jet.
1342 0 : if ( outgoing1[j] != 99
1343 0 : && outgoing1[j] != 2212
1344 0 : && ( abs(outgoing1[j]) < 10
1345 0 : || (abs(outgoing1[j]) > 1000000 && abs(outgoing1[j]) < 1000010)
1346 0 : || (abs(outgoing1[j]) > 2000000 && abs(outgoing1[j]) < 2000010)
1347 0 : || abs(outgoing2[j]) == 1000021 )
1348 0 : && event[i].isFinal()
1349 0 : && event[i].id() == outgoing1[j] ){
1350 0 : out1copy.insert(make_pair(j, i));
1351 0 : }
1352 :
1353 0 : if ( out1copy.size() > out2copy.size()){
1354 :
1355 : // In case the index of the multimap is filled twice, make sure not to
1356 : // arbitrarily overwrite set values.
1357 0 : vector<int> indexWasSet;
1358 0 : for ( multimap<int, int>::iterator it = out2copy.begin();
1359 0 : it != out2copy.end(); ++it ) {
1360 0 : if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1361 :
1362 : // Skip all particles that have already been identified
1363 : bool skip = false;
1364 0 : for(int k=0; k < int(iPosChecked.size()); ++k)
1365 0 : if (it->second == iPosChecked[k]) skip = true;
1366 : // Skip all indices that have already been identified
1367 0 : for(int k=0; k < int(indexWasSet.size()); ++k)
1368 0 : if (it->first == indexWasSet[k]) skip = true;
1369 0 : if (skip) continue;
1370 :
1371 : // Save parton
1372 0 : PosOutgoing2[it->first] = it->second;
1373 : // remove entry form lists
1374 0 : outgoing2[it->first] = 99;
1375 0 : iPosChecked.push_back(it->second);
1376 0 : indexWasSet.push_back(it->first);
1377 0 : }
1378 : }
1379 :
1380 0 : indexWasSet.resize(0);
1381 0 : for ( multimap<int, int>::iterator it = out1copy.begin();
1382 0 : it != out1copy.end(); ++it ) {
1383 0 : if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1384 :
1385 : // Skip all particles that have already been identified
1386 : bool skip = false;
1387 0 : for(int k=0; k < int(iPosChecked.size()); ++k)
1388 0 : if (it->second == iPosChecked[k]) skip = true;
1389 : // Skip all indices that have already been identified
1390 0 : for(int k=0; k < int(indexWasSet.size()); ++k)
1391 0 : if (it->first == indexWasSet[k]) skip = true;
1392 0 : if (skip) continue;
1393 :
1394 : // Save parton
1395 0 : PosOutgoing1[it->first] = it->second;
1396 : // remove entry form lists
1397 0 : outgoing1[it->first] = 99;
1398 0 : iPosChecked.push_back(it->second);
1399 0 : indexWasSet.push_back(it->first);
1400 0 : }
1401 : }
1402 :
1403 0 : } else {
1404 :
1405 : // In case the index of the multimap is filled twice, make sure not to
1406 : // arbitraryly overwrite set values.
1407 0 : vector<int> indexWasSet;
1408 0 : for ( multimap<int, int>::iterator it = out1copy.begin();
1409 0 : it != out1copy.end(); ++it ) {
1410 0 : if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1411 :
1412 : // Skip all particles that have already been identified
1413 : bool skip = false;
1414 0 : for(int k=0; k < int(iPosChecked.size()); ++k)
1415 0 : if (it->second == iPosChecked[k]) skip = true;
1416 : // Skip all indices that have already been identified
1417 0 : for(int k=0; k < int(indexWasSet.size()); ++k)
1418 0 : if (it->first == indexWasSet[k]) skip = true;
1419 0 : if (skip) continue;
1420 :
1421 : // Save parton
1422 0 : PosOutgoing1[it->first] = it->second;
1423 : // remove entry form lists
1424 0 : outgoing1[it->first] = 99;
1425 0 : iPosChecked.push_back(it->second);
1426 0 : indexWasSet.push_back(it->first);
1427 0 : }
1428 : }
1429 :
1430 0 : indexWasSet.resize(0);
1431 0 : for ( multimap<int, int>::iterator it = out2copy.begin();
1432 0 : it != out2copy.end(); ++it ) {
1433 0 : if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){
1434 :
1435 : // Skip all particles that have already been identified
1436 : bool skip = false;
1437 0 : for(int k=0; k < int(iPosChecked.size()); ++k)
1438 0 : if (it->second == iPosChecked[k]) skip = true;
1439 : // Skip all indices that have already been identified
1440 0 : for(int k=0; k < int(indexWasSet.size()); ++k)
1441 0 : if (it->first == indexWasSet[k]) skip = true;
1442 0 : if (skip) continue;
1443 :
1444 : // Save parton
1445 0 : PosOutgoing2[it->first] = it->second;
1446 : // remove entry form lists
1447 0 : outgoing2[it->first] = 99;
1448 0 : iPosChecked.push_back(it->second);
1449 0 : indexWasSet.push_back(it->first);
1450 0 : }
1451 : }
1452 0 : }
1453 :
1454 : // It sometimes happens that MadEvent does not put a
1455 : // heavy coloured resonance into the LHE file, even if requested.
1456 : // This means that the decay products of this resonance need to be
1457 : // found separately.
1458 : // Loop through event to find hard process (anti)quarks
1459 0 : for(int i=0; i < int(event.size()); ++i){
1460 :
1461 : // Skip non-final particles and final partons
1462 0 : if ( !event[i].isFinal() || event[i].colType() == 0)
1463 : continue;
1464 :
1465 : // Skip all particles that have already been identified
1466 : bool skip = false;
1467 0 : for(int k=0; k < int(iPosChecked.size()); ++k){
1468 0 : if (i == iPosChecked[k])
1469 0 : skip = true;
1470 : }
1471 0 : if (skip) continue;
1472 :
1473 : // Check if any hard outgoing quarks remain
1474 0 : for(int j=0; j < int(outgoing2.size()); ++j){
1475 : // Do nothing if this particle has already be found,
1476 : // or if this particle is a jet, lepton container or lepton
1477 0 : if ( outgoing2[j] == 99
1478 0 : || outgoing2[j] == 2212
1479 0 : || abs(outgoing2[j]) > 10)
1480 : continue;
1481 : // If the particle matches an outgoing quark, save it
1482 0 : if (event[i].id() == outgoing2[j]){
1483 : // Save parton
1484 0 : PosOutgoing2[j] = i;
1485 : // remove entry form lists
1486 0 : outgoing2[j] = 99;
1487 0 : iPosChecked.push_back(i);
1488 : }
1489 : }
1490 :
1491 : // Check if any hard outgoing antiquarks remain
1492 0 : for(int j=0; j < int(outgoing1.size()); ++j){
1493 : // Do nothing if this particle has already be found,
1494 : // or if this particle is a jet, lepton container or lepton
1495 0 : if ( outgoing1[j] == 99
1496 0 : || outgoing1[j] == 2212
1497 0 : || abs(outgoing1[j]) > 10)
1498 : continue;
1499 : // If the particle matches an outgoing antiquark, save it
1500 0 : if (event[i].id() == outgoing1[j]){
1501 : // Save parton
1502 0 : PosOutgoing1[j] = i;
1503 : // Remove parton from list
1504 0 : outgoing1[j] = 99;
1505 0 : iPosChecked.push_back(i);
1506 : }
1507 : }
1508 0 : }
1509 :
1510 : // Done
1511 0 : }
1512 :
1513 : //--------------------------------------------------------------------------
1514 :
1515 : // Function to check if the particle event[iPos] matches any of
1516 : // the stored outgoing particles of the hard subprocess
1517 :
1518 : bool HardProcess::matchesAnyOutgoing(int iPos, const Event& event){
1519 :
1520 : // Match quantum numbers of any first outgoing candidate
1521 : bool matchQN1 = false;
1522 : // Match quantum numbers of any second outgoing candidate
1523 : bool matchQN2 = false;
1524 : // Match parton in the hard process,
1525 : // or parton from decay of electroweak boson in hard process,
1526 : // or parton from decay of electroweak boson from decay of top
1527 : bool matchHP = false;
1528 :
1529 : // Check outgoing candidates
1530 0 : for(int i=0; i < int(PosOutgoing1.size()); ++i)
1531 : // Compare particle properties
1532 0 : if ( event[iPos].id() == state[PosOutgoing1[i]].id()
1533 0 : && event[iPos].colType() == state[PosOutgoing1[i]].colType()
1534 0 : && event[iPos].chargeType() == state[PosOutgoing1[i]].chargeType()
1535 0 : && ( ( event[iPos].col() > 0
1536 0 : && event[iPos].col() == state[PosOutgoing1[i]].col())
1537 0 : || ( event[iPos].acol() > 0
1538 0 : && event[iPos].acol() == state[PosOutgoing1[i]].acol()))
1539 0 : && event[iPos].charge() == state[PosOutgoing1[i]].charge() )
1540 0 : matchQN1 = true;
1541 :
1542 : // Check outgoing candidates
1543 0 : for(int i=0; i < int(PosOutgoing2.size()); ++i)
1544 : // Compare particle properties
1545 0 : if ( event[iPos].id() == state[PosOutgoing2[i]].id()
1546 0 : && event[iPos].colType() == state[PosOutgoing2[i]].colType()
1547 0 : && event[iPos].chargeType() == state[PosOutgoing2[i]].chargeType()
1548 0 : && ( ( event[iPos].col() > 0
1549 0 : && event[iPos].col() == state[PosOutgoing2[i]].col())
1550 0 : || ( event[iPos].acol() > 0
1551 0 : && event[iPos].acol() == state[PosOutgoing2[i]].acol()))
1552 0 : && event[iPos].charge() == state[PosOutgoing2[i]].charge() )
1553 0 : matchQN2 = true;
1554 :
1555 : // Check if maps to hard process:
1556 : // Check that particle is in hard process
1557 0 : if ( event[iPos].mother1()*event[iPos].mother2() == 12
1558 : // Or particle has taken recoil from first splitting
1559 0 : || ( event[iPos].status() == 44
1560 0 : && event[event[iPos].mother1()].mother1()
1561 0 : *event[event[iPos].mother1()].mother2() == 12 )
1562 : // Or particle has on-shell resonace as mother
1563 0 : || ( event[iPos].status() == 23
1564 0 : && event[event[iPos].mother1()].mother1()
1565 0 : *event[event[iPos].mother1()].mother2() == 12 )
1566 : // Or particle has on-shell resonace as mother,
1567 : // which again has and on-shell resonance as mother
1568 0 : || ( event[iPos].status() == 23
1569 0 : && event[event[iPos].mother1()].status() == -22
1570 0 : && event[event[event[iPos].mother1()].mother1()].status() == -22
1571 0 : && event[event[event[iPos].mother1()].mother1()].mother1()
1572 0 : *event[event[event[iPos].mother1()].mother1()].mother2() == 12 ) )
1573 0 : matchHP = true;
1574 :
1575 : // Done
1576 0 : return ( matchHP && (matchQN1 || matchQN2) );
1577 :
1578 : }
1579 :
1580 :
1581 : //--------------------------------------------------------------------------
1582 :
1583 : // Function to check if instead of the particle event[iCandidate], another
1584 : // particle could serve as part of the hard process. Assumes that iCandidate
1585 : // is already stored as part of the hard process.
1586 :
1587 : bool HardProcess::findOtherCandidates(int iPos, const Event& event,
1588 : bool doReplace){
1589 :
1590 : // Return value
1591 : bool foundCopy = false;
1592 :
1593 : // Save stored candidates' properties.
1594 0 : int id = event[iPos].id();
1595 0 : int col = event[iPos].col();
1596 0 : int acl = event[iPos].acol();
1597 :
1598 : // If the particle's mother is an identified intermediate resonance,
1599 : // then do not attempt any replacement.
1600 0 : int iMoth1 = event[iPos].mother1();
1601 0 : int iMoth2 = event[iPos].mother2();
1602 0 : if ( iMoth1 > 0 && iMoth2 == 0 ) {
1603 : bool hasIdentifiedMother = false;
1604 0 : for(int i=0; i < int(PosIntermediate.size()); ++i)
1605 : // Compare particle properties
1606 0 : if ( event[iMoth1].id() == state[PosIntermediate[i]].id()
1607 0 : && event[iMoth1].colType() == state[PosIntermediate[i]].colType()
1608 0 : && event[iMoth1].chargeType() == state[PosIntermediate[i]].chargeType()
1609 0 : && event[iMoth1].col() == state[PosIntermediate[i]].col()
1610 0 : && event[iMoth1].acol() == state[PosIntermediate[i]].acol()
1611 0 : && event[iMoth1].charge() == state[PosIntermediate[i]].charge() )
1612 0 : hasIdentifiedMother = true;
1613 0 : if(hasIdentifiedMother && event[iMoth1].id() != id) return false;
1614 0 : }
1615 :
1616 : // Find candidate amongst the already stored ME process candidates.
1617 0 : vector<int> candidates1;
1618 0 : vector<int> candidates2;
1619 : // Check outgoing candidates
1620 0 : for(int i=0; i < int(PosOutgoing1.size()); ++i)
1621 : // Compare particle properties
1622 0 : if ( id == state[PosOutgoing1[i]].id()
1623 0 : && col == state[PosOutgoing1[i]].col()
1624 0 : && acl == state[PosOutgoing1[i]].acol() )
1625 0 : candidates1.push_back(i);
1626 : // Check outgoing candidates
1627 0 : for(int i=0; i < int(PosOutgoing2.size()); ++i)
1628 : // Compare particle properties
1629 0 : if ( id == state[PosOutgoing2[i]].id()
1630 0 : && col == state[PosOutgoing2[i]].col()
1631 0 : && acl == state[PosOutgoing2[i]].acol() )
1632 0 : candidates2.push_back(i);
1633 :
1634 : // If more / less than one stored candidate for iPos has been found, exit.
1635 0 : if ( candidates1.size() + candidates2.size() != 1 ) return false;
1636 :
1637 : // Now check for other allowed candidates.
1638 0 : map<int,int> further1;
1639 0 : for(int i=0; i < int(state.size()); ++i)
1640 0 : for(int j=0; j < int(PosOutgoing1.size()); ++j)
1641 : // Do nothing if this particle has already be found,
1642 : // or if this particle is a jet, lepton container or lepton
1643 0 : if ( state[i].isFinal()
1644 0 : && i != PosOutgoing1[j]
1645 0 : && state[PosOutgoing1[j]].id() == id
1646 0 : && state[i].id() == id ){
1647 : // Declare vector of already existing candiates.
1648 0 : vector<int> newPosOutgoing1;
1649 0 : for(int k=0; k < int(PosOutgoing1.size()); ++k)
1650 0 : if ( k != j ) newPosOutgoing1.push_back( PosOutgoing1[k] );
1651 : // If allowed, remember replacment parton.
1652 0 : if ( allowCandidates(i, newPosOutgoing1, PosOutgoing2, state) )
1653 0 : further1.insert(make_pair(j, i));
1654 0 : }
1655 :
1656 : // Now check for other allowed candidates.
1657 0 : map<int,int> further2;
1658 0 : for(int i=0; i < int(state.size()); ++i)
1659 0 : for(int j=0; j < int(PosOutgoing2.size()); ++j)
1660 : // Do nothing if this particle has already be found,
1661 : // or if this particle is a jet, lepton container or lepton
1662 0 : if ( state[i].isFinal()
1663 0 : && i != PosOutgoing2[j]
1664 0 : && state[PosOutgoing2[j]].id() == id
1665 0 : && state[i].id() == id ){
1666 : // Declare vector of already existing candidates.
1667 0 : vector<int> newPosOutgoing2;
1668 0 : for(int k=0; k < int(PosOutgoing2.size()); ++k)
1669 0 : if ( k != j ) newPosOutgoing2.push_back( PosOutgoing2[k] );
1670 : // If allowed, remember replacment parton.
1671 0 : if ( allowCandidates(i, PosOutgoing1, newPosOutgoing2, state) )
1672 0 : further2.insert(make_pair(j, i));
1673 0 : }
1674 :
1675 : // Remove all hard process particles that would be counted twice.
1676 0 : map<int,int>::iterator it2 = further2.begin();
1677 0 : while(it2 != further2.end()) {
1678 : bool remove = false;
1679 0 : for(int j=0; j < int(PosOutgoing2.size()); ++j)
1680 0 : if (it2->second == PosOutgoing2[j] ) remove = true;
1681 0 : if ( remove ) further2.erase(it2++);
1682 0 : else ++it2;
1683 : }
1684 0 : map<int,int>::iterator it1 = further1.begin();
1685 0 : while(it1 != further1.end()) {
1686 : bool remove = false;
1687 0 : for(int j=0; j < int(PosOutgoing1.size()); ++j)
1688 0 : if (it1->second == PosOutgoing1[j] ) remove = true;
1689 0 : if ( remove ) further1.erase(it1++);
1690 0 : else ++it1;
1691 : }
1692 :
1693 : // Decide of a replacment candidate has been found.
1694 0 : foundCopy = (doReplace)
1695 0 : ? exchangeCandidates(candidates1, candidates2, further1, further2)
1696 0 : : (further1.size() + further2.size() > 0);
1697 :
1698 : // Done
1699 0 : return foundCopy;
1700 :
1701 0 : }
1702 :
1703 : //--------------------------------------------------------------------------
1704 :
1705 : // Function to exchange hard process candidates.
1706 :
1707 : bool HardProcess::exchangeCandidates( vector<int> candidates1,
1708 : vector<int> candidates2, map<int,int> further1, map<int,int> further2) {
1709 :
1710 0 : int nOld1 = candidates1.size();
1711 0 : int nOld2 = candidates2.size();
1712 0 : int nNew1 = further1.size();
1713 0 : int nNew2 = further2.size();
1714 : bool exchanged = false;
1715 : // Replace, if one-to-one correspondence exists.
1716 0 : if ( nOld1 == 1 && nOld2 == 0 && nNew1 == 1 && nNew2 == 0){
1717 0 : PosOutgoing1[further1.begin()->first] = further1.begin()->second;
1718 : exchanged = true;
1719 0 : } else if ( nOld1 == 0 && nOld2 == 1 && nNew1 == 0 && nNew2 == 1){
1720 0 : PosOutgoing2[further2.begin()->first] = further2.begin()->second;
1721 : exchanged = true;
1722 : // Else simply swap with the first candidate.
1723 0 : } else if ( nNew1 > 1 && nNew2 == 0 ) {
1724 0 : PosOutgoing1[further1.begin()->first] = further1.begin()->second;
1725 : exchanged = true;
1726 0 : } else if ( nNew1 == 0 && nNew2 > 0 ) {
1727 0 : PosOutgoing2[further2.begin()->first] = further2.begin()->second;
1728 : exchanged = true;
1729 0 : }
1730 :
1731 : // Done
1732 0 : return exchanged;
1733 :
1734 : }
1735 :
1736 : //--------------------------------------------------------------------------
1737 :
1738 : // Function to get the number of coloured final state partons in the
1739 : // hard process
1740 :
1741 : int HardProcess::nQuarksOut(){
1742 : int nFin =0;
1743 0 : for(int i =0; i< int(hardOutgoing1.size()); ++i){
1744 0 : if (hardOutgoing1[i] == 2212 || abs(hardOutgoing1[i]) < 10) nFin++;
1745 : }
1746 0 : for(int i =0; i< int(hardOutgoing2.size()); ++i){
1747 0 : if (hardOutgoing2[i] == 2212 || abs(hardOutgoing2[i]) < 10) nFin++;
1748 : }
1749 : // For very loose hard process definition, check number of hard process
1750 : // b-quarks explicitly.
1751 0 : for(int i =0; i< int(hardOutgoing1.size()); ++i)
1752 0 : if (hardOutgoing1[i] == 5000)
1753 0 : for(int j =0; j< int(PosOutgoing1.size()); ++j)
1754 0 : if (state[PosOutgoing1[j]].idAbs() == 5)
1755 0 : nFin++;
1756 0 : for(int i =0; i< int(hardOutgoing2.size()); ++i)
1757 0 : if (hardOutgoing2[i] == 5000)
1758 0 : for(int j =0; j< int(PosOutgoing2.size()); ++j)
1759 0 : if (state[PosOutgoing2[j]].idAbs() == 5)
1760 0 : nFin++;
1761 0 : return nFin;
1762 : }
1763 :
1764 : //--------------------------------------------------------------------------
1765 :
1766 : // Function to get the number of uncoloured final state particles in the
1767 : // hard process
1768 :
1769 : int HardProcess::nLeptonOut(){
1770 : int nFin =0;
1771 0 : for(int i =0; i< int(hardOutgoing1.size()); ++i){
1772 0 : if (abs(hardOutgoing1[i]) > 10 && abs(hardOutgoing1[i]) < 20) nFin++;
1773 : // Bookkeep MSSM neutralinos as leptons
1774 0 : if (abs(hardOutgoing1[i]) == 1000022) nFin++;
1775 : // Bookkeep sleptons as leptons
1776 0 : if ( abs(hardOutgoing1[i]) == 1000011 || abs(hardOutgoing1[i]) == 2000011
1777 0 : || abs(hardOutgoing1[i]) == 1000013 || abs(hardOutgoing1[i]) == 2000013
1778 0 : || abs(hardOutgoing1[i]) == 1000015 || abs(hardOutgoing1[i]) == 2000015)
1779 0 : nFin++;
1780 : }
1781 0 : for(int i =0; i< int(hardOutgoing2.size()); ++i){
1782 0 : if (abs(hardOutgoing2[i]) > 10 && abs(hardOutgoing2[i]) < 20) nFin++;
1783 : // Bookkeep MSSM neutralinos as leptons
1784 0 : if (abs(hardOutgoing2[i]) == 1000022) nFin++;
1785 : // Bookkeep sleptons as leptons
1786 0 : if ( abs(hardOutgoing2[i]) == 1000011 || abs(hardOutgoing2[i]) == 2000011
1787 0 : || abs(hardOutgoing2[i]) == 1000013 || abs(hardOutgoing2[i]) == 2000013
1788 0 : || abs(hardOutgoing2[i]) == 1000015 || abs(hardOutgoing2[i]) == 2000015)
1789 0 : nFin++;
1790 : }
1791 :
1792 : // For very loose hard process definition, check number of hard process
1793 : // lepton explicitly.
1794 : // Check lepton / neutrino containers as leptons
1795 0 : for(int i =0; i< int(hardOutgoing1.size()); ++i)
1796 0 : if (hardOutgoing1[i] == 1100)
1797 0 : for(int j =0; j< int(PosOutgoing1.size()); ++j)
1798 0 : if ( abs(state[PosOutgoing1[j]].id()) == 11
1799 0 : || abs(state[PosOutgoing1[j]].id()) == 13
1800 0 : || abs(state[PosOutgoing1[j]].id()) == 15 )
1801 0 : nFin++;
1802 0 : for(int i =0; i< int(hardOutgoing2.size()); ++i)
1803 0 : if (hardOutgoing2[i] == 1200)
1804 0 : for(int j =0; j< int(PosOutgoing2.size()); ++j)
1805 0 : if ( abs(state[PosOutgoing2[j]].id()) == 12
1806 0 : || abs(state[PosOutgoing2[j]].id()) == 14
1807 0 : || abs(state[PosOutgoing2[j]].id()) == 16 )
1808 0 : nFin++;
1809 0 : return nFin;
1810 : }
1811 :
1812 : //--------------------------------------------------------------------------
1813 :
1814 : // Function to get the number of uncoloured final state particles in the
1815 : // hard process
1816 :
1817 : int HardProcess::nBosonsOut(){
1818 : int nFin =0;
1819 0 : for(int i =0; i< int(hardOutgoing1.size()); ++i){
1820 0 : if (abs(hardOutgoing1[i]) > 20 && abs(hardOutgoing1[i]) <= 25) nFin++;
1821 : }
1822 0 : for(int i =0; i< int(hardOutgoing2.size()); ++i){
1823 0 : if (abs(hardOutgoing2[i]) > 20 && abs(hardOutgoing2[i]) <= 25) nFin++;
1824 0 : if ( hardOutgoing2[i] == 2400) nFin++;
1825 : }
1826 0 : return nFin;
1827 : }
1828 :
1829 : //--------------------------------------------------------------------------
1830 :
1831 : // Function to get the number of coloured initial state partons in the
1832 : // hard process
1833 :
1834 : int HardProcess::nQuarksIn(){
1835 : int nIn =0;
1836 0 : if (hardIncoming1 == 2212 || abs(hardIncoming1) < 10) nIn++;
1837 0 : if (hardIncoming2 == 2212 || abs(hardIncoming2) < 10) nIn++;
1838 0 : return nIn;
1839 : }
1840 :
1841 : //--------------------------------------------------------------------------
1842 :
1843 : // Function to get the number of uncoloured initial state particles in the
1844 : // hard process
1845 :
1846 : int HardProcess::nLeptonIn(){
1847 : int nIn =0;
1848 0 : if (abs(hardIncoming1) > 10 && abs(hardIncoming1) < 20) nIn++;
1849 0 : if (abs(hardIncoming2) > 10 && abs(hardIncoming2) < 20) nIn++;
1850 0 : return nIn;
1851 : }
1852 :
1853 : //--------------------------------------------------------------------------
1854 :
1855 : // Function to report if a resonace decay was found in the
1856 : // 2->2 hard sub-process in the current state
1857 :
1858 : int HardProcess::hasResInCurrent(){
1859 0 : for(int i =0; i< int(PosIntermediate.size()); ++i)
1860 0 : if (PosIntermediate[i] == 0) return 0;
1861 : // Do not count final state bosons as resonaces
1862 0 : for(int i =0; i< int(PosIntermediate.size()); ++i){
1863 0 : for(int j =0; j< int(PosOutgoing1.size()); ++j){
1864 0 : if ( PosIntermediate[i] == PosOutgoing1[j] )
1865 0 : return 0;
1866 : }
1867 0 : for(int j =0; j< int(PosOutgoing2.size()); ++j){
1868 0 : if ( PosIntermediate[i] == PosOutgoing2[j] )
1869 0 : return 0;
1870 : }
1871 : }
1872 0 : return 1;
1873 0 : }
1874 :
1875 : //--------------------------------------------------------------------------
1876 :
1877 : // Function to report the number of resonace decays in the 2->2 sub-process
1878 : // of the current state
1879 :
1880 : int HardProcess::nResInCurrent(){
1881 : int nRes = 0;
1882 0 : for(int i =0; i< int(PosIntermediate.size()); ++i){
1883 0 : if (PosIntermediate[i] != 0) {
1884 : bool matchesFinalBoson = false;
1885 0 : for(int j =0; j< int(PosOutgoing1.size()); ++j){
1886 0 : if ( PosIntermediate[i] == PosOutgoing1[j] )
1887 0 : matchesFinalBoson = true;
1888 : }
1889 0 : for(int j =0; j< int(PosOutgoing2.size()); ++j){
1890 0 : if ( PosIntermediate[i] == PosOutgoing2[j] )
1891 0 : matchesFinalBoson = true;
1892 : }
1893 0 : if (!matchesFinalBoson) nRes++;
1894 0 : }
1895 : }
1896 0 : return nRes;
1897 : }
1898 :
1899 : //--------------------------------------------------------------------------
1900 :
1901 : // Function to report if a resonace decay was found in the
1902 : // 2->2 hard core process
1903 :
1904 : bool HardProcess::hasResInProc(){
1905 :
1906 0 : for(int i =0; i< int(hardIntermediate.size()); ++i)
1907 0 : if (hardIntermediate[i] == 0) return false;
1908 : // Do not count final state bosons as resonaces
1909 0 : for(int i =0; i< int(hardIntermediate.size()); ++i){
1910 0 : for(int j =0; j< int(hardOutgoing1.size()); ++j){
1911 0 : if ( hardIntermediate[i] == hardOutgoing1[j] )
1912 0 : return false;
1913 : }
1914 0 : for(int j =0; j< int(hardOutgoing2.size()); ++j){
1915 0 : if ( hardIntermediate[i] == hardOutgoing2[j] )
1916 0 : return false;
1917 : }
1918 : }
1919 0 : return true;
1920 0 : }
1921 :
1922 : //--------------------------------------------------------------------------
1923 :
1924 : // Function to print the hard process (for debug)
1925 :
1926 : void HardProcess::list() const {
1927 0 : cout << " Hard Process: ";
1928 0 : cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
1929 0 : cout << " \t -----> \t ";
1930 0 : for(int i =0; i < int(hardIntermediate.size());++i)
1931 0 : cout << hardIntermediate[i] << " ";
1932 0 : cout << " \t -----> \t ";
1933 0 : for(int i =0; i < int(hardOutgoing1.size());++i)
1934 0 : cout << hardOutgoing1[i] << " ";
1935 0 : for(int i =0; i < int(hardOutgoing2.size());++i)
1936 0 : cout << hardOutgoing2[i] << " ";
1937 0 : cout << endl;
1938 0 : }
1939 :
1940 : //--------------------------------------------------------------------------
1941 :
1942 : // Function to list the hard process candiates in the matrix element state
1943 : // (for debug)
1944 :
1945 : void HardProcess::listCandidates() const {
1946 0 : cout << " Hard Process candidates: ";
1947 0 : cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
1948 0 : cout << " \t -----> \t ";
1949 0 : for(int i =0; i < int(PosIntermediate.size());++i)
1950 0 : cout << PosIntermediate[i] << " ";
1951 0 : cout << " \t -----> \t ";
1952 0 : for(int i =0; i < int(PosOutgoing1.size());++i)
1953 0 : cout << PosOutgoing1[i] << " ";
1954 0 : for(int i =0; i < int(PosOutgoing2.size());++i)
1955 0 : cout << PosOutgoing2[i] << " ";
1956 0 : cout << endl;
1957 0 : }
1958 :
1959 : //--------------------------------------------------------------------------
1960 :
1961 : // Function to clear hard process information
1962 :
1963 : void HardProcess::clear() {
1964 :
1965 : // Clear flavour of the first incoming particle
1966 0 : hardIncoming1 = hardIncoming2 = 0;
1967 : // Clear outgoing particles
1968 0 : hardOutgoing1.resize(0);
1969 0 : hardOutgoing2.resize(0);
1970 : // Clear intermediate bosons in the hard 2->2
1971 0 : hardIntermediate.resize(0);
1972 :
1973 : // Clear reference event
1974 0 : state.clear();
1975 :
1976 : // Clear potential positions of outgoing particles in reference event
1977 0 : PosOutgoing1.resize(0);
1978 0 : PosOutgoing2.resize(0);
1979 : // Clear potential positions of intermediate bosons in reference event
1980 0 : PosIntermediate.resize(0);
1981 :
1982 : // Clear merging scale read from LHE file
1983 0 : tms = 0.;
1984 :
1985 0 : }
1986 :
1987 : //==========================================================================
1988 :
1989 : // The MergingHooks class.
1990 :
1991 : //--------------------------------------------------------------------------
1992 :
1993 : // Initialise MergingHooks class
1994 :
1995 : void MergingHooks::init( Settings settings, Info* infoPtrIn,
1996 : ParticleData* particleDataPtrIn, PartonSystems* partonSystemsPtrIn,
1997 : ostream& os){
1998 :
1999 : // Save pointers
2000 0 : infoPtr = infoPtrIn;
2001 0 : particleDataPtr = particleDataPtrIn;
2002 0 : partonSystemsPtr = partonSystemsPtrIn;
2003 :
2004 : // Initialise AlphaS objects for reweighting
2005 0 : double alphaSvalueFSR = settings.parm("TimeShower:alphaSvalue");
2006 0 : int alphaSorderFSR = settings.mode("TimeShower:alphaSorder");
2007 0 : int alphaSnfmax = settings.mode("StandardModel:alphaSnfmax");
2008 0 : int alphaSuseCMWFSR= settings.flag("TimeShower:alphaSuseCMW");
2009 0 : AlphaS_FSRSave.init(alphaSvalueFSR, alphaSorderFSR, alphaSnfmax,
2010 : alphaSuseCMWFSR);
2011 0 : double alphaSvalueISR = settings.parm("SpaceShower:alphaSvalue");
2012 0 : int alphaSorderISR = settings.mode("SpaceShower:alphaSorder");
2013 0 : int alphaSuseCMWISR= settings.flag("SpaceShower:alphaSuseCMW");
2014 0 : AlphaS_ISRSave.init(alphaSvalueISR, alphaSorderISR, alphaSnfmax,
2015 : alphaSuseCMWISR);
2016 :
2017 : // Initialise AlphaEM objects for reweighting
2018 0 : int alphaEMFSRorder = settings.mode("TimeShower:alphaEMorder");
2019 0 : AlphaEM_FSRSave.init(alphaEMFSRorder, &settings);
2020 :
2021 : // Initialise merging switches
2022 0 : doUserMergingSave = settings.flag("Merging:doUserMerging");
2023 : // Initialise automated MadGraph kT merging
2024 0 : doMGMergingSave = settings.flag("Merging:doMGMerging");
2025 : // Initialise kT merging
2026 0 : doKTMergingSave = settings.flag("Merging:doKTMerging");
2027 : // Initialise evolution-pT merging
2028 0 : doPTLundMergingSave = settings.flag("Merging:doPTLundMerging");
2029 : // Initialise \Delta_R_{ij}, pT_i Q_{ij} merging
2030 0 : doCutBasedMergingSave = settings.flag("Merging:doCutBasedMerging");
2031 : // Initialise exact definition of kT
2032 0 : ktTypeSave = settings.mode("Merging:ktType");
2033 :
2034 : // Initialise NL3 switches.
2035 0 : doNL3TreeSave = settings.flag("Merging:doNL3Tree");
2036 0 : doNL3LoopSave = settings.flag("Merging:doNL3Loop");
2037 0 : doNL3SubtSave = settings.flag("Merging:doNL3Subt");
2038 0 : bool doNL3 = doNL3TreeSave || doNL3LoopSave || doNL3SubtSave;
2039 :
2040 : // Initialise UNLOPS switches.
2041 0 : doUNLOPSTreeSave = settings.flag("Merging:doUNLOPSTree");
2042 0 : doUNLOPSLoopSave = settings.flag("Merging:doUNLOPSLoop");
2043 0 : doUNLOPSSubtSave = settings.flag("Merging:doUNLOPSSubt");
2044 0 : doUNLOPSSubtNLOSave = settings.flag("Merging:doUNLOPSSubtNLO");
2045 0 : bool doUNLOPS = doUNLOPSTreeSave || doUNLOPSLoopSave
2046 0 : || doUNLOPSSubtSave || doUNLOPSSubtNLOSave;
2047 :
2048 : // Initialise UMEPS switches
2049 0 : doUMEPSTreeSave = settings.flag("Merging:doUMEPSTree");
2050 0 : doUMEPSSubtSave = settings.flag("Merging:doUMEPSSubt");
2051 0 : nReclusterSave = settings.mode("Merging:nRecluster");
2052 0 : nQuarksMergeSave = settings.mode("Merging:nQuarksMerge");
2053 0 : bool doUMEPS = doUMEPSTreeSave || doUMEPSSubtSave;
2054 :
2055 : // Flag to only do phase space cut.
2056 0 : doEstimateXSection = settings.flag("Merging:doXSectionEstimate");
2057 :
2058 : // Get core process from user input
2059 0 : processSave = settings.word("Merging:Process");
2060 :
2061 : // Clear hard process
2062 0 : hardProcess.clear();
2063 :
2064 : // Initialise input event.
2065 0 : inputEvent.init("(hard process)", particleDataPtr);
2066 0 : doRemoveDecayProducts = settings.flag("Merging:mayRemoveDecayProducts");
2067 :
2068 : // Initialise the hard process
2069 0 : if ( doMGMergingSave )
2070 0 : hardProcess.initOnLHEF(lheInputFile, particleDataPtr);
2071 : else
2072 0 : hardProcess.initOnProcess(processSave, particleDataPtr);
2073 :
2074 : // Parameters for reconstruction of evolution scales
2075 0 : includeMassiveSave = settings.flag("Merging:includeMassive");
2076 0 : enforceStrongOrderingSave = settings.flag("Merging:enforceStrongOrdering");
2077 0 : scaleSeparationFactorSave = settings.parm("Merging:scaleSeparationFactor");
2078 0 : orderInRapiditySave = settings.flag("Merging:orderInRapidity");
2079 :
2080 : // Parameters for choosing history probabilistically
2081 0 : nonJoinedNormSave = settings.parm("Merging:nonJoinedNorm");
2082 0 : fsrInRecNormSave = settings.parm("Merging:fsrInRecNorm");
2083 0 : pickByFullPSave = settings.flag("Merging:pickByFullP");
2084 0 : pickByPoPT2Save = settings.flag("Merging:pickByPoPT2");
2085 0 : includeRedundantSave = settings.flag("Merging:includeRedundant");
2086 :
2087 : // Parameters for scale choices
2088 0 : unorderedScalePrescipSave =
2089 0 : settings.mode("Merging:unorderedScalePrescrip");
2090 0 : unorderedASscalePrescipSave =
2091 0 : settings.mode("Merging:unorderedASscalePrescrip");
2092 0 : unorderedPDFscalePrescipSave =
2093 0 : settings.mode("Merging:unorderedPDFscalePrescrip");
2094 0 : incompleteScalePrescipSave =
2095 0 : settings.mode("Merging:incompleteScalePrescrip");
2096 :
2097 : // Parameter for allowing swapping of one colour index while reclustering
2098 0 : allowColourShufflingSave =
2099 0 : settings.flag("Merging:allowColourShuffling");
2100 :
2101 : // Parameters to allow setting hard process scales to default (dynamical)
2102 : // Pythia values.
2103 0 : resetHardQRenSave = settings.flag("Merging:usePythiaQRenHard");
2104 0 : resetHardQFacSave = settings.flag("Merging:usePythiaQFacHard");
2105 :
2106 : // Parameters for choosing history by sum(|pT|)
2107 0 : pickBySumPTSave = settings.flag("Merging:pickBySumPT");
2108 0 : herwigAcollFSRSave = settings.parm("Merging:aCollFSR");
2109 0 : herwigAcollISRSave = settings.parm("Merging:aCollISR");
2110 :
2111 : // Information on the shower cut-off scale
2112 0 : pT0ISRSave = settings.parm("SpaceShower:pT0Ref");
2113 0 : pTcutSave = settings.parm("SpaceShower:pTmin");
2114 0 : pTcutSave = max(pTcutSave,pT0ISRSave);
2115 :
2116 : // Initialise CKKWL weight
2117 0 : weightCKKWLSave = 1.;
2118 0 : weightFIRSTSave = 0.;
2119 0 : nMinMPISave = 100;
2120 0 : muMISave = -1.;
2121 :
2122 : // Initialise merging scale
2123 0 : tmsValueSave = 0.;
2124 0 : tmsListSave.resize(0);
2125 :
2126 0 : kFactor0jSave = settings.parm("Merging:kFactor0j");
2127 0 : kFactor1jSave = settings.parm("Merging:kFactor1j");
2128 0 : kFactor2jSave = settings.parm("Merging:kFactor2j");
2129 :
2130 0 : muFSave = settings.parm("Merging:muFac");
2131 0 : muRSave = settings.parm("Merging:muRen");
2132 0 : muFinMESave = settings.parm("Merging:muFacInME");
2133 0 : muRinMESave = settings.parm("Merging:muRenInME");
2134 :
2135 0 : doWClusteringSave = settings.flag("Merging:allowWClustering");
2136 0 : doSQCDClusteringSave = settings.flag("Merging:allowSQCDClustering");
2137 0 : DparameterSave = settings.parm("Merging:Dparameter");
2138 :
2139 : // Save merging scale on maximal number of jets
2140 0 : if ( doKTMergingSave || doUserMergingSave || doPTLundMergingSave
2141 0 : || doUMEPS ) {
2142 : // Read merging scale (defined in kT) from input parameter.
2143 0 : tmsValueSave = settings.parm("Merging:TMS");
2144 0 : nJetMaxSave = settings.mode("Merging:nJetMax");
2145 0 : nJetMaxNLOSave = -1;
2146 0 : } else if (doMGMergingSave) {
2147 : // Read merging scale (defined in kT) from LHE file.
2148 0 : tmsValueSave = hardProcess.tms;
2149 0 : nJetMaxSave = settings.mode("Merging:nJetMax");
2150 0 : nJetMaxNLOSave = -1;
2151 0 : } else if (doCutBasedMergingSave) {
2152 :
2153 : // Save list of cuts defining the merging scale.
2154 0 : nJetMaxSave = settings.mode("Merging:nJetMax");
2155 0 : nJetMaxNLOSave = -1;
2156 : // Write tms cut values to list of cut values,
2157 : // ordered by DeltaR_{ij}, pT_{i}, Q_{ij}.
2158 0 : tmsListSave.resize(0);
2159 0 : double drms = settings.parm("Merging:dRijMS");
2160 0 : double ptms = settings.parm("Merging:pTiMS");
2161 0 : double qms = settings.parm("Merging:QijMS");
2162 0 : tmsListSave.push_back(drms);
2163 0 : tmsListSave.push_back(ptms);
2164 0 : tmsListSave.push_back(qms);
2165 :
2166 0 : }
2167 :
2168 0 : if ( doNL3 || doUNLOPS || doEstimateXSection ) {
2169 0 : tmsValueSave = settings.parm("Merging:TMS");
2170 0 : nJetMaxSave = settings.mode("Merging:nJetMax");
2171 0 : nJetMaxNLOSave = settings.mode("Merging:nJetMaxNLO");
2172 0 : }
2173 :
2174 0 : bool writeBanner = doKTMergingSave || doMGMergingSave
2175 0 : || doUserMergingSave
2176 0 : || doNL3 || doUNLOPS || doUMEPS
2177 0 : || doPTLundMergingSave || doCutBasedMergingSave;
2178 :
2179 0 : if (!writeBanner) return;
2180 :
2181 : // Write banner
2182 0 : os << "\n *------------------ MEPS Merging Initialization ---------------"
2183 0 : << "---*";
2184 0 : os << "\n | "
2185 0 : << " |\n";
2186 0 : if ( doKTMergingSave || doMGMergingSave || doUserMergingSave
2187 0 : || doPTLundMergingSave || doCutBasedMergingSave )
2188 0 : os << " | CKKW-L merge "
2189 0 : << " |\n"
2190 0 : << " |"<< setw(34) << processSave << " with up to"
2191 0 : << setw(3) << nJetMaxSave << " additional jets |\n";
2192 0 : else if ( doNL3 )
2193 0 : os << " | NL3 merge "
2194 0 : << " |\n"
2195 0 : << " |" << setw(31) << processSave << " with jets up to"
2196 0 : << setw(3) << nJetMaxNLOSave << " correct to NLO |\n"
2197 0 : << " | and up to" << setw(3) << nJetMaxSave
2198 0 : << " additional jets included by CKKW-L merging at LO |\n";
2199 0 : else if ( doUNLOPS )
2200 0 : os << " | UNLOPS merge "
2201 0 : << " |\n"
2202 0 : << " |" << setw(31) << processSave << " with jets up to"
2203 0 : << setw(3)<< nJetMaxNLOSave << " correct to NLO |\n"
2204 0 : << " | and up to" << setw(3) << nJetMaxSave
2205 0 : << " additional jets included by UMEPS merging at LO |\n";
2206 0 : else if ( doUMEPS )
2207 0 : os << " | UMEPS merge "
2208 0 : << " |\n"
2209 0 : << " |" << setw(34) << processSave << " with up to"
2210 0 : << setw(3) << nJetMaxSave << " additional jets |\n";
2211 :
2212 0 : if ( doKTMergingSave )
2213 0 : os << " | Merging scale is defined in kT, with value ktMS = "
2214 0 : << tmsValueSave << " GeV";
2215 0 : else if ( doMGMergingSave )
2216 0 : os << " | Perform automanted MG/ME merging \n"
2217 0 : << " | Merging scale is defined in kT, with value ktMS = "
2218 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2219 0 : else if ( doUserMergingSave )
2220 0 : os << " | Merging scale is defined by the user, with value tMS = "
2221 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " |";
2222 0 : else if ( doPTLundMergingSave )
2223 0 : os << " | Merging scale is defined by Lund pT, with value tMS = "
2224 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2225 0 : else if ( doCutBasedMergingSave )
2226 0 : os << " | Merging scale is defined by combination of Delta R_{ij}, pT_i "
2227 0 : << " |\n"
2228 0 : << " | and Q_{ij} cut, with values "
2229 0 : << " |\n"
2230 0 : << " | Delta R_{ij,min} = "
2231 0 : << setw(7) << scientific << setprecision(2) << tmsListSave[0]
2232 0 : << " |\n"
2233 0 : << " | pT_{i,min} = "
2234 0 : << setw(6) << fixed << setprecision(1) << tmsListSave[1]
2235 0 : << " GeV |\n"
2236 0 : << " | Q_{ij,min} = "
2237 0 : << setw(6) << fixed << setprecision(1) << tmsListSave[2]
2238 0 : << " GeV |";
2239 0 : else if ( doNL3TreeSave )
2240 0 : os << " | Generate tree-level O(alpha_s)-subtracted events "
2241 0 : << " |\n"
2242 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2243 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2244 0 : else if ( doNL3LoopSave )
2245 0 : os << " | Generate virtual correction unit-weight events "
2246 0 : << " |\n"
2247 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2248 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2249 0 : else if ( doNL3SubtSave )
2250 0 : os << " | Generate reclustered tree-level events "
2251 0 : << " |\n"
2252 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2253 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2254 0 : else if ( doUNLOPSTreeSave )
2255 0 : os << " | Generate tree-level O(alpha_s)-subtracted events "
2256 0 : << " |\n"
2257 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2258 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2259 0 : else if ( doUNLOPSLoopSave )
2260 0 : os << " | Generate virtual correction unit-weight events "
2261 0 : << " |\n"
2262 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2263 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2264 0 : else if ( doUNLOPSSubtSave )
2265 0 : os << " | Generate reclustered tree-level events "
2266 0 : << " |\n"
2267 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2268 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2269 0 : else if ( doUNLOPSSubtNLOSave )
2270 0 : os << " | Generate reclustered loop-level events "
2271 0 : << " |\n"
2272 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2273 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2274 0 : else if ( doUMEPSTreeSave )
2275 0 : os << " | Generate tree-level events "
2276 0 : << " |\n"
2277 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2278 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2279 0 : else if ( doUMEPSSubtSave )
2280 0 : os << " | Generate reclustered tree-level events "
2281 0 : << " |\n"
2282 0 : << " | Merging scale is defined by Lund pT, with value tMS = "
2283 0 : << setw(6) << fixed << setprecision(1) << tmsValueSave << " GeV |";
2284 :
2285 0 : os << "\n | "
2286 0 : << " |";
2287 0 : os << "\n *-------------- END MEPS Merging Initialization ---------------"
2288 0 : << "---*\n\n";
2289 :
2290 0 : }
2291 :
2292 : //--------------------------------------------------------------------------
2293 :
2294 : // Function to check if emission should be rejected.
2295 :
2296 : bool MergingHooks::doVetoEmission( const Event& event) {
2297 :
2298 : // Do nothing in trial showers, or after first step.
2299 0 : if ( doIgnoreEmissionsSave ) return false;
2300 :
2301 : // Do nothing in CKKW-L
2302 0 : if ( doUserMerging() || doMGMerging() || doKTMerging()
2303 0 : || doPTLundMerging() || doCutBasedMerging() )
2304 0 : return false;
2305 :
2306 : // For NLO merging, count and veto emissions above the merging scale
2307 : bool veto = false;
2308 : // Get number of clustering steps
2309 0 : int nSteps = getNumberOfClusteringSteps(event);
2310 : // Get merging scale in current event
2311 0 : double tnow = rhoms( event, false);
2312 :
2313 : // Get maximal number of additional jets
2314 0 : int nJetMax = nMaxJets();
2315 : // Always remove emissions above the merging scale for
2316 : // samples containing reclusterings!
2317 0 : if ( nRecluster() > 0 ) nSteps = max(1, min(nJetMax-2, 1));
2318 : // Check veto condition
2319 0 : if ( nSteps - 1 < nJetMax && nSteps >= 1 && tnow > tms() ) veto = true;
2320 :
2321 : // Do not veto if state already includes MPI.
2322 0 : if ( infoPtr->nMPI() > 1 ) veto = false;
2323 :
2324 : // When performing NL3 merging of tree-level events, reset the
2325 : // CKKWL weight.
2326 0 : if ( veto && doNL3Tree() ) setWeightCKKWL(0.);
2327 :
2328 : // If the emission is allowed, do not check any further emissions
2329 0 : if ( !veto ) doIgnoreEmissionsSave = true;
2330 :
2331 : // Done
2332 : return veto;
2333 :
2334 0 : }
2335 :
2336 : //--------------------------------------------------------------------------
2337 :
2338 : // Function to check if emission should be rejected.
2339 :
2340 : bool MergingHooks::doVetoStep( const Event& process, const Event& event,
2341 : bool doResonance ) {
2342 :
2343 : // Do nothing in trial showers, or after first step.
2344 0 : if ( doIgnoreStepSave && !doResonance ) return false;
2345 :
2346 : // Do nothing in UMEPS or UNLOPS
2347 0 : if ( doUMEPSTree() || doUMEPSSubt() || doUMEPSMerging() || doUNLOPSTree()
2348 0 : || doUNLOPSLoop() || doUNLOPSSubt() || doUNLOPSSubtNLO()
2349 0 : || doUNLOPSMerging() )
2350 0 : return false;
2351 :
2352 : // Get number of clustering steps. If necessary, remove resonance
2353 : // decay products first.
2354 0 : int nSteps = (doResonance) ? getNumberOfClusteringSteps(process)
2355 0 : : getNumberOfClusteringSteps( bareEvent( process, false) );
2356 : // Get maximal number of additional jets.
2357 0 : int nJetMax = nMaxJets();
2358 : // Get merging scale in current event.
2359 0 : double tnow = tmsNow( event );
2360 :
2361 : // For non-resonant showers, simply check veto. If event should indeed be
2362 : // vetoed, save the current pT and weights in case the veto needs to be
2363 : // revoked.
2364 0 : if ( !doResonance ) {
2365 :
2366 : // Store pT to check if veto needs to be revoked later.
2367 0 : pTsave = infoPtr->pTnow();
2368 0 : if ( nRecluster() == 1) nSteps--;
2369 :
2370 : // Check merging veto condition.
2371 : bool veto = false;
2372 0 : if ( nSteps > nMaxJetsNLO() && nSteps < nJetMax && tnow > tms() ) {
2373 : // Set weight to zero if event should be vetoed.
2374 0 : weightCKKWL1Save = 0.;
2375 : // Save weight before veto, in case veto needs to be revoked.
2376 0 : weightCKKWL2Save = getWeightCKKWL();
2377 : // Reset stored weights.
2378 0 : setWeightCKKWL(0.);
2379 : veto = true;
2380 0 : }
2381 :
2382 : // Done
2383 0 : return veto;
2384 :
2385 : // For resonant showers, check if any previous veto should be revoked.
2386 : // This means we treat showers off resonance decay products identical to
2387 : // MPI: If a hard resonance emission has been produced, the event should
2388 : // have been kept. Following this reasoning, it might be necessary to revoke
2389 : // any previous veto.
2390 : } else {
2391 :
2392 : // Initialise switch to revoke vetoing.
2393 : bool revokeVeto = false;
2394 :
2395 : // Nothing to check if pTsave was not stored, i.e. no emission to
2396 : // possibly veto was recorded.
2397 : // Only allow revoking the veto for diboson processes, with resonant
2398 : // electroweak bosons
2399 0 : bool check = (nHardInLeptons() == 0)&& (nHardOutLeptons() == 2)
2400 0 : && (nHardOutPartons() == 2);
2401 :
2402 : // For current purpose only!!!
2403 : check = false;
2404 :
2405 : // For hadronic resonance decays at hadron colliders, do not veto
2406 : // events with a hard emission of the resonance decay products,
2407 : // since such states are not included in the matrix element
2408 0 : if ( pTsave > 0. && check ) {
2409 :
2410 : // Check how many resonance decay systems are allowed
2411 0 : int nResNow = nResInCurrent();
2412 :
2413 : // Find systems really containing only emission off a resonance
2414 : // decay
2415 0 : vector<int>goodSys;
2416 : // Resonance decay systems are considered last, thus at the end of
2417 : // the list
2418 0 : int sysSize = partonSystemsPtr->sizeSys();
2419 0 : for ( int i=0; i < nResNow; ++i ) {
2420 0 : if ( partonSystemsPtr->sizeOut(sysSize - 1 - i) == 3 )
2421 0 : goodSys.push_back(sysSize - 1 - i);
2422 : }
2423 :
2424 : // Check the members of the resonance decay systems
2425 0 : for ( int i=0; i < int(goodSys.size()); ++i ) {
2426 :
2427 : // Save the (three) members of the resonance decay system
2428 0 : int iMem1 = partonSystemsPtr->getOut(goodSys[i],0);
2429 0 : int iMem2 = partonSystemsPtr->getOut(goodSys[i],1);
2430 0 : int iMem3 = partonSystemsPtr->getOut(goodSys[i],2);
2431 :
2432 : // Now find emitted gluon or gamma
2433 0 : int iEmtGlue = ((event[iMem1].id() == 21) ? iMem1
2434 0 : : ((event[iMem2].id() == 21) ? iMem2
2435 0 : : ((event[iMem3].id() == 21) ? iMem3: 0)));
2436 0 : int iEmtGamm = ((event[iMem1].id() == 22) ? iMem1
2437 0 : : ((event[iMem2].id() == 22) ? iMem2
2438 0 : : ((event[iMem3].id() == 22) ? iMem3: 0)));
2439 : // Per system, only one emission
2440 0 : int iEmt = (iEmtGlue != 0) ? iEmtGlue : iEmtGamm;
2441 :
2442 : int iRad = 0;
2443 : int iRec = 0;
2444 0 : if ( iEmt == iMem1 ) {
2445 0 : iRad = (event[iMem2].mother1() != event[iMem2].mother2())
2446 : ? iMem2 : iMem3;
2447 0 : iRec = (event[iMem3].mother1() == event[iMem3].mother2())
2448 : ? iMem3 : iMem2;
2449 0 : } else if ( iEmt == iMem2 ) {
2450 0 : iRad = (event[iMem1].mother1() != event[iMem1].mother2())
2451 : ? iMem1 : iMem3;
2452 0 : iRec = (event[iMem3].mother1() == event[iMem3].mother2())
2453 : ? iMem3 : iMem1;
2454 0 : } else {
2455 0 : iRad = (event[iMem1].mother1() != event[iMem1].mother2())
2456 : ? iMem1 : iMem2;
2457 0 : iRec = (event[iMem2].mother1() == event[iMem2].mother2())
2458 : ? iMem2 : iMem1;
2459 : }
2460 :
2461 0 : double pTres = rhoPythia(event[iRad], event[iEmt], event[iRec], 1);
2462 :
2463 : // Revoke previous veto of last emission if a splitting of the
2464 : // resonance produced a harder parton, i.e. we are inside the
2465 : // PS region
2466 0 : if ( pTres > pTsave ) {
2467 : revokeVeto = true;
2468 : // Do nothing (i.e. allow other first emission veto) for soft
2469 : // splitting
2470 0 : } else {
2471 : revokeVeto = false;
2472 : }
2473 : // Done with one system
2474 : }
2475 : // Done with all systems
2476 0 : }
2477 :
2478 : // Check veto condition
2479 : bool veto = false;
2480 0 : if ( revokeVeto && check ) {
2481 0 : setWeightCKKWL(weightCKKWL2Save);
2482 0 : } else if ( check ) {
2483 0 : setWeightCKKWL(weightCKKWL1Save);
2484 0 : if ( weightCKKWL1Save == 0. ) veto = true;
2485 : }
2486 :
2487 : // Check veto condition.
2488 0 : if ( !check && nSteps > nMaxJetsNLO() && nSteps < nJetMax && tnow > tms()){
2489 : // Set stored weights to zero.
2490 0 : setWeightCKKWL(0.);
2491 : // Now allow veto.
2492 : veto = true;
2493 0 : }
2494 :
2495 : // If the emission is allowed, do not check any further emissions
2496 0 : if ( !veto || !doIgnoreStepSave ) doIgnoreStepSave = true;
2497 :
2498 : // Done
2499 : return veto;
2500 :
2501 : }
2502 :
2503 : // Done
2504 : return false;
2505 :
2506 0 : }
2507 :
2508 : //--------------------------------------------------------------------------
2509 :
2510 : // Return event stripped from decay products.
2511 :
2512 : Event MergingHooks::bareEvent(const Event& inputEventIn,
2513 : bool storeInputEvent ) {
2514 :
2515 : // Find and detach decay products.
2516 0 : Event newProcess = Event();
2517 0 : newProcess.init("(hard process-modified)", particleDataPtr);
2518 :
2519 : // If desired, store input event.
2520 0 : if ( storeInputEvent ) {
2521 0 : resonances.resize(0);
2522 0 : inputEvent.clear();
2523 0 : for (int i = 0; i < inputEventIn.size(); ++i)
2524 0 : inputEvent.append( inputEventIn[i] );
2525 0 : for (int i = 0; i < inputEventIn.sizeJunction(); ++i)
2526 0 : inputEvent.appendJunction( inputEventIn.getJunction(i) );
2527 0 : inputEvent.saveSize();
2528 0 : inputEvent.saveJunctionSize();
2529 0 : }
2530 :
2531 : // Now remove decay products.
2532 0 : if ( doRemoveDecayProducts ) {
2533 :
2534 : // Add the beam and initial partons to the event record.
2535 0 : for (int i = 0; i < inputEventIn.size(); ++ i) {
2536 0 : if ( inputEventIn[i].mother1() > 4
2537 0 : || inputEventIn[i].statusAbs() == 22
2538 0 : || inputEventIn[i].statusAbs() == 23)
2539 0 : break;
2540 0 : newProcess.append(inputEventIn[i]);
2541 : }
2542 :
2543 : // Add the intermediate particles to the event record.
2544 0 : for (int i = 0; i < inputEventIn.size(); ++ i) {
2545 0 : if (inputEventIn[i].mother1() > 4) break;
2546 0 : if ( inputEventIn[i].status() == -22) {
2547 0 : int j = newProcess.append(inputEventIn[i]);
2548 0 : newProcess[j].statusPos();
2549 0 : if ( storeInputEvent ) resonances.push_back( make_pair(j, i) );
2550 0 : newProcess[j].daughters(0, 0);
2551 0 : }
2552 : }
2553 :
2554 : // Add remaining outgoing particles to the event record.
2555 0 : for (int i = 0; i < inputEventIn.size(); ++ i) {
2556 0 : if (inputEventIn[i].mother1() > 4) break;
2557 0 : if ( inputEventIn[i].statusAbs() != 11
2558 0 : && inputEventIn[i].statusAbs() != 12
2559 0 : && inputEventIn[i].statusAbs() != 21
2560 0 : && inputEventIn[i].statusAbs() != 22)
2561 0 : newProcess.append(inputEventIn[i]);
2562 : }
2563 :
2564 : // Update event colour tag to maximum in whole process.
2565 : int maxColTag = 0;
2566 0 : for (int i = 0; i < inputEventIn.size(); ++ i) {
2567 0 : if ( inputEventIn[i].col() > maxColTag )
2568 0 : maxColTag = inputEventIn[i].col();
2569 0 : if ( inputEventIn[i].acol() > maxColTag )
2570 0 : maxColTag = inputEventIn[i].acol();
2571 : }
2572 0 : newProcess.initColTag(maxColTag);
2573 :
2574 : // Copy junctions from process to newProcess.
2575 0 : for (int i = 0; i < inputEventIn.sizeJunction(); ++i)
2576 0 : newProcess.appendJunction( inputEventIn.getJunction(i));
2577 :
2578 0 : newProcess.saveSize();
2579 0 : newProcess.saveJunctionSize();
2580 :
2581 0 : } else {
2582 0 : newProcess = inputEventIn;
2583 : }
2584 :
2585 : // Remember scale
2586 0 : newProcess.scale( inputEventIn.scale() );
2587 :
2588 : // Done
2589 : return newProcess;
2590 :
2591 0 : }
2592 :
2593 : //--------------------------------------------------------------------------
2594 :
2595 : // Write event with decay products attached to argument. Only possible if an
2596 : // input event with decay producs had been stored before.
2597 :
2598 : bool MergingHooks::reattachResonanceDecays(Event& process ) {
2599 :
2600 : // Now reattach the decay products.
2601 0 : if ( doRemoveDecayProducts && inputEvent.size() > 0 ) {
2602 :
2603 0 : int sizeBef = process.size();
2604 : // Vector of resonances for which the decay products were already attached.
2605 0 : vector<int> iAftChecked;
2606 : // Reset daughters and status of intermediate particles.
2607 0 : for ( int i = 0; i < int(resonances.size()); ++i ) {
2608 0 : for (int j = 0; j < sizeBef; ++j ) {
2609 0 : if ( j != resonances[i].first ) continue;
2610 :
2611 0 : int iOldDaughter1 = inputEvent[resonances[i].second].daughter1();
2612 0 : int iOldDaughter2 = inputEvent[resonances[i].second].daughter2();
2613 :
2614 : // Get momenta in case of reclustering.
2615 0 : int iHardMother = resonances[i].second;
2616 0 : Particle& hardMother = inputEvent[iHardMother];
2617 : // Find current mother copy (after clustering).
2618 : int iAftMother = 0;
2619 0 : for ( int k = 0; k < process.size(); ++k )
2620 0 : if ( process[k].id() == inputEvent[resonances[i].second].id() ) {
2621 : // Only attempt if decays of this resonance were not attached.
2622 : bool checked = false;
2623 0 : for ( int l = 0; l < int(iAftChecked.size()); ++l )
2624 0 : if ( k == iAftChecked[l] )
2625 0 : checked = true;
2626 0 : if ( !checked ) {
2627 0 : iAftChecked.push_back(k);
2628 0 : iAftMother = k;
2629 0 : break;
2630 : }
2631 0 : }
2632 :
2633 0 : Particle& aftMother = process[iAftMother];
2634 :
2635 : // Resonance can have been moved by clustering,
2636 : // so prepare to update colour and momentum information for system.
2637 0 : int colBef = hardMother.col();
2638 0 : int acolBef = hardMother.acol();
2639 0 : int colAft = aftMother.col();
2640 0 : int acolAft = aftMother.acol();
2641 0 : RotBstMatrix M;
2642 0 : M.bst( hardMother.p(), aftMother.p());
2643 :
2644 : // Attach resonance decay products.
2645 : int iNewDaughter1 = 0;
2646 : int iNewDaughter2 = 0;
2647 0 : for ( int k = iOldDaughter1; k <= iOldDaughter2; ++k ) {
2648 0 : if ( k == iOldDaughter1 )
2649 0 : iNewDaughter1 = process.append(inputEvent[k] );
2650 : else
2651 0 : iNewDaughter2 = process.append(inputEvent[k] );
2652 0 : process.back().statusPos();
2653 0 : Particle& now = process.back();
2654 : // Update colour and momentum information.
2655 0 : if (now.col() != 0 && now.col() == colBef ) now.col(colAft);
2656 0 : if (now.acol() != 0 && now.acol() == acolBef) now.acol(acolAft);
2657 0 : now.rotbst( M);
2658 : // Update vertex information.
2659 0 : if (now.hasVertex()) now.vProd( aftMother.vDec() );
2660 : // Update mothers.
2661 0 : now.mothers(iAftMother,0);
2662 : }
2663 :
2664 0 : process[iAftMother].daughters( iNewDaughter1, iNewDaughter2 );
2665 0 : process[iAftMother].statusNeg();
2666 :
2667 : // Loop through event and attach remaining decays
2668 : int iDec = 0;
2669 0 : do {
2670 0 : if ( process[iDec].isFinal() && process[iDec].canDecay()
2671 0 : && process[iDec].mayDecay() && process[iDec].isResonance() ) {
2672 :
2673 0 : int iD1 = process[iDec].daughter1();
2674 0 : int iD2 = process[iDec].daughter2();
2675 :
2676 : // Done if no daughters exist.
2677 0 : if ( iD1 == 0 || iD2 == 0 ) continue;
2678 :
2679 : // Attach daughters.
2680 : int iNewDaughter12 = 0;
2681 : int iNewDaughter22 = 0;
2682 0 : for ( int k = iD1; k <= iD2; ++k ) {
2683 0 : if ( k == iD1 )
2684 0 : iNewDaughter12 = process.append(inputEvent[k] );
2685 : else
2686 0 : iNewDaughter22 = process.append(inputEvent[k] );
2687 0 : process.back().statusPos();
2688 0 : Particle& now = process.back();
2689 : // Update colour and momentum information.
2690 0 : if (now.col() != 0 && now.col() == colBef ) now.col(colAft);
2691 0 : if (now.acol()!= 0 && now.acol()== acolBef) now.acol(acolAft);
2692 0 : now.rotbst( M);
2693 : // Update vertex information.
2694 0 : if (now.hasVertex()) now.vProd( process[iDec].vDec() );
2695 : // Update mothers.
2696 0 : now.mothers(iDec,0);
2697 : }
2698 :
2699 : // Modify mother status and daughters.
2700 0 : process[iDec].status(-22);
2701 0 : process[iDec].daughters(iNewDaughter12, iNewDaughter22);
2702 :
2703 : // End of loop over all entries.
2704 0 : }
2705 0 : } while (++iDec < process.size());
2706 0 : } // End loop over process entries.
2707 : } // End loop over resonances.
2708 :
2709 : // Update event colour tag to maximum in whole process.
2710 : int maxColTag = 0;
2711 0 : for (int i = 0; i < process.size(); ++ i) {
2712 0 : if (process[i].col() > maxColTag) maxColTag = process[i].col();
2713 0 : if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
2714 : }
2715 0 : process.initColTag(maxColTag);
2716 :
2717 0 : }
2718 :
2719 : // Done.
2720 0 : return (doRemoveDecayProducts) ? inputEvent.size() > 0 : true;
2721 :
2722 0 : }
2723 :
2724 : //--------------------------------------------------------------------------
2725 :
2726 : bool MergingHooks::isInHard( int iPos, const Event& event){
2727 :
2728 : // MPI not part of hard process
2729 0 : if ( event[iPos].statusAbs() > 30 && event[iPos].statusAbs() < 40 )
2730 0 : return false;
2731 : // Beam remnants and hadronisation not part of hard process
2732 0 : if ( event[iPos].statusAbs() > 60 )
2733 0 : return false;
2734 :
2735 : // Still MPI: Check that the particle is not due to radiation off MPIs.
2736 : // First get all intermediate MPI partons in the state.
2737 0 : vector<int> mpiParticlePos;
2738 0 : mpiParticlePos.clear();
2739 0 : for ( int i=0; i < event.size(); ++i )
2740 0 : if ( event[i].statusAbs() > 30
2741 0 : && event[i].statusAbs() < 40 )
2742 0 : mpiParticlePos.push_back(i);
2743 : // Disregard any parton iPos that has MPI ancestors.
2744 0 : for ( int i=0; i < int(mpiParticlePos.size()); ++i)
2745 0 : if ( event[iPos].isAncestor( mpiParticlePos[i]) )
2746 0 : return false;
2747 :
2748 : // Disallow other systems.
2749 : // Get sub-system of particles for iPos
2750 0 : int iSys = partonSystemsPtr->getSystemOf(iPos, !event[iPos].isFinal() );
2751 0 : if ( iSys > 0 ) {
2752 : // Check all partons belonging to the same system as iPos. If any is
2753 : // produced in MPI or has MPI ancestors, the whole system is not the
2754 : // hard subprocess, i.e. iPos is not in the hard subprocess.
2755 0 : int sysSize = partonSystemsPtr->sizeAll(iSys);
2756 0 : for ( int i = 0; i < sysSize; ++i ) {
2757 0 : int iPosNow = partonSystemsPtr->getAll( iSys, i );
2758 : // MPI not part of hard process
2759 0 : if ( event[iPosNow].statusAbs() > 30
2760 0 : && event[iPosNow].statusAbs() < 40 )
2761 0 : return false;
2762 : // Disregard any parton iPos that has MPI ancestors.
2763 0 : for ( int j=0; j < int(mpiParticlePos.size()); ++j)
2764 0 : if ( event[iPosNow].isAncestor( mpiParticlePos[j]) )
2765 0 : return false;
2766 : // Beam remnants and hadronisation not part of hard process
2767 0 : if ( event[iPosNow].statusAbs() > 60 )
2768 0 : return false;
2769 0 : }
2770 0 : }
2771 :
2772 : // Check if any ancestor contains the hard incoming partons as daughters.
2773 : // Begin loop to trace upwards from the daughter.
2774 : bool containsInitialParton = false;
2775 : int iUp = iPos;
2776 0 : for ( ; ; ) {
2777 : // If out of range then failed to find match.
2778 0 : if (iUp <= 0 || iUp > event.size()) break;
2779 : // If positive match then done.
2780 0 : if ( iUp == 3 || iUp == 4 ) {
2781 : containsInitialParton = true;
2782 0 : break;
2783 : }
2784 0 : if ( event[iUp].mother1() == 1
2785 0 : && (event[iUp].daughter1() == 3 || event[iUp].daughter2() == 3) ) {
2786 : containsInitialParton = true;
2787 0 : break;
2788 : }
2789 0 : if ( event[iUp].mother1() == 2
2790 0 : && (event[iUp].daughter1() == 4 || event[iUp].daughter2() == 4) ) {
2791 : containsInitialParton = true;
2792 0 : break;
2793 : }
2794 : // If unique mother then keep on moving up the chain.
2795 0 : iUp = event[iUp].mother1();
2796 : }
2797 :
2798 0 : if ( !containsInitialParton ) return false;
2799 :
2800 : // Done
2801 0 : return true;
2802 :
2803 0 : }
2804 :
2805 : //--------------------------------------------------------------------------
2806 :
2807 : // Function to return the number of clustering steps for the current event
2808 :
2809 : int MergingHooks::getNumberOfClusteringSteps(const Event& event){
2810 :
2811 : // Count the number of final state partons
2812 : int nFinalPartons = 0;
2813 0 : for ( int i=0; i < event.size(); ++i)
2814 0 : if ( event[i].isFinal()
2815 0 : && isInHard( i, event)
2816 0 : && (event[i].isQuark() || event[i].isGluon()) )
2817 0 : nFinalPartons++;
2818 :
2819 : // Count the number of final state leptons
2820 : int nFinalLeptons = 0;
2821 0 : for( int i=0; i < event.size(); ++i)
2822 0 : if ( event[i].isFinal() && isInHard( i, event) && event[i].isLepton())
2823 0 : nFinalLeptons++;
2824 :
2825 : // Add neutralinos to number of leptons
2826 0 : for( int i=0; i < event.size(); ++i)
2827 0 : if ( event[i].isFinal() && isInHard( i, event)
2828 0 : && event[i].idAbs() == 1000022)
2829 0 : nFinalLeptons++;
2830 :
2831 : // Add sleptons to number of leptons
2832 0 : for( int i=0; i < event.size(); ++i)
2833 0 : if ( event[i].isFinal() && isInHard( i, event)
2834 0 : && (event[i].idAbs() == 1000011
2835 0 : || event[i].idAbs() == 2000011
2836 0 : || event[i].idAbs() == 1000013
2837 0 : || event[i].idAbs() == 2000013
2838 0 : || event[i].idAbs() == 1000015
2839 0 : || event[i].idAbs() == 2000015) )
2840 0 : nFinalLeptons++;
2841 :
2842 : // Count the number of final state electroweak bosons
2843 : int nFinalBosons = 0;
2844 0 : for( int i=0; i < event.size(); ++i)
2845 0 : if ( event[i].isFinal() && isInHard( i, event)
2846 0 : && ( event[i].idAbs() == 22
2847 0 : || event[i].idAbs() == 23
2848 0 : || event[i].idAbs() == 24
2849 0 : || event[i].idAbs() == 25 ) )
2850 0 : nFinalBosons++;
2851 :
2852 : // Save sum of all final state particles
2853 0 : int nFinal = nFinalPartons + nFinalLeptons
2854 0 : + 2*(nFinalBosons - nHardOutBosons() );
2855 :
2856 : // Return the difference to the core process outgoing particles
2857 0 : return (nFinal - nHardOutPartons() - nHardOutLeptons() );
2858 :
2859 : }
2860 :
2861 : //--------------------------------------------------------------------------
2862 :
2863 : // Function to check if event contains an emission not present in the hard
2864 : // process.
2865 :
2866 : bool MergingHooks::isFirstEmission(const Event& event ) {
2867 :
2868 : // If the beam remnant treatment or hadronisation has already started, do
2869 : // no veto.
2870 0 : for ( int i=0; i < event.size(); ++i)
2871 0 : if ( event[i].statusAbs() > 60 ) return false;
2872 :
2873 : // Count particle types
2874 : int nFinalQuarks = 0;
2875 : int nFinalGluons = 0;
2876 : int nFinalLeptons = 0;
2877 : int nFinalBosons = 0;
2878 : int nFinalPhotons = 0;
2879 : int nFinal = 0;
2880 0 : for( int i=0; i < event.size(); ++i) {
2881 0 : if (event[i].isFinal() && isInHard(i, event) ){
2882 0 : if ( event[i].spinType() == 2 && event[i].colType() == 0)
2883 0 : nFinalLeptons++;
2884 0 : if ( event[i].id() == 23
2885 0 : || event[i].idAbs() == 24
2886 0 : || event[i].id() == 25)
2887 0 : nFinalBosons++;
2888 0 : if ( event[i].id() == 22)
2889 0 : nFinalPhotons++;
2890 0 : if ( event[i].isQuark())
2891 0 : nFinalQuarks++;
2892 0 : if ( event[i].isGluon())
2893 0 : nFinalGluons++;
2894 0 : if ( !event[i].isDiquark() )
2895 0 : nFinal++;
2896 : }
2897 : }
2898 :
2899 : // Return highest value if the event does not contain any final state
2900 : // coloured particles.
2901 0 : if (nFinalQuarks + nFinalGluons == 0) return false;
2902 :
2903 : // Use MergingHooks functions to get information on the hard process.
2904 0 : int nLeptons = nHardOutLeptons();
2905 :
2906 : // The state is already in the PS region if the number of leptons had been
2907 : // increased bt QED splittings.
2908 0 : if (nFinalLeptons > nLeptons) return false;
2909 :
2910 : // If the mumber of photons if larger than in the hard process, QED
2911 : // radiation has pushed the state into the PS region.
2912 : int nPhotons = 0;
2913 0 : for(int i =0; i< int(hardProcess.hardOutgoing1.size()); ++i)
2914 0 : if (hardProcess.hardOutgoing1[i] == 22)
2915 0 : nPhotons++;
2916 0 : for(int i =0; i< int(hardProcess.hardOutgoing2.size()); ++i)
2917 0 : if (hardProcess.hardOutgoing2[i] == 22)
2918 0 : nPhotons++;
2919 0 : if (nFinalPhotons > nPhotons) return false;
2920 :
2921 : // Done
2922 0 : return true;
2923 0 : }
2924 :
2925 : //--------------------------------------------------------------------------
2926 :
2927 : // Function to set the correct starting scales of the shower
2928 :
2929 : // Set starting scales
2930 : bool MergingHooks::setShowerStartingScales( bool isTrial,
2931 : bool doMergeFirstEmm, double& pTscaleIn, const Event& event,
2932 : double& pTmaxFSRIn, bool& limitPTmaxFSRIn,
2933 : double& pTmaxISRIn, bool& limitPTmaxISRIn,
2934 : double& pTmaxMPIIn, bool& limitPTmaxMPIIn ) {
2935 :
2936 : // MPI treated differently in case of qcd djiet merging, since hard MPI
2937 : // can be misidentified as hard process.
2938 0 : bool isPureQCD = ( getProcessString().compare("pp>jj") == 0 );
2939 0 : int nSteps = getNumberOfClusteringSteps( bareEvent(event, false) );
2940 0 : double scale = event.scale();
2941 :
2942 0 : bool doRecluster = doUMEPSSubt() || doNL3Subt() || doUNLOPSSubt()
2943 0 : || doUNLOPSSubtNLO();
2944 :
2945 : // Set correct starting scales for trial showers.
2946 0 : if ( isTrial ) {
2947 : // Reset shower scales.
2948 0 : pTmaxISRIn = pTmaxFSRIn = scale;
2949 0 : if ( limitPTmaxISRIn ) pTmaxISRIn = min(scale,muF());
2950 0 : if ( limitPTmaxFSRIn ) pTmaxFSRIn = min(scale,muF());
2951 : // Reset MPI scale.
2952 0 : if ( !isPureQCD ) pTmaxMPIIn = scale;
2953 0 : else pTmaxMPIIn = infoPtr->eCM();
2954 : // Reset phase space limitation flags
2955 0 : if ( pTscaleIn < infoPtr->eCM() ) {
2956 0 : limitPTmaxISRIn = limitPTmaxFSRIn = true;
2957 0 : if ( !isPureQCD ) limitPTmaxMPIIn = true;
2958 0 : else limitPTmaxMPIIn = false;
2959 : }
2960 : }
2961 :
2962 : // Reset starting scales.
2963 0 : if ( isPureQCD && doMergeFirstEmm ) {
2964 : // Set correct shower starting scales.
2965 0 : if ( nSteps == 0 && !doRecluster && pTscaleIn < infoPtr->eCM() ) {
2966 0 : pTmaxMPIIn = infoPtr->eCM();
2967 0 : limitPTmaxMPIIn = false;
2968 0 : }
2969 : }
2970 :
2971 : // Reset starting scales.
2972 0 : if ( !isPureQCD && doMergeFirstEmm ) {
2973 : // Set correct shower starting scales.
2974 0 : if ( nSteps == 0 && !doRecluster ) {
2975 0 : pTscaleIn = infoPtr->eCM();
2976 0 : limitPTmaxMPIIn = false;
2977 0 : }
2978 0 : if ( limitPTmaxISRIn ) pTscaleIn = min(scale,muF());
2979 0 : if ( limitPTmaxFSRIn ) pTscaleIn = min(scale,muF());
2980 0 : pTmaxISRIn = pTmaxFSRIn = pTscaleIn;
2981 : // Reset MPI starting scale. Standard treatment
2982 0 : if ( pTscaleIn < infoPtr->eCM()
2983 0 : && !(nSteps == 0 && !doRecluster
2984 0 : && (limitPTmaxISRIn || limitPTmaxFSRIn)) ) {
2985 0 : pTmaxMPIIn = pTscaleIn;
2986 0 : limitPTmaxMPIIn = true;
2987 0 : }
2988 0 : if (doRecluster) {
2989 0 : pTmaxMPIIn = muMI();
2990 0 : limitPTmaxMPIIn = true;
2991 0 : }
2992 : }
2993 :
2994 : // Done
2995 0 : return true;
2996 :
2997 0 : }
2998 :
2999 : //--------------------------------------------------------------------------
3000 :
3001 : // Function to return the value of the merging scale function in the current
3002 : // event.
3003 :
3004 : double MergingHooks::tmsNow( const Event& event ) {
3005 :
3006 : // Get merging scale in current event.
3007 : double tnow = 0.;
3008 : // Use KT/Durham merging scale definition.
3009 0 : if ( doKTMerging() || doMGMerging() )
3010 0 : tnow = kTms(event);
3011 : // Use Lund PT merging scale definition.
3012 0 : else if ( doPTLundMerging() )
3013 0 : tnow = rhoms(event, false);
3014 : // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition.
3015 0 : else if ( doCutBasedMerging() )
3016 0 : tnow = cutbasedms(event);
3017 : // Use NLO merging (Lund PT) merging scale definition.
3018 0 : else if ( doNL3Merging() )
3019 0 : tnow = rhoms(event, false);
3020 : // Use NLO merging (Lund PT) merging scale definition.
3021 0 : else if ( doUNLOPSMerging() )
3022 0 : tnow = rhoms(event, false);
3023 : // Use UMEPS (Lund PT) merging scale definition.
3024 0 : else if ( doUMEPSMerging() )
3025 0 : tnow = rhoms(event, false);
3026 : // Use user-defined merging scale.
3027 : else
3028 0 : tnow = tmsDefinition(event);
3029 : // Return merging scale value. Done
3030 0 : return tnow;
3031 : }
3032 :
3033 : //--------------------------------------------------------------------------
3034 :
3035 : // Function to check if the properties of the input particle should be
3036 : // checked against the cut-based merging scale defintion.
3037 :
3038 : bool MergingHooks::checkAgainstCut( const Particle& particle){
3039 :
3040 : // Do not check uncoloured particles.
3041 0 : if (particle.colType() == 0) return false;
3042 : // By default, use u-, d-, c-, s- and b-quarks.
3043 0 : if ( particle.idAbs() != 21 && particle.idAbs() > nQuarksMergeSave )
3044 0 : return false;
3045 : // Done
3046 0 : return true;
3047 :
3048 0 : }
3049 :
3050 : //--------------------------------------------------------------------------
3051 :
3052 : // Function to return the minimal kT in the event. If doKTMerging = true, this
3053 : // function will be used as a merging scale definition.
3054 :
3055 : double MergingHooks::kTms(const Event& event) {
3056 :
3057 : // Only check first emission.
3058 0 : if (!isFirstEmission(event)) return 0.;
3059 :
3060 : // Find all electroweak decayed bosons in the state.
3061 0 : vector<int> ewResonancePos;
3062 0 : ewResonancePos.clear();
3063 0 : for (int i=0; i < event.size(); ++i)
3064 0 : if ( abs(event[i].status()) == 22
3065 0 : && ( event[i].idAbs() == 22
3066 0 : || event[i].idAbs() == 23
3067 0 : || event[i].idAbs() == 24
3068 0 : || event[i].idAbs() == 25
3069 0 : || event[i].idAbs() == 6 ) )
3070 0 : ewResonancePos.push_back(i);
3071 :
3072 : // Declare final parton vectors
3073 0 : vector <int> FinalPartPos;
3074 0 : FinalPartPos.clear();
3075 : // Search inEvent record for final state partons.
3076 : // Exclude decay products of ew resonance.
3077 0 : for (int i=0; i < event.size(); ++i){
3078 0 : if ( event[i].isFinal()
3079 0 : && isInHard( i, event )
3080 0 : && event[i].colType() != 0
3081 0 : && checkAgainstCut(event[i]) ){
3082 : bool isDecayProduct = false;
3083 0 : for(int j=0; j < int(ewResonancePos.size()); ++j)
3084 0 : if ( event[i].isAncestor( ewResonancePos[j]) )
3085 0 : isDecayProduct = true;
3086 : // Except for e+e- -> jets, do not check radiation in resonance decays.
3087 0 : if ( !isDecayProduct
3088 0 : || getProcessString().compare("e+e->jj") == 0
3089 0 : || getProcessString().compare("e+e->(z>jj)") == 0 )
3090 0 : FinalPartPos.push_back(i);
3091 0 : }
3092 : }
3093 :
3094 : // Find minimal Durham kT in event, using own function: Check
3095 : // definition of separation
3096 0 : int type = (event[3].colType() == 0
3097 0 : && event[4].colType() == 0) ? -1 : ktTypeSave;
3098 : // Find minimal kT
3099 0 : double ktmin = event[0].e();
3100 0 : for(int i=0; i < int(FinalPartPos.size()); ++i){
3101 0 : double kt12 = ktmin;
3102 : // Compute separation to the beam axis for hadronic collisions
3103 0 : if (type == 1 || type == 2) {
3104 0 : double temp = event[FinalPartPos[i]].pT();
3105 0 : kt12 = min(kt12, temp);
3106 0 : }
3107 : // Compute separation to other final state jets
3108 0 : for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
3109 0 : double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
3110 0 : type, DparameterSave);
3111 0 : kt12 = min(kt12, temp);
3112 0 : }
3113 : // Keep the minimal Durham separation
3114 0 : ktmin = min(ktmin,kt12);
3115 0 : }
3116 :
3117 : // Return minimal Durham kT
3118 0 : return ktmin;
3119 :
3120 0 : }
3121 :
3122 : //--------------------------------------------------------------------------
3123 :
3124 : // Function to compute durham y separation from Particle input.
3125 :
3126 : double MergingHooks::kTdurham(const Particle& RadAfterBranch,
3127 : const Particle& EmtAfterBranch, int Type, double D ){
3128 :
3129 : // Declare return variable
3130 : double ktdur;
3131 : // Save 4-momenta of final state particles
3132 0 : Vec4 jet1 = RadAfterBranch.p();
3133 0 : Vec4 jet2 = EmtAfterBranch.p();
3134 :
3135 0 : if ( Type == -1) {
3136 : // Get angle between jets for e+e- collisions, make sure that
3137 : // -1 <= cos(theta) <= 1
3138 : double costh;
3139 0 : if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
3140 : else {
3141 0 : costh = costheta(jet1,jet2);
3142 : }
3143 : // Calculate kt durham separation between jets for e+e- collisions
3144 0 : ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
3145 0 : } else if ( Type == 1 ){
3146 : // Get delta_y for hadronic collisions:
3147 : // Get mT of first jet
3148 0 : double mT1sq = jet1.m2Calc() + jet1.pT2();
3149 : double mT1 = 0.;
3150 0 : if (mT1sq < 0) mT1 = - sqrt(-mT1sq);
3151 0 : else mT1 = sqrt(mT1sq);
3152 : // Get mT of second jet
3153 0 : double mT2sq = jet2.m2Calc() + jet2.pT2();
3154 : double mT2 = 0.;
3155 0 : if (mT2sq < 0) mT2 = - sqrt(-mT2sq);
3156 0 : else mT2 = sqrt(mT2sq);
3157 : // Get rapidity of first jet
3158 0 : double y1 = log( ( jet1.e() + abs(jet1.pz()) ) / mT1 );
3159 0 : if (jet1.pz() < 0) y1 *= -1.;
3160 : // Get rapidity of second jet
3161 0 : double y2 = log( ( jet2.e() + abs(jet2.pz()) ) / mT2 );
3162 0 : if (jet2.pz() < 0) y2 *= -1.;
3163 : // Get delta_phi for hadronic collisions
3164 0 : double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
3165 0 : double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
3166 0 : double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
3167 0 : double dPhi = acos( cosdPhi );
3168 : // Calculate kT durham like fastjet,
3169 : // but with rapidity instead of pseudo-rapidity
3170 0 : ktdur = min( pow(pt1,2),pow(pt2,2) )
3171 0 : * ( pow(y1-y2,2) + pow(dPhi,2) ) / pow(D,2);
3172 0 : } else if ( Type == 2 ){
3173 :
3174 : // Get mT of first jet
3175 0 : double mT1sq = jet1.m2Calc() + jet1.pT2();
3176 : double mT1 = 0.;
3177 0 : if (mT1sq < 0) mT1 = - sqrt(-mT1sq);
3178 0 : else mT1 = sqrt(mT1sq);
3179 : // Get mT of second jet
3180 0 : double mT2sq = jet2.m2Calc() + jet2.pT2();
3181 : double mT2 = 0.;
3182 0 : if (mT2sq < 0) mT2 = - sqrt(-mT2sq);
3183 0 : else mT2 = sqrt(mT2sq);
3184 : // Get pseudo-rapidity of first jet
3185 0 : double eta1 = log( ( sqrt(jet1.px()*jet1.px() + jet1.py()*jet1.py()
3186 0 : + jet1.pz()*jet1.pz()) + abs(jet1.pz()) ) / mT1);
3187 0 : if (jet1.pz() < 0) eta1 *= -1.;
3188 : // Get pseudo-rapidity of second jet
3189 0 : double eta2 = log( ( sqrt(jet2.px()*jet2.px() + jet2.py()*jet2.py()
3190 0 : + jet2.pz()*jet2.pz()) + abs(jet2.pz()) ) / mT2);
3191 0 : if (jet2.pz() < 0) eta2 *= -1.;
3192 :
3193 : // Get delta_phi and cos(Delta_phi) for hadronic collisions
3194 0 : double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
3195 0 : double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
3196 0 : double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
3197 0 : double dPhi = acos( cosdPhi );
3198 : // Calculate kT durham like fastjet
3199 0 : ktdur = min( pow(pt1,2),pow(pt2,2) )
3200 0 : * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
3201 0 : } else if ( Type == 3 ){
3202 : // Get cosh(Delta_eta) for hadronic collisions
3203 0 : double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
3204 0 : double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
3205 0 : double coshdEta = cosh( eta1 - eta2 );
3206 : // Get delta_phi and cos(Delta_phi) for hadronic collisions
3207 0 : double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
3208 0 : double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
3209 0 : double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
3210 : // Calculate kT durham separation "SHERPA-like"
3211 0 : ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
3212 0 : * ( coshdEta - cosdPhi ) / pow(D,2);
3213 0 : } else {
3214 : ktdur = 0.0;
3215 : }
3216 : // Return kT
3217 0 : return sqrt(ktdur);
3218 0 : }
3219 :
3220 : //--------------------------------------------------------------------------
3221 :
3222 : // Find the minimal Lund pT between coloured partons in the input
3223 : // event. If doPTLundMerging = true, this function will be used as a merging
3224 : // scale definition.
3225 :
3226 : double MergingHooks::rhoms( const Event& event, bool withColour){
3227 :
3228 : // Only check first emission.
3229 0 : if (!isFirstEmission(event)) return 0.;
3230 :
3231 : // Find all electroweak decayed bosons in the state.
3232 0 : vector<int> ewResonancePos;
3233 0 : ewResonancePos.clear();
3234 0 : for (int i=0; i < event.size(); ++i)
3235 0 : if ( abs(event[i].status()) == 22
3236 0 : && ( event[i].idAbs() == 22
3237 0 : || event[i].idAbs() == 23
3238 0 : || event[i].idAbs() == 24
3239 0 : || event[i].idAbs() == 25
3240 0 : || event[i].idAbs() == 6 ) )
3241 0 : ewResonancePos.push_back(i);
3242 :
3243 : // Declare final parton vectors
3244 0 : vector <int> FinalPartPos;
3245 0 : FinalPartPos.clear();
3246 : // Search inEvent record for final state partons.
3247 : // Exclude decay products of ew resonance.
3248 0 : for (int i=0; i < event.size(); ++i){
3249 :
3250 0 : if ( event[i].isFinal()
3251 0 : && isInHard( i, event )
3252 0 : && event[i].colType() != 0
3253 0 : && checkAgainstCut(event[i]) ){
3254 : bool isDecayProduct = false;
3255 0 : for(int j=0; j < int(ewResonancePos.size()); ++j)
3256 0 : if ( event[i].isAncestor( ewResonancePos[j]) )
3257 0 : isDecayProduct = true;
3258 : // Except for e+e- -> jets, do not check radiation in resonance decays.
3259 0 : if ( !isDecayProduct
3260 0 : || getProcessString().compare("e+e->jj") == 0
3261 0 : || getProcessString().compare("e+e->(z>jj)") == 0 )
3262 0 : FinalPartPos.push_back(i);
3263 0 : }
3264 : }
3265 :
3266 : // Get index of first incoming
3267 : int in1 = 0;
3268 0 : for (int i=0; i < event.size(); ++i)
3269 0 : if (abs(event[i].status()) == 41 ){
3270 : in1 = i;
3271 0 : break;
3272 : }
3273 :
3274 : // Get index of second incoming
3275 : int in2 = 0;
3276 0 : for (int i=0; i < event.size(); ++i)
3277 0 : if (abs(event[i].status()) == 42 ){
3278 : in2 = i;
3279 0 : break;
3280 : }
3281 :
3282 : // If no incoming of the cascade are found, try incoming
3283 0 : if (in1 == 0 || in2 == 0){
3284 : // Find current incoming partons
3285 0 : for(int i=3; i < int(event.size()); ++i){
3286 0 : if ( !isInHard( i, event ) ) continue;
3287 0 : if (event[i].mother1() == 1) in1 = i;
3288 0 : if (event[i].mother1() == 2) in2 = i;
3289 : }
3290 0 : }
3291 :
3292 : // Find minimal pythia pt in event
3293 0 : double ptmin = event[0].e();
3294 0 : for(int i=0; i < int(FinalPartPos.size()); ++i){
3295 0 : double pt12 = ptmin;
3296 : // Compute pythia ISR separation i-jet and first incoming
3297 0 : if (event[in1].colType() != 0) {
3298 0 : double temp = rhoPythia( event[in1],
3299 0 : event[FinalPartPos[i]], event[in2], -1 );
3300 0 : pt12 = min(pt12, temp);
3301 0 : }
3302 :
3303 : // Compute pythia ISR separation i-jet and second incoming
3304 0 : if ( event[in2].colType() != 0) {
3305 0 : double temp = rhoPythia( event[in2],
3306 0 : event[FinalPartPos[i]], event[in1], -1 );
3307 0 : pt12 = min(pt12, temp);
3308 0 : }
3309 :
3310 0 : if (withColour) {
3311 : // Compute pythia FSR separation between two jets,
3312 : // with knowledge of colour connections
3313 0 : for(int j=0; j < int(FinalPartPos.size()); ++j) {
3314 :
3315 : // Find recoiler in event
3316 0 : if ( i!=j ){
3317 : bool isHard = false;
3318 0 : int radAcl = event[FinalPartPos[i]].acol();
3319 0 : int radCol = event[FinalPartPos[i]].col();
3320 0 : int emtAcl = event[FinalPartPos[j]].acol();
3321 0 : int emtCol = event[FinalPartPos[j]].col();
3322 : int iRec = -1;
3323 : // Check in final state
3324 0 : if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
3325 0 : iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
3326 : event,1, isHard);
3327 0 : if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
3328 0 : iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
3329 : event,1, isHard);
3330 0 : if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
3331 0 : iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
3332 : event,1, isHard);
3333 0 : if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
3334 0 : iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
3335 : event,1, isHard);
3336 :
3337 : // Check in initial state
3338 0 : if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
3339 0 : iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
3340 : event,2, isHard);
3341 0 : if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
3342 0 : iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
3343 : event,2, isHard);
3344 0 : if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
3345 0 : iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
3346 : event,2, isHard);
3347 0 : if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
3348 0 : iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
3349 : event,2, isHard);
3350 :
3351 0 : if (iRec > 0) {
3352 0 : double temp = rhoPythia( event[FinalPartPos[i]],
3353 0 : event[FinalPartPos[j]],
3354 0 : event[iRec], 1 );
3355 0 : pt12 = min(pt12, temp);
3356 0 : }
3357 0 : }
3358 : }
3359 :
3360 0 : } else {
3361 : // Compute pythia FSR separation between two jets,
3362 : // without any knowledge of colour connections
3363 0 : for(int j=0; j < int(FinalPartPos.size()); ++j) {
3364 0 : for(int k=0; k < int(FinalPartPos.size()); ++k) {
3365 : // Allow any parton as recoiler
3366 0 : if ( (i != j) && (i != k) && (j != k) ){
3367 :
3368 0 : double temp = 0.;
3369 : // Only check splittings allowed by flavour, e.g.
3370 : // only q -> qg and g -> qqbar
3371 0 : temp = rhoPythia( event[FinalPartPos[i]],
3372 0 : event[FinalPartPos[j]],
3373 0 : event[FinalPartPos[k]], 1 );
3374 0 : pt12 = min(pt12, temp);
3375 0 : }
3376 : }
3377 : }
3378 :
3379 : // Compute pythia FSR separation between two jets, with initial recoiler
3380 : // without any knowledge of colour connections
3381 0 : if ( event[in1].colType() != 0 && event[in2].colType() != 0) {
3382 0 : for(int j=0; j < int(FinalPartPos.size()); ++j) {
3383 : // Allow both initial partons as recoiler
3384 0 : if ( i != j ){
3385 : // Check with first initial as recoiler
3386 0 : double temp = rhoPythia( event[FinalPartPos[i]],
3387 0 : event[FinalPartPos[j]],
3388 0 : event[in1], 1 );
3389 0 : pt12 = min(pt12, temp);
3390 : // Check with second initial as recoiler
3391 0 : temp = rhoPythia( event[FinalPartPos[i]],
3392 0 : event[FinalPartPos[j]],
3393 0 : event[in2], 1 );
3394 0 : pt12 = min(pt12, temp);
3395 0 : }
3396 : }
3397 0 : }
3398 :
3399 : }
3400 : // Reset minimal y separation
3401 0 : ptmin = min(ptmin,pt12);
3402 0 : }
3403 :
3404 0 : return ptmin;
3405 :
3406 0 : }
3407 :
3408 : //--------------------------------------------------------------------------
3409 :
3410 : // Function to compute "pythia pT separation" from Particle input, as a helper
3411 : // for rhoms(...).
3412 :
3413 : double MergingHooks::rhoPythia(const Particle& RadAfterBranch,
3414 : const Particle& EmtAfterBranch,
3415 : const Particle& RecAfterBranch, int ShowerType){
3416 :
3417 : // Save type: 1 = FSR pT definition, else ISR definition
3418 : int Type = ShowerType;
3419 : // Calculate virtuality of splitting
3420 0 : int sign = (Type==1) ? 1 : -1;
3421 0 : Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p());
3422 0 : double Qsq = sign * Q.m2Calc();
3423 : // Mass term of radiator
3424 0 : double m2Rad = ( includeMassive()
3425 0 : && abs(RadAfterBranch.id()) >= 4
3426 0 : && abs(RadAfterBranch.id()) < 7)
3427 0 : ? pow(particleDataPtr->m0(RadAfterBranch.id()), 2)
3428 : : 0.;
3429 : // Construct 2->3 variables for FSR
3430 0 : Vec4 sum = RadAfterBranch.p() + RecAfterBranch.p()
3431 0 : + EmtAfterBranch.p();
3432 0 : double m2Dip = sum.m2Calc();
3433 0 : double x1 = 2. * (sum * RadAfterBranch.p()) / m2Dip;
3434 0 : double x3 = 2. * (sum * EmtAfterBranch.p()) / m2Dip;
3435 : // Construct momenta of dipole before/after splitting for ISR
3436 0 : Vec4 qBR(RadAfterBranch.p() - EmtAfterBranch.p() + RecAfterBranch.p());
3437 0 : Vec4 qAR(RadAfterBranch.p() + RecAfterBranch.p());
3438 : // Calculate z of splitting, different for FSR and ISR
3439 0 : double z = (Type==1) ? x1/(x1+x3)
3440 0 : : (qBR.m2Calc())/( qAR.m2Calc());
3441 : // Separation of splitting, different for FSR and ISR
3442 0 : double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
3443 : // pT^2 = separation*virtuality
3444 0 : pTpyth *= (Qsq - sign*m2Rad);
3445 0 : if (pTpyth < 0.) pTpyth = 0.;
3446 : // Return pT
3447 0 : return sqrt(pTpyth);
3448 0 : }
3449 :
3450 : //--------------------------------------------------------------------------
3451 :
3452 : // Function to find a colour (anticolour) index in the input event.
3453 : // Helper for rhoms
3454 : // IN int col : Colour tag to be investigated
3455 : // int iExclude1 : Identifier of first particle to be excluded
3456 : // from search
3457 : // int iExclude2 : Identifier of second particle to be excluded
3458 : // from search
3459 : // Event event : event to be searched for colour tag
3460 : // int type : Tag to define if col should be counted as
3461 : // colour (type = 1) [->find anti-colour index
3462 : // contracted with col]
3463 : // anticolour (type = 2) [->find colour index
3464 : // contracted with col]
3465 : // OUT int : Position of particle in event record
3466 : // contraced with col [0 if col is free tag]
3467 :
3468 : int MergingHooks::findColour(int col, int iExclude1, int iExclude2,
3469 : const Event& event, int type, bool isHardIn){
3470 :
3471 0 : bool isHard = isHardIn;
3472 : int index = 0;
3473 :
3474 0 : if (isHard){
3475 : // Search event record for matching colour & anticolour
3476 0 : for(int n = 0; n < event.size(); ++n) {
3477 0 : if ( n != iExclude1 && n != iExclude2
3478 0 : && event[n].colType() != 0
3479 0 : &&( event[n].status() > 0 // Check outgoing
3480 0 : || event[n].status() == -21) ) { // Check incoming
3481 0 : if ( event[n].acol() == col ) {
3482 0 : index = -n;
3483 0 : break;
3484 : }
3485 0 : if ( event[n].col() == col ){
3486 : index = n;
3487 0 : break;
3488 : }
3489 : }
3490 : }
3491 0 : } else {
3492 :
3493 : // Search event record for matching colour & anticolour
3494 0 : for(int n = 0; n < event.size(); ++n) {
3495 0 : if ( n != iExclude1 && n != iExclude2
3496 0 : && event[n].colType() != 0
3497 0 : &&( event[n].status() == 43 // Check outgoing from ISR
3498 0 : || event[n].status() == 51 // Check outgoing from FSR
3499 0 : || event[n].status() == 52 // Check outgoing from FSR
3500 0 : || event[n].status() == -41 // first initial
3501 0 : || event[n].status() == -42) ) { // second initial
3502 0 : if ( event[n].acol() == col ) {
3503 0 : index = -n;
3504 0 : break;
3505 : }
3506 0 : if ( event[n].col() == col ){
3507 : index = n;
3508 0 : break;
3509 : }
3510 : }
3511 : }
3512 : }
3513 : // if no matching colour / anticolour has been found, return false
3514 0 : if ( type == 1 && index < 0) return abs(index);
3515 0 : if ( type == 2 && index > 0) return abs(index);
3516 :
3517 0 : return 0;
3518 0 : }
3519 :
3520 : //--------------------------------------------------------------------------
3521 :
3522 : // Find the if the event passes the Delta R_{ij}, pT_{i} and Q_{ij} cuts on
3523 : // the matrix element, as needed for cut-based merging scale definition.
3524 :
3525 : double MergingHooks::cutbasedms( const Event& event ){
3526 :
3527 : // Only check first emission.
3528 0 : if (!isFirstEmission(event)) return -1.;
3529 :
3530 : // Save allowed final state particles
3531 0 : vector<int> partons;
3532 0 : for( int i=0; i < event.size(); ++i) {
3533 0 : if ( event[i].isFinal()
3534 0 : && isInHard( i, event )
3535 0 : && checkAgainstCut(event[i]) )
3536 0 : partons.push_back(i);
3537 : }
3538 :
3539 : // Declare overall veto
3540 : bool doVeto = false;
3541 : // Declare vetoes
3542 : bool vetoPT = false;
3543 : bool vetoRjj = false;
3544 : bool vetoMjj = false;
3545 : // Declare cuts used in matrix element
3546 0 : double pTjmin = pTiMS();
3547 0 : double mjjmin = QijMS();
3548 0 : double rjjmin = dRijMS();
3549 :
3550 : // Declare minimum values
3551 0 : double minPT = event[0].e();
3552 0 : double minRJJ = 10.;
3553 0 : double minMJJ = event[0].e();
3554 :
3555 : // Check matrix element cuts
3556 0 : for( int i=0; i < int(partons.size()); ++i){
3557 : // Save pT value
3558 0 : minPT = min(minPT,event[partons[i]].pT());
3559 :
3560 : // Check two-parton cuts
3561 0 : for( int j=0; j < int(partons.size()); ++j){
3562 0 : if (i!=j){
3563 :
3564 : // Save delta R value
3565 0 : minRJJ = min(minRJJ, deltaRij( event[partons[i]].p(),
3566 0 : event[partons[j]].p()) );
3567 : // Save delta R value
3568 0 : minMJJ = min(minMJJ, ( event[partons[i]].p()
3569 0 : +event[partons[j]].p() ).mCalc() );
3570 :
3571 0 : }
3572 : }
3573 : // Done with cut evaluation
3574 : }
3575 :
3576 : // Check if all partons are in the matrix element region
3577 0 : if (minPT > pTjmin) vetoPT = true;
3578 0 : if (minRJJ > rjjmin) vetoRjj = true;
3579 0 : if (minMJJ > mjjmin) vetoMjj = true;
3580 :
3581 : // In the matrix element calculation, we have rejected the events if any of
3582 : // the cuts had not been fulfilled,
3583 : // i.e. we are in the matrix element domain if all cuts are fulfilled.
3584 : // We veto any emission in the ME region.
3585 : // Disregard the two-parton cuts if only one parton in the final state
3586 0 : if (int(partons.size() == 1))
3587 0 : doVeto = vetoPT;
3588 : else
3589 : // Veto if the combination of cuts would be in the ME region
3590 0 : doVeto = vetoPT && vetoRjj && vetoMjj;
3591 :
3592 : // If event is above merging scale, veto
3593 0 : if (doVeto) return 1.;
3594 :
3595 : // Else, do nothing
3596 0 : return -1.;
3597 :
3598 0 : }
3599 :
3600 : //--------------------------------------------------------------------------
3601 :
3602 : // Function to compute Delta R separation from 4-vector input.
3603 :
3604 : double MergingHooks::deltaRij(Vec4 jet1, Vec4 jet2){
3605 :
3606 : // Declare return variable
3607 : double deltaR = 0.;
3608 : // Get delta_eta and cosh(Delta_eta) for hadronic collisions
3609 0 : double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
3610 0 : double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
3611 : // Get delta_phi and cos(Delta_phi) for hadronic collisions
3612 0 : double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
3613 0 : double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );
3614 0 : double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
3615 0 : double dPhi = acos( cosdPhi );
3616 : // Calculate kT durham like fastjet
3617 0 : deltaR = sqrt(pow(eta1-eta2,2) + pow(dPhi,2));
3618 : // Return kT
3619 0 : return deltaR;
3620 : }
3621 :
3622 : //==========================================================================
3623 :
3624 : } // end namespace Pythia8
|
