StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireMerging.cc
1 // DireMerging.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Stefan Prestel, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) used for Dire merging.
7 
8 #include "Pythia8/DireMerging.h"
9 #include "Pythia8/DireSpace.h"
10 #include "Pythia8/DireTimes.h"
11 #include <ctime>
12 
13 namespace Pythia8 {
14 
15 //==========================================================================
16 
17 // The Merging class.
18 
19 //--------------------------------------------------------------------------
20 
21 
22 // Check colour/flavour correctness of state.
23 
24 bool validEvent( const Event& event ) {
25 
26  bool validColour = true;
27  bool validCharge = true;
28  bool validMomenta = true;
29  double mTolErr=1e-2;
30 
31  // Check charge sum in initial and final state
32  double initCharge = event[3].charge() + event[4].charge();
33  double finalCharge = 0.0;
34  for(int i = 0; i < event.size(); ++i)
35  if (event[i].isFinal()) finalCharge += event[i].charge();
36  if (abs(initCharge-finalCharge) > 1e-12) validCharge = false;
37 
38  // Check that overall pT is vanishing.
39  Vec4 pSum(0.,0.,0.,0.);
40  for ( int i = 0; i < event.size(); ++i) {
41  if ( event[i].status() == -21 ) pSum -= event[i].p();
42  if ( event[i].isFinal() ) pSum += event[i].p();
43  }
44  if ( abs(pSum.px()) > mTolErr || abs(pSum.py()) > mTolErr) {
45  validMomenta = false;
46  }
47 
48  if ( event[3].status() == -21
49  && (abs(event[3].px()) > mTolErr || abs(event[3].py()) > mTolErr)){
50  validMomenta = false;
51  }
52  if ( event[4].status() == -21
53  && (abs(event[4].px()) > mTolErr || abs(event[4].py()) > mTolErr)){
54  validMomenta = false;
55  }
56 
57  return (validColour && validCharge && validMomenta);
58 
59 }
60 
61 //--------------------------------------------------------------------------
62 
63 // Initialise Merging class
64 
65 void DireMerging::init(){
66  // Reset minimal tms value.
67  tmsNowMin = infoPtr->eCM();
68 
69  enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
70  doMOPS = settingsPtr->flag("Dire:doMOPS");
71  applyTMSCut = settingsPtr->flag("Merging:doXSectionEstimate");
72  doMerging = settingsPtr->flag("Dire:doMerging");
73  usePDF = settingsPtr->flag("ShowerPDF:usePDF");
74  allowReject = settingsPtr->flag("Merging:applyVeto");
75  doMECs = settingsPtr->flag("Dire:doMECs");
76  doMEM = settingsPtr->flag("Dire:doMEM");
77  doGenerateSubtractions
78  = settingsPtr->flag("Dire:doGenerateSubtractions");
79  doGenerateMergingWeights
80  = settingsPtr->flag("Dire:doGenerateMergingWeights");
81  doExitAfterMerging
82  = settingsPtr->flag("Dire:doExitAfterMerging");
83  allowIncompleteReal
84  = settingsPtr->flag("Merging:allowIncompleteHistoriesInReal");
85  nQuarksMerge = settingsPtr->mode("Merging:nQuarksMerge");
86 
87  first = true;
88 
89 }
90 
91 //--------------------------------------------------------------------------
92 
93 // Function to print information.
94 void DireMerging::statistics() {
95 
96  // Recall merging scale value.
97  double tmsval = mergingHooksPtr->tms();
98  bool printBanner = enforceCutOnLHE && tmsNowMin > TMSMISMATCH*tmsval
99  && tmsval > 0.;
100  // Reset minimal tms value.
101  tmsNowMin = infoPtr->eCM();
102 
103  if (doMOPS) printBanner = false;
104  if (doMEM) printBanner = false;
105  if (doMECs) printBanner = false;
106 
107  if (!printBanner) return;
108 
109  // Header.
110  cout << "\n *------- PYTHIA Matrix Element Merging Information ------"
111  << "-------------------------------------------------------*\n"
112  << " | "
113  << " |\n";
114  // Print warning if the minimal tms value of any event was significantly
115  // above the desired merging scale value.
116  cout << " | Warning in DireMerging::statistics: All Les Houches events"
117  << " significantly above Merging:TMS cut. Please check. |\n";
118 
119  // Listing finished.
120  cout << " | "
121  << " |\n"
122  << " *------- End PYTHIA Matrix Element Merging Information -----"
123  << "-----------------------------------------------------*" << endl;
124 }
125 
126 //--------------------------------------------------------------------------
127 
128 void DireMerging::storeInfos() {
129 
130  // Clear previous information.
131  clearInfos();
132 
133  // Store information on every possible last clustering.
134  for ( int i = 0 ; i < int(myHistory->children.size()); ++i) {
135  // Just store pT and mass for now.
136  stoppingScalesSave.push_back(myHistory->children[i]->clusterIn.pT());
137  radSave.push_back(myHistory->children[i]->clusterIn.radPos());
138  emtSave.push_back(myHistory->children[i]->clusterIn.emtPos());
139  recSave.push_back(myHistory->children[i]->clusterIn.recPos());
140  mDipSave.push_back(myHistory->children[i]->clusterIn.mass());
141  }
142 
143 }
144 
145 //--------------------------------------------------------------------------
146 
147 void DireMerging::getStoppingInfo(double scales [100][100],
148  double masses [100][100]) {
149 
150  int posOffest=2;
151  for (unsigned int i=0; i < radSave.size(); ++i){
152  scales[radSave[i]-posOffest][recSave[i]-posOffest] = stoppingScalesSave[i];
153  masses[radSave[i]-posOffest][recSave[i]-posOffest] = mDipSave[i];
154  }
155 
156 }
157 
158 double DireMerging::generateSingleSudakov ( double pTbegAll,
159  double pTendAll, double m2dip, int idA, int type, double s, double x) {
160  return isr->noEmissionProbability( pTbegAll, pTendAll, m2dip, idA,
161  type, s, x);
162 }
163 
164 //--------------------------------------------------------------------------
165 
166 void DireMerging::reset() {
167  partonSystemsPtr->clear();
168  isr->clear();
169  fsr->clear();
170  beamAPtr->clear();
171  beamBPtr->clear();
172 }
173 
174 // Function to steer different merging prescriptions.
175 
176 bool DireMerging::generateUnorderedPoint(Event& process){
177 
178  cout << "generate new unordered point" << endl;
179 
180  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
181 
182  bool found = false;
183 
184  // Dryrun to find overestimate enhancements if necessary.
185  if (first) {
186 
187  isr->dryrun = true;
188  fsr->dryrun = true;
189 
190  for (int iii=0; iii<50; ++iii) {
191 
192  int sizeOld = newProcess.size();
193  int sizeNew = newProcess.size();
194  int nSplit = 2;
195 
196  while (sizeNew < sizeOld+nSplit) {
197 
198  // Incoming partons to hard process are stored in slots 3 and 4.
199  int inS = 0;
200  int inP = 3;
201  int inM = 4;
202  int sizeProcess = newProcess.size();
203  int sizeEvent = newProcess.size();
204  int nOffset = sizeEvent - sizeProcess;
205  int nHardDone = sizeProcess;
206  double scale = newProcess.scale();
207  // Store participating partons as first set in list of all systems.
208  partonSystemsPtr->addSys();
209  partonSystemsPtr->setInA(0, inP + nOffset);
210  partonSystemsPtr->setInB(0, inM + nOffset);
211  for (int i = inM + 1; i < nHardDone; ++i)
212  partonSystemsPtr->addOut(0, i + nOffset);
213  partonSystemsPtr->setSHat( 0,
214  (newProcess[inP + nOffset].p() +
215  newProcess[inM + nOffset].p()).m2Calc() );
216  partonSystemsPtr->setPTHat( 0, scale);
217 
218  // Add incoming hard-scattering partons to list in beam remnants.
219  double x1 = newProcess[inP].pPos() / newProcess[inS].m();
220  double x2 = newProcess[inM].pNeg() / newProcess[inS].m();
221  beamAPtr->append( inP + nOffset, newProcess[inP].id(), x1);
222  beamBPtr->append( inM + nOffset, newProcess[inM].id(), x2);
223  // Scale. Find whether incoming partons are valence or sea. Store.
224  beamAPtr->xfISR( 0, newProcess[inP].id(), x1, scale*scale);
225  int vsc1 = beamAPtr->pickValSeaComp();
226  beamBPtr->xfISR( 0, newProcess[inM].id(), x2, scale*scale);
227  int vsc2 = beamBPtr->pickValSeaComp();
228  bool isVal1 = (vsc1 == -3);
229  bool isVal2 = (vsc2 == -3);
230  infoPtr->setValence( isVal1, isVal2);
231 
232  isr->prepare(0, newProcess, true);
233  fsr->prepare(0, newProcess, true);
234 
235  Event in(newProcess);
236  double pTnowFSR = fsr->newPoint(in);
237  double pTnowISR = isr->newPoint(in);
238  if (pTnowFSR==0. && pTnowISR==0.) { reset(); continue; }
239  bool branched=false;
240  if (pTnowFSR > pTnowISR) branched = fsr->branch(in);
241  else branched = isr->branch(in);
242  if (!branched) { reset(); continue; }
243  if (pTnowFSR > pTnowISR) isr->update(0, in, false);
244  else fsr->update(0, in, false);
245 
246  pTnowISR = isr->newPoint(in);
247  pTnowFSR = fsr->newPoint(in);
248  if (pTnowFSR==0. && pTnowISR==0.) { reset(); continue; }
249  branched=false;
250  if (pTnowFSR > pTnowISR) branched = fsr->branch(in);
251  else branched = isr->branch(in);
252  if (!branched) { reset(); continue; }
253  if (pTnowFSR > pTnowISR) isr->update(0, in, false);
254  else fsr->update(0, in, false);
255 
256  // Done. Clean up event and return.
257  Event bla = fsr->makeHardEvent(0, in, false);
258  sizeNew = bla.size();
259  partonSystemsPtr->clear();
260  isr->clear();
261  fsr->clear();
262  beamAPtr->clear();
263  beamBPtr->clear();
264  }
265  }
266 
267  // Generate arbitrary phase-space point with two additional particles.
268  } else {
269 
270  int sizeOld = newProcess.size();
271  int sizeNew = newProcess.size();
272  int nSplit = 2;
273 
274  while (sizeNew < sizeOld+nSplit) {
275 
276  // Incoming partons to hard process are stored in slots 3 and 4.
277  int inS = 0;
278  int inP = 3;
279  int inM = 4;
280  int sizeProcess = newProcess.size();
281  int sizeEvent = newProcess.size();
282  int nOffset = sizeEvent - sizeProcess;
283  int nHardDone = sizeProcess;
284  double scale = newProcess.scale();
285  // Store participating partons as first set in list of all systems.
286  partonSystemsPtr->addSys();
287  partonSystemsPtr->setInA(0, inP + nOffset);
288  partonSystemsPtr->setInB(0, inM + nOffset);
289  for (int i = inM + 1; i < nHardDone; ++i)
290  partonSystemsPtr->addOut(0, i + nOffset);
291  partonSystemsPtr->setSHat( 0,
292  (newProcess[inP + nOffset].p() +
293  newProcess[inM + nOffset].p()).m2Calc() );
294  partonSystemsPtr->setPTHat( 0, scale);
295 
296  // Add incoming hard-scattering partons to list in beam remnants.
297  double x1 = newProcess[inP].pPos() / newProcess[inS].m();
298  double x2 = newProcess[inM].pNeg() / newProcess[inS].m();
299  beamAPtr->append( inP + nOffset, newProcess[inP].id(), x1);
300  beamBPtr->append( inM + nOffset, newProcess[inM].id(), x2);
301 
302  // Scale. Find whether incoming partons are valence or sea. Store.
303  // When an x-dependent matter profile is used with nonDiffractive,
304  // trial interactions mean that the valence/sea choice has already
305  // been made and should be restored here.
306  beamAPtr->xfISR( 0, newProcess[inP].id(), x1, scale*scale);
307  int vsc1 = beamAPtr->pickValSeaComp();
308  beamBPtr->xfISR( 0, newProcess[inM].id(), x2, scale*scale);
309  int vsc2 = beamBPtr->pickValSeaComp();
310  bool isVal1 = (vsc1 == -3);
311  bool isVal2 = (vsc2 == -3);
312  infoPtr->setValence( isVal1, isVal2);
313 
314  isr->prepare(0, newProcess, true);
315  fsr->prepare(0, newProcess, true);
316 
317  Event in(newProcess);
318 
319  double pTnowFSR = fsr->newPoint(in);
320  double pTnowISR = isr->newPoint(in);
321 
322  if (pTnowFSR==0. && pTnowISR==0.) { reset(); continue; }
323  bool branched=false;
324  if (pTnowFSR > pTnowISR) branched = fsr->branch(in);
325  else branched = isr->branch(in);
326  if (!branched) { reset(); continue; }
327  if (pTnowFSR > pTnowISR) isr->update(0, in, false);
328  else fsr->update(0, in, false);
329 
330  pTnowISR = isr->newPoint(in);
331  pTnowFSR = fsr->newPoint(in);
332  if (pTnowFSR==0. && pTnowISR==0.) { reset(); continue; }
333  branched=false;
334  if (pTnowFSR > pTnowISR) branched = fsr->branch(in);
335  else branched = isr->branch(in);
336  if (!branched) { reset(); continue; }
337  if (pTnowFSR > pTnowISR) isr->update(0, in, false);
338  else fsr->update(0, in, false);
339 
340  // Done. Clean up event and return.
341  in = fsr->makeHardEvent(0, in, false);
342  reset();
343 
344  // Loop through event and count.
345  int nPartons(0);
346  for(int i=0; i < int(in.size()); ++i)
347  if ( in[i].isFinal()
348  && in[i].colType()!= 0
349  && ( in[i].id() == 21 || in[i].idAbs() <= nQuarksMerge))
350  nPartons++;
351  nPartons -= mergingHooksPtr->hardProcess->nQuarksOut();
352 
353  // Set number of requested partons.
354  settingsPtr->mode("Merging:nRequested", nPartons);
355  mergingHooksPtr->nRequestedSave
356  = settingsPtr->mode("Merging:nRequested");
357  mergingHooksPtr->reattachResonanceDecays(in);
358 
359  generateHistories(in, false);
360  reset();
361  if (myHistory->foundAnyOrderedPaths()) continue;
362  newProcess = in;
363  sizeNew = newProcess.size();
364  found = true;
365 
366  }
367  }
368 
369  if (!found) {
370  // Loop through event and count.
371  int nPartons(0);
372  for(int i=0; i < int(newProcess.size()); ++i)
373  if ( newProcess[i].isFinal()
374  && newProcess[i].colType()!= 0
375  && ( newProcess[i].id() == 21
376  || newProcess[i].idAbs() <= nQuarksMerge))
377  nPartons++;
378  nPartons -= mergingHooksPtr->hardProcess->nQuarksOut();
379 
380  // Set number of requested partons.
381  settingsPtr->mode("Merging:nRequested", nPartons);
382  mergingHooksPtr->nRequestedSave
383  = settingsPtr->mode("Merging:nRequested");
384  mergingHooksPtr->reattachResonanceDecays(newProcess);
385  }
386  process = newProcess;
387 
388  cout << "found unordered point" << endl;
389 
390  first=false;
391  isr->dryrun=false;
392  fsr->dryrun=false;
393 
394  return true;
395 
396 }
397 
398 //--------------------------------------------------------------------------
399 
400 // Function to steer different merging prescriptions.
401 
402 int DireMerging::mergeProcess(Event& process){
403 
404  // Clear all previous event-by-event information.
405  clearInfos();
406 
407  int vetoCode = 1;
408 
409  // Reinitialise hard process.
410  mergingHooksPtr->hardProcess->clear();
411  string processNow = settingsPtr->word("Merging:Process");
412  mergingHooksPtr->hardProcess->initOnProcess(processNow, particleDataPtr);
413 
414  // Remove whitespace from process string
415  while(processNow.find(" ", 0) != string::npos)
416  processNow.erase(processNow.begin()+processNow.find(" ",0));
417  mergingHooksPtr->processSave = processNow;
418 
419  mergingHooksPtr->doUserMergingSave
420  = settingsPtr->flag("Merging:doUserMerging");
421  mergingHooksPtr->doMGMergingSave
422  = settingsPtr->flag("Merging:doMGMerging");
423  mergingHooksPtr->doKTMergingSave
424  = settingsPtr->flag("Merging:doKTMerging");
425  mergingHooksPtr->doPTLundMergingSave
426  = settingsPtr->flag("Merging:doPTLundMerging");
427  mergingHooksPtr->doCutBasedMergingSave
428  = settingsPtr->flag("Merging:doCutBasedMerging");
429  mergingHooksPtr->doNL3TreeSave
430  = settingsPtr->flag("Merging:doNL3Tree");
431  mergingHooksPtr->doNL3LoopSave
432  = settingsPtr->flag("Merging:doNL3Loop");
433  mergingHooksPtr->doNL3SubtSave
434  = settingsPtr->flag("Merging:doNL3Subt");
435  mergingHooksPtr->doUNLOPSTreeSave
436  = settingsPtr->flag("Merging:doUNLOPSTree");
437  mergingHooksPtr->doUNLOPSLoopSave
438  = settingsPtr->flag("Merging:doUNLOPSLoop");
439  mergingHooksPtr->doUNLOPSSubtSave
440  = settingsPtr->flag("Merging:doUNLOPSSubt");
441  mergingHooksPtr->doUNLOPSSubtNLOSave
442  = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
443  mergingHooksPtr->doUMEPSTreeSave
444  = settingsPtr->flag("Merging:doUMEPSTree");
445  mergingHooksPtr->doUMEPSSubtSave
446  = settingsPtr->flag("Merging:doUMEPSSubt");
447  mergingHooksPtr->nReclusterSave
448  = settingsPtr->mode("Merging:nRecluster");
449 
450  mergingHooksPtr->hasJetMaxLocal = false;
451  mergingHooksPtr->nJetMaxLocal
452  = mergingHooksPtr->nJetMaxSave;
453  mergingHooksPtr->nJetMaxNLOLocal
454  = mergingHooksPtr->nJetMaxNLOSave;
455  mergingHooksPtr->nRequestedSave
456  = settingsPtr->mode("Merging:nRequested");
457 
458  // Reset to default merging scale.
459  mergingHooksPtr->tms(mergingHooksPtr->tmsCut());
460 
461  // Ensure that merging weight is not counted twice.
462  bool includeWGT = mergingHooksPtr->includeWGTinXSEC();
463 
464  // Directly retrive Sudakov w/o PDF factors from showers and exit.
465  //if (doMcAtNloDelta && !usePDF) return genSud(process);
466 
467  // Possibility to apply merging scale to an input event.
468  if ( applyTMSCut && cutOnProcess(process) ) {
469  if (includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
470  return -1;
471  }
472 
473  // Done if only a cut should be applied.
474  if ( applyTMSCut ) return 1;
475 
476  if (doMerging){
477 
478  int nPartons = 0;
479 
480  // Do not include resonance decay products in the counting.
481  Event newp( mergingHooksPtr->bareEvent( process, false) );
482  // Loop through event and count.
483  for(int i=0; i < int(newp.size()); ++i)
484  if ( newp[i].isFinal()
485  && newp[i].colType()!= 0
486  && ( newp[i].id() == 21 || newp[i].idAbs() <= nQuarksMerge))
487  nPartons++;
488 
489  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newp, false);
490 
491  nPartons -= mergingHooksPtr->hardProcess->nQuarksOut();
492 
493  // Set number of requested partons.
494  settingsPtr->mode("Merging:nRequested", nPartons);
495 
496  mergingHooksPtr->hasJetMaxLocal = false;
497  mergingHooksPtr->nJetMaxLocal
498  = mergingHooksPtr->nJetMaxSave;
499  mergingHooksPtr->nJetMaxNLOLocal
500  = mergingHooksPtr->nJetMaxNLOSave;
501  mergingHooksPtr->nRequestedSave
502  = settingsPtr->mode("Merging:nRequested");
503 
504  if (doMEM) {
505  int nFinal(0), nQuarks(0), nGammas(0);
506  for (int i=0; i < newp.size(); ++i) {
507  if (newp[i].idAbs() < 7) nQuarks++;
508  if (newp[i].idAbs() == 22) nGammas++;
509  if (newp[i].isFinal()) nFinal++;
510  }
511  settingsPtr->mode("DireSpace:nFinalMax",nFinal-1);
512  settingsPtr->mode("DireTimes:nFinalMax",nFinal-1);
513  if (nQuarks > 4) return 1;
514  }
515 
516  // Reset to default merging scale.
517  mergingHooksPtr->tms(mergingHooksPtr->tmsCut());
518 
519  // For ME corrections, only do mergingHooksPtr reinitialization here,
520  // and do not perform any veto.
521  if (doMECs) return 1;
522 
523  if (doMEM) mergingHooksPtr->orderHistories(false);
524 
525  clock_t begin = clock();
526 
527  if (doMOPS) generateUnorderedPoint(process);
528 
529  bool foundHistories = generateHistories(process, false);
530  int returnCode = (foundHistories) ? 1 : 0;
531 
532  if (doMOPS && myHistory->foundAnyOrderedPaths() && nSteps > 0)
533  returnCode = 0;
534 
535  std::clock_t end = std::clock();
536  double elapsed_secs_1 = double(end - begin) / CLOCKS_PER_SEC;
537  sum_paths += myHistory->goodBranches.size();
538  sum_time_1 += elapsed_secs_1;
539 
540  // For interative resummed matrix element method, tag histories and exit.
541  if (doMEM) {
542  tagHistories();
543  return 1;
544  }
545 
546  if (doGenerateSubtractions) calculateSubtractions();
547  bool useAll = doMOPS;
548 
549  double RNpath = getPathIndex(useAll);
550  if ((doMOPS && returnCode > 0) || doGenerateMergingWeights)
551  returnCode = calculateWeights(RNpath, useAll);
552 
553  end = std::clock();
554  double elapsed_secs_2 = double(end - begin) / CLOCKS_PER_SEC;
555  sum_time_2 += elapsed_secs_2;
556 
557  int tmp_code = getStartingConditions( RNpath, process);
558  if (returnCode > 0) returnCode = tmp_code;
559 
560  // Ensure that merging weight is not counted twice.
561  if (returnCode == 0) {
562  mergingHooksPtr->setWeightCKKWL({0.});
563  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
564  }
565 
566  if (!allowReject && returnCode < 1) returnCode=1;
567 
568  // Store information before leaving.
569  if (foundHistories) storeInfos();
570 
571  if (doMOPS) {
572  if (returnCode < 1) mergingHooksPtr->setWeightCKKWL({0.});
573  return returnCode;
574  }
575 
576  // Veto if we do not want to do event generation.
577  if (doExitAfterMerging) return -1;
578 
579  return 1;
580  }
581 
582  // Possibility to perform CKKW-L merging on this event.
583  if ( mergingHooksPtr->doCKKWLMerging() )
584  vetoCode = mergeProcessCKKWL(process);
585 
586  // Possibility to perform UMEPS merging on this event.
587  if ( mergingHooksPtr->doUMEPSMerging() )
588  vetoCode = mergeProcessUMEPS(process);
589 
590  // Possibility to perform NL3 NLO merging on this event.
591  if ( mergingHooksPtr->doNL3Merging() )
592  vetoCode = mergeProcessNL3(process);
593 
594  // Possibility to perform UNLOPS merging on this event.
595  if ( mergingHooksPtr->doUNLOPSMerging() )
596  vetoCode = mergeProcessUNLOPS(process);
597 
598  return vetoCode;
599 
600 }
601 
602 //--------------------------------------------------------------------------
603 
604 // Function to perform CKKW-L merging on this event.
605 
606 int DireMerging::mergeProcessCKKWL( Event& process) {
607 
608  // Ensure that merging hooks to not veto events in the trial showers.
609  mergingHooksPtr->doIgnoreStep(true);
610  // For pp > h, allow cut on state, so that underlying processes
611  // can be clustered to gg > h
612  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0 )
613  mergingHooksPtr->allowCutOnRecState(true);
614 
615  // Construct all histories.
616  // This needs to be done because MECs can depend on unordered paths if
617  // these unordered paths are ordered up to some point.
618  mergingHooksPtr->orderHistories(false);
619 
620  // Ensure that merging weight is not counted twice.
621  bool includeWGT = mergingHooksPtr->includeWGTinXSEC();
622 
623  // Reset weight of the event.
624  double wgt = 1.0;
625  mergingHooksPtr->setWeightCKKWL({1.});
626  mergingHooksPtr->muMI(-1.);
627 
628  // Prepare process record for merging. If Pythia has already decayed
629  // resonances used to define the hard process, remove resonance decay
630  // products.
631  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
632  // Reset any incoming spins for W+-.
633  if (mergingHooksPtr->doWeakClustering())
634  for (int i = 0;i < newProcess.size();++i)
635  newProcess[i].pol(9);
636  // Store candidates for the splitting V -> qqbar'.
637  mergingHooksPtr->storeHardProcessCandidates( newProcess);
638 
639  // Check if event passes the merging scale cut.
640  // Get merging scale in current event.
641  // Calculate number of clustering steps.
642  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
643 
644  double tmsnow = mergingHooksPtr->tmsNow( newProcess );
645 
646  // Too few steps can be possible if a chain of resonance decays has been
647  // removed. In this case, reject this event, since it will be handled in
648  // lower-multiplicity samples.
649  int nRequested = mergingHooksPtr->nRequested();
650 
651  // Store hard event cut information, reset veto information.
652  mergingHooksPtr->setHardProcessInfo(nSteps, tmsnow);
653  mergingHooksPtr->setEventVetoInfo(-1, -1.);
654 
655  if (nSteps < nRequested && allowReject) {
656  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({0.});
657  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
658  return -1;
659  }
660 
661  // Reset the minimal tms value, if necessary.
662  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
663 
664  // Set dummy process scale.
665  newProcess.scale(0.0);
666  // Generate all histories.
667  DireHistory FullHistory( nSteps, 0.0, newProcess, DireClustering(),
668  mergingHooksPtr, (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
669  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
670  1.0, 1.0, 1.0, 1.0, 0);
671 
672  // Project histories onto desired branches, e.g. only ordered paths.
673  FullHistory.projectOntoDesiredHistories();
674 
675  // Setup to choose shower starting conditions randomly.
676  double sumAll(0.), sumFullAll(0.);
677  for ( map<double, DireHistory*>::iterator
678  it = FullHistory.goodBranches.begin();
679  it != FullHistory.goodBranches.end(); ++it ) {
680  sumAll += it->second->prodOfProbs;
681  sumFullAll += it->second->prodOfProbsFull;
682  }
683  // Store a double with which to access each of the paths.
684  double lastp(0.);
685  vector<double> path_index;
686  for ( map<double, DireHistory*>::iterator
687  it = FullHistory.goodBranches.begin();
688  it != FullHistory.goodBranches.end(); ++it ) {
689  // Double to access path.
690  double indexNow = (lastp + 0.5*(it->first - lastp))/sumAll;
691  path_index.push_back(indexNow);
692  lastp = it->first;
693  }
694  // Randomly pick path.
695  int sizeBranches = FullHistory.goodBranches.size();
696  int iPosRN = (sizeBranches > 0)
697  ? rndmPtr->pick(
698  vector<double>(sizeBranches, 1./double(sizeBranches)) )
699  : 0;
700  double RN = (sizeBranches > 0) ? path_index[iPosRN] : rndmPtr->flat();
701 
702  // Setup the selected path. Needed for
703  FullHistory.select(RN)->setSelectedChild();
704 
705  // Do not apply cut if the configuration could not be projected onto an
706  // underlying born configuration.
707  bool applyCut = allowReject
708  && nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
709 
710  Event core( FullHistory.lowestMultProc(RN) );
711  // Set event-specific merging scale cut. Q2-dependent for DIS.
712  if ( mergingHooksPtr->getProcessString().compare("e+p>e+j") == 0
713  || mergingHooksPtr->getProcessString().compare("e-p>e-j") == 0 ) {
714 
715  // Set dynamical merging scale for DIS
716  if (FullHistory.isDIS2to2(core)) {
717  int iInEl(0), iOutEl(0);
718  for ( int i=0; i < core.size(); ++i )
719  if ( core[i].idAbs() == 11 ) {
720  if ( core[i].status() == -21 ) iInEl = i;
721  if ( core[i].isFinal() ) iOutEl = i;
722  }
723  double Q = sqrt( -(core[iInEl].p() - core[iOutEl].p() ).m2Calc());
724  double tmsCut = mergingHooksPtr->tmsCut();
725  double tmsEvt = tmsCut / sqrt( 1. + pow( tmsCut/ ( 0.5*Q ), 2) );
726  mergingHooksPtr->tms(tmsEvt);
727 
728  } else if (FullHistory.isMassless2to2(core)) {
729  double mT(1.);
730  for ( int i=0; i < core.size(); ++i )
731  if ( core[i].isFinal() ) mT *= core[i].mT();
732  double Q = sqrt(mT);
733  double tmsCut = mergingHooksPtr->tmsCut();
734  double tmsEvt = tmsCut / sqrt( 1. + pow( tmsCut/ ( 0.5*Q ), 2) );
735  mergingHooksPtr->tms(tmsEvt);
736  }
737  }
738  double tmsval = mergingHooksPtr->tms();
739 
740  // Enfore merging scale cut if the event did not pass the merging scale
741  // criterion.
742  if ( enforceCutOnLHE && applyCut && tmsnow < tmsval ) {
743  string message="Warning in DireMerging::mergeProcessCKKWL: "
744  "Les Houches Event";
745  message+=" fails merging scale cut. Reject event.";
746  infoPtr->errorMsg(message);
747  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({0.});
748  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
749  return -1;
750  }
751 
752  // Check if more steps should be taken.
753  int nFinalP(0), nFinalW(0), nFinalZ(0);
754  for ( int i = 0; i < core.size(); ++i )
755  if ( core[i].isFinal() ) {
756  if ( core[i].colType() != 0 ) nFinalP++;
757  if ( core[i].idAbs() == 24 ) nFinalW++;
758  if ( core[i].idAbs() == 23 ) nFinalZ++;
759  }
760  bool complete = (FullHistory.select(RN)->nClusterings() == nSteps) ||
761  ( mergingHooksPtr->doWeakClustering() && nFinalP == 2
762  && nFinalW+nFinalZ == 0);
763  if ( !complete ) {
764  string message="Warning in DireMerging::mergeProcessCKKWL: No clusterings";
765  message+=" found. History incomplete.";
766  infoPtr->errorMsg(message);
767  }
768 
769  // Calculate CKKWL reweighting for all paths.
770  double wgtsum(0.);
771  lastp = 0.;
772  for ( map<double, DireHistory*>::iterator it =
773  FullHistory.goodBranches.begin();
774  it != FullHistory.goodBranches.end(); ++it ) {
775 
776  // Double to access path.
777  double indexNow = (lastp + 0.5*(it->first - lastp))/sumAll;
778  lastp = it->first;
779 
780  // Probability of path.
781  double probPath = it->second->prodOfProbsFull/sumFullAll;
782 
783  FullHistory.select(indexNow)->setSelectedChild();
784 
785  // Calculate CKKWL weight:
786  double w = FullHistory.weightTREE( trialPartonLevelPtr,
787  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
788  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
789  indexNow);
790 
791  wgtsum += probPath*w;
792  }
793 
794  wgt = wgtsum;
795 
796  // Event with production scales set for further (trial) showering
797  // and starting conditions for the shower.
798  FullHistory.getStartingConditions( RN, process );
799  // If necessary, reattach resonance decay products.
800  mergingHooksPtr->reattachResonanceDecays(process);
801 
802  // Allow to dampen histories in which the lowest multiplicity reclustered
803  // state does not pass the lowest multiplicity cut of the matrix element.
804  double dampWeight = mergingHooksPtr->dampenIfFailCuts(
805  FullHistory.lowestMultProc(RN) );
806  // Save the weight of the event for histogramming. Only change the
807  // event weight after trial shower on the matrix element
808  // multiplicity event (= in doVetoStep).
809  wgt *= dampWeight;
810 
811  // Save the weight of the event for histogramming.
812  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({wgt});
813 
814  // Update the event weight.
815  if ( includeWGT) infoPtr->weightContainerPtr->
816  setWeightNominal(infoPtr->weight()*wgt);
817 
818  // Allow merging hooks to veto events from now on.
819  mergingHooksPtr->doIgnoreStep(false);
820 
821  // If no-emission probability is zero.
822  if ( allowReject && wgt == 0. ) return 0;
823 
824  // Done
825  return 1;
826 
827 }
828 
829 //--------------------------------------------------------------------------
830 
831 // Function to perform UMEPS merging on this event.
832 
833 int DireMerging::mergeProcessUMEPS( Event& process) {
834 
835  // Initialise which part of UMEPS merging is applied.
836  bool doUMEPSTree = settingsPtr->flag("Merging:doUMEPSTree");
837  bool doUMEPSSubt = settingsPtr->flag("Merging:doUMEPSSubt");
838  // Save number of looping steps
839  mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
840  int nRecluster = settingsPtr->mode("Merging:nRecluster");
841 
842  // Ensure that merging hooks does not remove emissions.
843  mergingHooksPtr->doIgnoreEmissions(true);
844  // For pp > h, allow cut on state, so that underlying processes
845  // can be clustered to gg > h
846  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0 )
847  mergingHooksPtr->allowCutOnRecState(true);
848  // For now, prefer construction of ordered histories.
849  mergingHooksPtr->orderHistories(true);
850 
851  // Ensure that merging weight is not counted twice.
852  bool includeWGT = mergingHooksPtr->includeWGTinXSEC();
853 
854  // Reset any incoming spins for W+-.
855  if (mergingHooksPtr->doWeakClustering())
856  for (int i = 0;i < process.size();++i)
857  process[i].pol(9);
858 
859  // Reset weights of the event.
860  double wgt = 1.;
861  mergingHooksPtr->setWeightCKKWL({1.});
862  mergingHooksPtr->muMI(-1.);
863 
864  // Prepare process record for merging. If Pythia has already decayed
865  // resonances used to define the hard process, remove resonance decay
866  // products.
867  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
868  // Store candidates for the splitting V -> qqbar'.
869  mergingHooksPtr->storeHardProcessCandidates( newProcess );
870 
871  // Check if event passes the merging scale cut.
872  double tmsval = mergingHooksPtr->tms();
873  // Get merging scale in current event.
874  double tmsnow = mergingHooksPtr->tmsNow( newProcess );
875  // Calculate number of clustering steps.
876  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
877  int nRequested = mergingHooksPtr->nRequested();
878 
879  // Too few steps can be possible if a chain of resonance decays has been
880  // removed. In this case, reject this event, since it will be handled in
881  // lower-multiplicity samples.
882  if (nSteps < nRequested) {
883  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({0.});
884  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
885  return -1;
886  }
887 
888  // Reset the minimal tms value, if necessary.
889  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
890 
891  // Get random number to choose a path.
892  double RN = rndmPtr->flat();
893  // Set dummy process scale.
894  newProcess.scale(0.0);
895  // Generate all histories.
896  DireHistory FullHistory( nSteps, 0.0, newProcess, DireClustering(),
897  mergingHooksPtr,
898  (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
899  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
900  1.0, 1.0, 1.0, 1.0, 0);
901  // Project histories onto desired branches, e.g. only ordered paths.
902  FullHistory.projectOntoDesiredHistories();
903 
904  // Do not apply cut if the configuration could not be projected onto an
905  // underlying born configuration.
906  bool applyCut = nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
907 
908  // Enfore merging scale cut if the event did not pass the merging scale
909  // criterion.
910  if ( enforceCutOnLHE && applyCut && tmsnow < tmsval ) {
911  string message="Warning in DireMerging::mergeProcessUMEPS: "
912  "Les Houches Event";
913  message+=" fails merging scale cut. Reject event.";
914  infoPtr->errorMsg(message);
915  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({0.});
916  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
917  return -1;
918  }
919 
920  // Check reclustering steps to correctly apply MPI.
921  int nPerformed = 0;
922  if ( nSteps > 0 && doUMEPSSubt
923  && !FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster,
924  newProcess, nPerformed, false ) ) {
925  // Discard if the state could not be reclustered to a state above TMS.
926  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({0.});
927  if ( includeWGT) infoPtr->weightContainerPtr->setWeightNominal(0.);
928  return -1;
929  }
930 
931  mergingHooksPtr->nMinMPI(nSteps - nPerformed);
932 
933  // Calculate CKKWL weight:
934  // Perform reweighting with Sudakov factors, save alpha_s ratios and
935  // PDF ratio weights.
936  if ( doUMEPSTree ) {
937  wgt = FullHistory.weight_UMEPS_TREE( trialPartonLevelPtr,
938  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
939  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(), RN);
940  } else {
941  wgt = FullHistory.weight_UMEPS_SUBT( trialPartonLevelPtr,
942  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
943  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(), RN);
944  }
945 
946  // Event with production scales set for further (trial) showering
947  // and starting conditions for the shower.
948  if ( doUMEPSTree ) FullHistory.getStartingConditions( RN, process );
949  // Do reclustering (looping) steps.
950  else FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster, process,
951  nPerformed, true );
952 
953  // Allow to dampen histories in which the lowest multiplicity reclustered
954  // state does not pass the lowest multiplicity cut of the matrix element
955  double dampWeight = mergingHooksPtr->dampenIfFailCuts(
956  FullHistory.lowestMultProc(RN) );
957  // Save the weight of the event for histogramming. Only change the
958  // event weight after trial shower on the matrix element
959  // multiplicity event (= in doVetoStep)
960  wgt *= dampWeight;
961 
962  // Save the weight of the event for histogramming.
963  if (!includeWGT) mergingHooksPtr->setWeightCKKWL({wgt});
964 
965  // Update the event weight.
966  if ( includeWGT) infoPtr->weightContainerPtr->
967  setWeightNominal(infoPtr->weight()*wgt);
968 
969  // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
970  // --> Set to minimal mT of partons.
971  int nFinal = 0;
972  double muf = process[0].e();
973  for ( int i=0; i < process.size(); ++i )
974  if ( process[i].isFinal()
975  && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
976  nFinal++;
977  muf = min( muf, abs(process[i].mT()) );
978  }
979 
980  // For pure QCD dijet events (only!), set the process scale to the
981  // transverse momentum of the outgoing partons.
982  // Calculate number of clustering steps.
983  int nStepsNew = mergingHooksPtr->getNumberOfClusteringSteps( process );
984  if ( nStepsNew == 0
985  && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
986  || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
987  process.scale(muf);
988 
989  // Reset hard process candidates (changed after clustering a parton).
990  mergingHooksPtr->storeHardProcessCandidates( process );
991  // If necessary, reattach resonance decay products.
992  mergingHooksPtr->reattachResonanceDecays(process);
993 
994  // Allow merging hooks to remove emissions from now on.
995  mergingHooksPtr->doIgnoreEmissions(false);
996 
997  // If no-emission probability is zero.
998  if ( wgt == 0. ) return 0;
999 
1000  // Done
1001  return 1;
1002 
1003 }
1004 
1005 //--------------------------------------------------------------------------
1006 
1007 // Function to perform NL3 NLO merging on this event.
1008 
1009 int DireMerging::mergeProcessNL3( Event& process) {
1010 
1011  // Initialise which part of NL3 merging is applied.
1012  bool doNL3Tree = settingsPtr->flag("Merging:doNL3Tree");
1013  bool doNL3Loop = settingsPtr->flag("Merging:doNL3Loop");
1014  bool doNL3Subt = settingsPtr->flag("Merging:doNL3Subt");
1015 
1016  // Ensure that hooks (NL3 part) to not remove emissions.
1017  mergingHooksPtr->doIgnoreEmissions(true);
1018  // Ensure that hooks (CKKWL part) to not veto events in trial showers.
1019  mergingHooksPtr->doIgnoreStep(true);
1020  // For pp > h, allow cut on state, so that underlying processes
1021  // can be clustered to gg > h
1022  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
1023  mergingHooksPtr->allowCutOnRecState(true);
1024  // For now, prefer construction of ordered histories.
1025  mergingHooksPtr->orderHistories(true);
1026 
1027  // Reset weight of the event
1028  double wgt = 1.;
1029  mergingHooksPtr->setWeightCKKWL({1.});
1030  // Reset the O(alphaS)-term of the CKKW-L weight.
1031  double wgtFIRST = 0.;
1032  mergingHooksPtr->setWeightFIRST({0.});
1033  mergingHooksPtr->muMI(-1.);
1034 
1035  // Prepare process record for merging. If Pythia has already decayed
1036  // resonances used to define the hard process, remove resonance decay
1037  // products.
1038  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
1039  // Store candidates for the splitting V -> qqbar'
1040  mergingHooksPtr->storeHardProcessCandidates( newProcess);
1041 
1042  // Check if event passes the merging scale cut.
1043  double tmsval = mergingHooksPtr->tms();
1044  // Get merging scale in current event.
1045  double tmsnow = mergingHooksPtr->tmsNow( newProcess );
1046  // Calculate number of clustering steps
1047  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
1048  int nRequested = mergingHooksPtr->nRequested();
1049 
1050  // Too few steps can be possible if a chain of resonance decays has been
1051  // removed. In this case, reject this event, since it will be handled in
1052  // lower-multiplicity samples.
1053  if (nSteps < nRequested) {
1054  mergingHooksPtr->setWeightCKKWL({0.});
1055  mergingHooksPtr->setWeightFIRST({0.});
1056  return -1;
1057  }
1058 
1059  // Reset the minimal tms value, if necessary.
1060  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
1061 
1062  // Enfore merging scale cut if the event did not pass the merging scale
1063  // criterion.
1064  if ( enforceCutOnLHE && nSteps > 0 && nSteps == nRequested
1065  && tmsnow < tmsval ) {
1066  string message = "Warning in DireMerging::mergeProcessNL3: Les Houches";
1067  message += " Event fails merging scale cut. Reject event.";
1068  infoPtr->errorMsg(message);
1069  mergingHooksPtr->setWeightCKKWL({0.});
1070  mergingHooksPtr->setWeightFIRST({0.});
1071  return -1;
1072  }
1073 
1074  // Get random number to choose a path.
1075  double RN = rndmPtr->flat();
1076  // Set dummy process scale.
1077  newProcess.scale(0.0);
1078  // Generate all histories
1079  DireHistory FullHistory( nSteps, 0.0, newProcess, DireClustering(),
1080  mergingHooksPtr,
1081  (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
1082  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
1083  1.0, 1.0, 1.0, 1.0, 0);
1084  // Project histories onto desired branches, e.g. only ordered paths.
1085  FullHistory.projectOntoDesiredHistories();
1086 
1087  // Discard states that cannot be projected unto a state with one less jet.
1088  if ( nSteps > 0 && doNL3Subt
1089  && FullHistory.select(RN)->nClusterings() == 0 ){
1090  mergingHooksPtr->setWeightCKKWL({0.});
1091  mergingHooksPtr->setWeightFIRST({0.});
1092  return -1;
1093  }
1094 
1095  // Potentially recluster real emission jets for powheg input containing
1096  // "too many" jets, i.e. real-emission kinematics.
1097  bool containsRealKin = nSteps > nRequested && nSteps > 0;
1098 
1099  // Perform one reclustering for real emission kinematics, then apply merging
1100  // scale cut on underlying Born kinematics.
1101  if ( containsRealKin ) {
1102  Event dummy = Event();
1103  // Initialise temporary output of reclustering.
1104  dummy.clear();
1105  dummy.init( "(hard process-modified)", particleDataPtr );
1106  dummy.clear();
1107  // Recluster once.
1108  if ( !FullHistory.getClusteredEvent( RN, nSteps, dummy )) {
1109  mergingHooksPtr->setWeightCKKWL({0.});
1110  mergingHooksPtr->setWeightFIRST({0.});
1111  return -1;
1112  }
1113  double tnowNew = mergingHooksPtr->tmsNow( dummy );
1114  // Veto if underlying Born kinematics do not pass merging scale cut.
1115  if ( enforceCutOnLHE && nRequested > 0 && tnowNew < tmsval ) {
1116  mergingHooksPtr->setWeightCKKWL({0.});
1117  mergingHooksPtr->setWeightFIRST({0.});
1118  return -1;
1119  }
1120  }
1121 
1122  // Remember number of jets, to include correct MPI no-emission probabilities.
1123  if ( doNL3Subt || containsRealKin ) mergingHooksPtr->nMinMPI(nSteps - 1);
1124  else mergingHooksPtr->nMinMPI(nSteps);
1125 
1126  // Calculate weight
1127  // Do LO or first part of NLO tree-level reweighting
1128  if( doNL3Tree ) {
1129  // Perform reweighting with Sudakov factors, save as ratios and
1130  // PDF ratio weights
1131  wgt = FullHistory.weightTREE( trialPartonLevelPtr,
1132  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1133  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(), RN);
1134  } else if( doNL3Loop || doNL3Subt ) {
1135  // No reweighting, just set event scales properly and incorporate MPI
1136  // no-emission probabilities.
1137  wgt = FullHistory.weightLOOP( trialPartonLevelPtr, RN);
1138  }
1139 
1140  // Event with production scales set for further (trial) showering
1141  // and starting conditions for the shower
1142  if ( !doNL3Subt && !containsRealKin )
1143  FullHistory.getStartingConditions(RN, process);
1144  // For sutraction of nSteps-additional resolved partons from
1145  // the nSteps-1 parton phase space, recluster the last parton
1146  // in nSteps-parton events, and sutract later
1147  else {
1148  // Function to return the reclustered event
1149  if ( !FullHistory.getClusteredEvent( RN, nSteps, process )) {
1150  mergingHooksPtr->setWeightCKKWL({0.});
1151  mergingHooksPtr->setWeightFIRST({0.});
1152  return -1;
1153  }
1154  }
1155 
1156  // Allow to dampen histories in which the lowest multiplicity reclustered
1157  // state does not pass the lowest multiplicity cut of the matrix element
1158  double dampWeight = mergingHooksPtr->dampenIfFailCuts(
1159  FullHistory.lowestMultProc(RN) );
1160  // Save the weight of the event for histogramming. Only change the
1161  // event weight after trial shower on the matrix element
1162  // multiplicity event (= in doVetoStep)
1163  wgt *= dampWeight;
1164 
1165  // For tree level samples in NL3, rescale with k-Factor
1166  if (doNL3Tree ){
1167  // Find k-factor
1168  double kFactor = 1.;
1169  if( nSteps > mergingHooksPtr->nMaxJetsNLO() )
1170  kFactor = mergingHooksPtr->kFactor( mergingHooksPtr->nMaxJetsNLO() );
1171  else kFactor = mergingHooksPtr->kFactor(nSteps);
1172  // For NLO merging, rescale CKKW-L weight with k-factor
1173  wgt *= kFactor;
1174  }
1175 
1176  // Save the weight of the event for histogramming
1177  mergingHooksPtr->setWeightCKKWL({wgt});
1178 
1179  // Check if we need to subtract the O(\alpha_s)-term. If the number
1180  // of additional partons is larger than the number of jets for
1181  // which loop matrix elements are available, do standard CKKW-L
1182  bool doOASTree = doNL3Tree && nSteps <= mergingHooksPtr->nMaxJetsNLO();
1183 
1184  // Now begin NLO part for tree-level events
1185  if ( doOASTree ) {
1186  // Calculate the O(\alpha_s)-term of the CKKWL weight
1187  wgtFIRST = FullHistory.weightFIRST( trialPartonLevelPtr,
1188  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1189  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(), RN,
1190  rndmPtr );
1191  // If necessary, also dampen the O(\alpha_s)-term
1192  wgtFIRST *= dampWeight;
1193  // Set the subtractive weight to the value calculated so far
1194  mergingHooksPtr->setWeightFIRST({wgtFIRST});
1195  // Subtract the O(\alpha_s)-term from the CKKW-L weight
1196  // If PDF contributions have not been included, subtract these later
1197  wgt = wgt - wgtFIRST;
1198  }
1199 
1200  // Set qcd 2->2 starting scale different from arbirtrary scale in LHEF!
1201  // --> Set to pT of partons
1202  double pT = 0.;
1203  for( int i=0; i < process.size(); ++i)
1204  if(process[i].isFinal() && process[i].colType() != 0) {
1205  pT = sqrt(pow(process[i].px(),2) + pow(process[i].py(),2));
1206  break;
1207  }
1208  // For pure QCD dijet events (only!), set the process scale to the
1209  // transverse momentum of the outgoing partons.
1210  if ( nSteps == 0
1211  && mergingHooksPtr->getProcessString().compare("pp>jj") == 0)
1212  process.scale(pT);
1213 
1214  // Reset hard process candidates (changed after clustering a parton).
1215  mergingHooksPtr->storeHardProcessCandidates( process );
1216  // If necessary, reattach resonance decay products.
1217  mergingHooksPtr->reattachResonanceDecays(process);
1218 
1219  // Allow merging hooks (NL3 part) to remove emissions from now on.
1220  mergingHooksPtr->doIgnoreEmissions(false);
1221  // Allow merging hooks (CKKWL part) to veto events from now on.
1222  mergingHooksPtr->doIgnoreStep(false);
1223 
1224  // Done
1225  return 1;
1226 
1227 }
1228 
1229 //--------------------------------------------------------------------------
1230 
1231 // Function to perform UNLOPS merging on this event.
1232 
1233 int DireMerging::mergeProcessUNLOPS( Event& process) {
1234 
1235  // Initialise which part of UNLOPS merging is applied.
1236  bool nloTilde = settingsPtr->flag("Merging:doUNLOPSTilde");
1237  bool doUNLOPSTree = settingsPtr->flag("Merging:doUNLOPSTree");
1238  bool doUNLOPSLoop = settingsPtr->flag("Merging:doUNLOPSLoop");
1239  bool doUNLOPSSubt = settingsPtr->flag("Merging:doUNLOPSSubt");
1240  bool doUNLOPSSubtNLO = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
1241  // Save number of looping steps
1242  mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
1243  int nRecluster = settingsPtr->mode("Merging:nRecluster");
1244 
1245  // Ensure that merging hooks to not remove emissions
1246  mergingHooksPtr->doIgnoreEmissions(true);
1247  // For now, prefer construction of ordered histories.
1248  mergingHooksPtr->orderHistories(true);
1249  // For pp > h, allow cut on state, so that underlying processes
1250  // can be clustered to gg > h
1251  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
1252  mergingHooksPtr->allowCutOnRecState(true);
1253 
1254  // Reset weight of the event.
1255  double wgt = 1.;
1256  mergingHooksPtr->setWeightCKKWL({1.});
1257  // Reset the O(alphaS)-term of the UMEPS weight.
1258  double wgtFIRST = 0.;
1259  mergingHooksPtr->setWeightFIRST({0.});
1260  mergingHooksPtr->muMI(-1.);
1261 
1262  // Prepare process record for merging. If Pythia has already decayed
1263  // resonances used to define the hard process, remove resonance decay
1264  // products.
1265  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
1266  // Store candidates for the splitting V -> qqbar'
1267  mergingHooksPtr->storeHardProcessCandidates( newProcess );
1268 
1269  // Check if event passes the merging scale cut.
1270  double tmsval = mergingHooksPtr->tms();
1271  // Get merging scale in current event.
1272  double tmsnow = mergingHooksPtr->tmsNow( newProcess );
1273  // Calculate number of clustering steps
1274  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
1275  int nRequested = mergingHooksPtr->nRequested();
1276 
1277  // Too few steps can be possible if a chain of resonance decays has been
1278  // removed. In this case, reject this event, since it will be handled in
1279  // lower-multiplicity samples.
1280  if (nSteps < nRequested) {
1281  string message="Warning in DireMerging::mergeProcessUNLOPS: "
1282  "Les Houches Event";
1283  message+=" after removing decay products does not contain enough partons.";
1284  infoPtr->errorMsg(message);
1285  mergingHooksPtr->setWeightCKKWL({0.});
1286  mergingHooksPtr->setWeightFIRST({0.});
1287  return -1;
1288  }
1289 
1290  // Reset the minimal tms value, if necessary.
1291  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
1292 
1293  // Get random number to choose a path.
1294  double RN = rndmPtr->flat();
1295  // Set dummy process scale.
1296  newProcess.scale(0.0);
1297  // Generate all histories
1298  DireHistory FullHistory( nSteps, 0.0, newProcess, DireClustering(),
1299  mergingHooksPtr,
1300  (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
1301  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
1302  1.0, 1.0, 1.0, 1.0, 0);
1303  // Project histories onto desired branches, e.g. only ordered paths.
1304  FullHistory.projectOntoDesiredHistories();
1305 
1306  // Do not apply cut if the configuration could not be projected onto an
1307  // underlying born configuration.
1308  bool applyCut = nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
1309 
1310  // Enfore merging scale cut if the event did not pass the merging scale
1311  // criterion.
1312  if ( enforceCutOnLHE && applyCut && nSteps == nRequested
1313  && tmsnow < tmsval && tmsval > 0.) {
1314  string message="Warning in DireMerging::mergeProcessUNLOPS: Les Houches";
1315  message+=" Event fails merging scale cut. Reject event.";
1316  infoPtr->errorMsg(message);
1317  mergingHooksPtr->setWeightCKKWL({0.});
1318  mergingHooksPtr->setWeightFIRST({0.});
1319  return -1;
1320  }
1321 
1322  // Potentially recluster real emission jets for powheg input containing
1323  // "too many" jets, i.e. real-emission kinematics.
1324  bool containsRealKin = nSteps > nRequested && nSteps > 0;
1325  if ( containsRealKin ) nRecluster += nSteps - nRequested;
1326 
1327  // Remove real emission events without underlying Born configuration from
1328  // the loop sample, since such states will be taken care of by tree-level
1329  // samples.
1330  if ( doUNLOPSLoop && containsRealKin && !allowIncompleteReal
1331  && FullHistory.select(RN)->nClusterings() == 0 ) {
1332  mergingHooksPtr->setWeightCKKWL({0.});
1333  mergingHooksPtr->setWeightFIRST({0.});
1334  return -1;
1335  }
1336 
1337  // Discard if the state could not be reclustered to any state above TMS.
1338  int nPerformed = 0;
1339  if ( nSteps > 0 && !allowIncompleteReal
1340  && ( doUNLOPSSubt || doUNLOPSSubtNLO || containsRealKin )
1341  && !FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster,
1342  newProcess, nPerformed, false ) ) {
1343  mergingHooksPtr->setWeightCKKWL({0.});
1344  mergingHooksPtr->setWeightFIRST({0.});
1345  return -1;
1346  }
1347 
1348  // Check reclustering steps to correctly apply MPI.
1349  mergingHooksPtr->nMinMPI(nSteps - nPerformed);
1350 
1351  // Perform one reclustering for real emission kinematics, then apply
1352  // merging scale cut on underlying Born kinematics.
1353  if ( containsRealKin ) {
1354  Event dummy = Event();
1355  // Initialise temporary output of reclustering.
1356  dummy.clear();
1357  dummy.init( "(hard process-modified)", particleDataPtr );
1358  dummy.clear();
1359  // Recluster once.
1360  FullHistory.getClusteredEvent( RN, nSteps, dummy );
1361  double tnowNew = mergingHooksPtr->tmsNow( dummy );
1362  // Veto if underlying Born kinematics do not pass merging scale cut.
1363  if (enforceCutOnLHE && nRequested > 0 && tnowNew < tmsval && tmsval > 0.) {
1364  string message="Warning in DireMerging::mergeProcessUNLOPS: Les Houches";
1365  message+=" Event fails merging scale cut. Reject event.";
1366  infoPtr->errorMsg(message);
1367  mergingHooksPtr->setWeightCKKWL({0.});
1368  mergingHooksPtr->setWeightFIRST({0.});
1369  return -1;
1370  }
1371  }
1372 
1373  // New UNLOPS strategy based on UN2LOPS.
1374  bool doUNLOPS2 = false;
1375  int depth = -1;
1376 
1377  // Calculate weights.
1378  // Do LO or first part of NLO tree-level reweighting
1379  if( doUNLOPSTree ) {
1380  // Perform reweighting with Sudakov factors, save as ratios and
1381  // PDF ratio weights
1382  wgt = FullHistory.weight_UNLOPS_TREE( trialPartonLevelPtr,
1383  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1384  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1385  RN, depth);
1386  } else if( doUNLOPSLoop ) {
1387  // Set event scales properly, reweight for new UNLOPS
1388  wgt = FullHistory.weight_UNLOPS_LOOP( trialPartonLevelPtr,
1389  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1390  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1391  RN, depth);
1392  } else if( doUNLOPSSubtNLO ) {
1393  // Set event scales properly, reweight for new UNLOPS
1394  wgt = FullHistory.weight_UNLOPS_SUBTNLO( trialPartonLevelPtr,
1395  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1396  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1397  RN, depth);
1398  } else if( doUNLOPSSubt ) {
1399  // Perform reweighting with Sudakov factors, save as ratios and
1400  // PDF ratio weights
1401  wgt = FullHistory.weight_UNLOPS_SUBT( trialPartonLevelPtr,
1402  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1403  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1404  RN, depth);
1405  }
1406 
1407  // Event with production scales set for further (trial) showering
1408  // and starting conditions for the shower.
1409  if (!doUNLOPSSubt && !doUNLOPSSubtNLO && !containsRealKin )
1410  FullHistory.getStartingConditions(RN, process);
1411  // Do reclustering (looping) steps.
1412  else FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster, process,
1413  nPerformed, true );
1414 
1415  // Allow to dampen histories in which the lowest multiplicity reclustered
1416  // state does not pass the lowest multiplicity cut of the matrix element
1417  double dampWeight = mergingHooksPtr->dampenIfFailCuts(
1418  FullHistory.lowestMultProc(RN) );
1419  // Save the weight of the event for histogramming. Only change the
1420  // event weight after trial shower on the matrix element
1421  // multiplicity event (= in doVetoStep)
1422  wgt *= dampWeight;
1423 
1424  // For tree-level or subtractive sammples, rescale with k-Factor
1425  if ( doUNLOPSTree || doUNLOPSSubt ){
1426  // Find k-factor
1427  double kFactor = 1.;
1428  if ( nSteps > mergingHooksPtr->nMaxJetsNLO() )
1429  kFactor = mergingHooksPtr->kFactor( mergingHooksPtr->nMaxJetsNLO() );
1430  else kFactor = mergingHooksPtr->kFactor(nSteps);
1431  // For NLO merging, rescale CKKW-L weight with k-factor
1432  wgt *= (nRecluster == 2 && nloTilde) ? 1. : kFactor;
1433  }
1434 
1435  // Save the weight of the event for histogramming
1436  mergingHooksPtr->setWeightCKKWL({wgt});
1437 
1438  // Check if we need to subtract the O(\alpha_s)-term. If the number
1439  // of additional partons is larger than the number of jets for
1440  // which loop matrix elements are available, do standard UMEPS.
1441  int nMaxNLO = mergingHooksPtr->nMaxJetsNLO();
1442  bool doOASTree = doUNLOPSTree && nSteps <= nMaxNLO;
1443  bool doOASSubt = doUNLOPSSubt && nSteps <= nMaxNLO+1 && nSteps > 0;
1444 
1445  // Now begin NLO part for tree-level events
1446  if ( doOASTree || doOASSubt ) {
1447 
1448  // Decide on which order to expand to.
1449  int order = ( nSteps > 0 && nSteps <= nMaxNLO) ? 1 : -1;
1450 
1451  // Exclusive inputs:
1452  // Subtract only the O(\alpha_s^{n+0})-term from the tree-level
1453  // subtraction, if we're at the highest NLO multiplicity (nMaxNLO).
1454  if ( nloTilde && doUNLOPSSubt && nRecluster == 1
1455  && nSteps == nMaxNLO+1 ) order = 0;
1456 
1457  // Exclusive inputs:
1458  // Do not remove the O(as)-term if the number of reclusterings
1459  // exceeds the number of NLO jets, or if more clusterings have
1460  // been performed.
1461  if (nloTilde && doUNLOPSSubt && ( nSteps > nMaxNLO+1
1462  || (nSteps == nMaxNLO+1 && nPerformed != nRecluster) ))
1463  order = -1;
1464 
1465  // Calculate terms in expansion of the CKKW-L weight.
1466  wgtFIRST = FullHistory.weight_UNLOPS_CORRECTION( order,
1467  trialPartonLevelPtr, mergingHooksPtr->AlphaS_FSR(),
1468  mergingHooksPtr->AlphaS_ISR(), mergingHooksPtr->AlphaEM_FSR(),
1469  mergingHooksPtr->AlphaEM_ISR(), RN, rndmPtr );
1470 
1471  // Exclusive inputs:
1472  // Subtract the O(\alpha_s^{n+1})-term from the tree-level
1473  // subtraction, not the O(\alpha_s^{n+0})-terms.
1474  if ( nloTilde && doUNLOPSSubt && nRecluster == 1
1475  && nPerformed == nRecluster && nSteps <= nMaxNLO )
1476  wgtFIRST += 1.;
1477 
1478  // If necessary, also dampen the O(\alpha_s)-term
1479  wgtFIRST *= dampWeight;
1480 
1481  // Set the subtractive weight to the value calculated so far
1482  mergingHooksPtr->setWeightFIRST({wgtFIRST});
1483  // Subtract the O(\alpha_s)-term from the CKKW-L weight
1484  // If PDF contributions have not been included, subtract these later
1485  // New UNLOPS based on UN2LOPS.
1486  if (doUNLOPS2 && order > -1) wgt = -wgt*(wgtFIRST-1.);
1487  else if (order > -1) wgt = wgt - wgtFIRST;
1488 
1489  }
1490 
1491  // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
1492  // --> Set to minimal mT of partons.
1493  int nFinal = 0;
1494  double muf = process[0].e();
1495  for ( int i=0; i < process.size(); ++i )
1496  if ( process[i].isFinal()
1497  && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
1498  nFinal++;
1499  muf = min( muf, abs(process[i].mT()) );
1500  }
1501  // For pure QCD dijet events (only!), set the process scale to the
1502  // transverse momentum of the outgoing partons.
1503  if ( nSteps == 0 && nFinal == 2
1504  && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
1505  || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
1506  process.scale(muf);
1507 
1508  // Reset hard process candidates (changed after clustering a parton).
1509  mergingHooksPtr->storeHardProcessCandidates( process );
1510 
1511  // Check if resonance structure has been changed
1512  // (e.g. because of clustering W/Z/gluino)
1513  vector <int> oldResonance;
1514  for ( int i=0; i < newProcess.size(); ++i )
1515  if ( newProcess[i].status() == 22 )
1516  oldResonance.push_back(newProcess[i].id());
1517  vector <int> newResonance;
1518  for ( int i=0; i < process.size(); ++i )
1519  if ( process[i].status() == 22 )
1520  newResonance.push_back(process[i].id());
1521  // Compare old and new resonances
1522  for ( int i=0; i < int(oldResonance.size()); ++i )
1523  for ( int j=0; j < int(newResonance.size()); ++j )
1524  if ( newResonance[j] == oldResonance[i] ) {
1525  oldResonance[i] = 99;
1526  break;
1527  }
1528  bool hasNewResonances = (newResonance.size() != oldResonance.size());
1529  for ( int i=0; i < int(oldResonance.size()); ++i )
1530  hasNewResonances = (oldResonance[i] != 99);
1531 
1532  // If necessary, reattach resonance decay products.
1533  if (!hasNewResonances) mergingHooksPtr->reattachResonanceDecays(process);
1534 
1535  // Allow merging hooks to remove emissions from now on.
1536  mergingHooksPtr->doIgnoreEmissions(false);
1537 
1538  // If no-emission probability is zero.
1539  if ( wgt == 0. ) return 0;
1540 
1541  // If the resonance structure of the process has changed due to reclustering,
1542  // redo the resonance decays in Pythia::next()
1543  if (hasNewResonances) return 2;
1544 
1545  // Done
1546  return 1;
1547 
1548 }
1549 
1550 //--------------------------------------------------------------------------
1551 
1552 // Function to set up all histories for an event.
1553 
1554 bool DireMerging::generateHistories( const Event& process, bool orderedOnly) {
1555 
1556  // Input not valid.
1557  if (!validEvent(process)) {
1558  cout << "Warning in DireMerging::generateHistories: Input event "
1559  << "has invalid flavour or momentum structure, thus reject. " << endl;
1560  return false;
1561  }
1562 
1563  // Clear previous history.
1564  if (myHistory) delete myHistory;
1565 
1566  // For now, prefer construction of ordered histories.
1567  mergingHooksPtr->orderHistories(orderedOnly);
1568 
1569  if (doMOPS) mergingHooksPtr->orderHistories(false);
1570 
1571  // For pp > h, allow cut on state, so that underlying processes
1572  // can be clustered to gg > h
1573  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
1574  mergingHooksPtr->allowCutOnRecState(true);
1575 
1576  // Prepare process record for merging. If Pythia has already decayed
1577  // resonances used to define the hard process, remove resonance decay
1578  // products.
1579  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
1580  // Store candidates for the splitting V -> qqbar'
1581  mergingHooksPtr->storeHardProcessCandidates( newProcess );
1582 
1583  // Calculate number of clustering steps
1584  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
1585 
1586  nSteps++;
1587 
1588  // Set dummy process scale.
1589  newProcess.scale(0.0);
1590  // Generate all histories
1591  myHistory = new DireHistory( nSteps, 0.0, newProcess, DireClustering(),
1592  mergingHooksPtr,
1593  (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
1594  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
1595  1.0, 1.0, 1.0, 1.0, 0);
1596  // Project histories onto desired branches, e.g. only ordered paths.
1597  bool foundHistories = myHistory->projectOntoDesiredHistories();
1598 
1599  // Done
1600  return (doMOPS ? foundHistories : true);
1601 
1602 }
1603 
1604 //--------------------------------------------------------------------------
1605 
1606 void DireMerging::tagHistories() {
1607 
1608  // Tag history paths as "signal" or "background"
1609  for ( map<double, DireHistory*>::iterator it =
1610  myHistory->goodBranches.begin();
1611  it != myHistory->goodBranches.end(); ++it )
1612  it->second->tagPath(it->second);
1613 
1614  double sumAll(0.), sumFullAll(0.), lastp(0.);
1615  for ( map<double, DireHistory*>::iterator it =
1616  myHistory->goodBranches.begin();
1617  it != myHistory->goodBranches.end(); ++it ) {
1618  sumAll += it->second->prodOfProbs;
1619  sumFullAll += it->second->prodOfProbsFull;
1620  }
1621 
1622  // Sum up signal and background probabilities.
1623  vector<double> sumSignalProb(createvector<double>(0.)(0.)(0.)),
1624  sumBkgrndProb(createvector<double>(0.)(0.)(0.));
1625 
1626  for ( map<double, DireHistory*>::iterator it =
1627  myHistory->goodBranches.begin();
1628  it != myHistory->goodBranches.end(); ++it ) {
1629 
1630  if (it->second == myHistory) continue;
1631 
1632  // Get ME weight.
1633  double prob = it->second->prodOfProbsFull;
1634  // Reweight with Sudakovs, couplings and PDFs.
1635  double indexNow = (lastp + 0.5*(it->first - lastp))/sumAll;
1636  lastp = it->first;
1637  myHistory->select(indexNow)->setSelectedChild();
1638  vector<double> w = myHistory->weightMEM( trialPartonLevelPtr,
1639  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaEM_FSR(), indexNow);
1640  for (unsigned int i=0; i < w.size(); ++i) {
1641  totalProbSave[i] += prob*w[i];
1642  if (it->second->hasTag("higgs")) signalProbSave["higgs"][i] += prob*w[i];
1643  else bkgrndProbSave["higgs"][i] += prob*w[i];
1644  if (it->second->hasTag("qed")) signalProbSave["qed"][i] += prob*w[i];
1645  else bkgrndProbSave["qed"][i] += prob*w[i];
1646  if (it->second->hasTag("qcd")) signalProbSave["qcd"][i] += prob*w[i];
1647  else bkgrndProbSave["qcd"][i] += prob*w[i];
1648 
1649  if (it->second->hasTag("higgs") ) signalProbSave["higgs-subt"][i]
1650  += prob*(w[i]-1.);
1651  else bkgrndProbSave["higgs-subt"][i]
1652  += prob*(w[i]-1.);
1653  if (it->second->hasTag("higgs") ) signalProbSave["higgs-nosud"][i]
1654  += prob;
1655  else bkgrndProbSave["higgs-nosud"][i]
1656  += prob;
1657  }
1658  }
1659 
1660 }
1661 
1662 //--------------------------------------------------------------------------
1663 
1664 double DireMerging::getPathIndex( bool useAll) {
1665 
1666  if (!useAll) return rndmPtr->flat();
1667 
1668  // Setup to choose shower starting conditions randomly.
1669  double sumAll(0.);
1670  for ( map<double, DireHistory*>::iterator it =
1671  myHistory->goodBranches.begin();
1672  it != myHistory->goodBranches.end(); ++it ) {
1673  sumAll += it->second->prodOfProbs;
1674  }
1675  // Store a double with which to access each of the paths.
1676  double lastp(0.);
1677  vector<double> path_index;
1678  for ( map<double, DireHistory*>::iterator it =
1679  myHistory->goodBranches.begin();
1680  it != myHistory->goodBranches.end(); ++it ) {
1681  // Double to access path.
1682  double indexNow = (lastp + 0.5*(it->first - lastp))/sumAll;
1683  path_index.push_back(indexNow);
1684  lastp = it->first;
1685  }
1686  // Randomly pick path.
1687  int sizeBranches = myHistory->goodBranches.size();
1688  int iPosRN = (sizeBranches > 0)
1689  ? rndmPtr->pick(
1690  vector<double>(sizeBranches, 1./double(sizeBranches)) )
1691  : 0;
1692  double RN = (sizeBranches > 0) ? path_index[iPosRN] : rndmPtr->flat();
1693  return RN;
1694 }
1695 
1696 //--------------------------------------------------------------------------
1697 
1698 // Function to set up all histories for an event.
1699 
1700 bool DireMerging::calculateSubtractions() {
1701 
1702  // Store shower subtractions.
1703  clearSubtractions();
1704  for ( int i = 0 ; i < int(myHistory->children.size()); ++i) {
1705 
1706  // Need to reattach resonance decays, if necessary.
1707  Event psppoint = myHistory->children[i]->state;
1708  // Reset hard process candidates (changed after clustering a parton).
1709  mergingHooksPtr->storeHardProcessCandidates( psppoint );
1710 
1711  // Check if resonance structure has been changed
1712  // (e.g. because of clustering W/Z/gluino)
1713  vector <int> oldResonance;
1714  for ( int n=0; n < myHistory->state.size(); ++n )
1715  if ( myHistory->state[n].status() == 22 )
1716  oldResonance.push_back(myHistory->state[n].id());
1717  vector <int> newResonance;
1718  for ( int n=0; n < psppoint.size(); ++n )
1719  if ( psppoint[n].status() == 22 )
1720  newResonance.push_back(psppoint[n].id());
1721  // Compare old and new resonances
1722  for ( int n=0; n < int(oldResonance.size()); ++n )
1723  for ( int m=0; m < int(newResonance.size()); ++m )
1724  if ( newResonance[m] == oldResonance[n] ) {
1725  oldResonance[n] = 99;
1726  break;
1727  }
1728  bool hasNewResonances = (newResonance.size() != oldResonance.size());
1729  for ( int n=0; n < int(oldResonance.size()); ++n )
1730  hasNewResonances = (oldResonance[n] != 99);
1731 
1732  // If necessary, reattach resonance decay products.
1733  if (!hasNewResonances) mergingHooksPtr->reattachResonanceDecays(psppoint);
1734  else {
1735  cout << "Warning in DireMerging::generateHistories: Resonance "
1736  << "structure changed due to clustering. Cannot attach decay "
1737  << "products correctly." << endl;
1738  }
1739 
1740  double prob = myHistory->children[i]->clusterProb;
1741 
1742  // Switch from 4pi to 8pi convention
1743  prob *= 2.;
1744 
1745  // Get clustering variables.
1746  map<string,double> stateVars;
1747  int rad = myHistory->children[i]->clusterIn.radPos();
1748  int emt = myHistory->children[i]->clusterIn.emtPos();
1749  int rec = myHistory->children[i]->clusterIn.recPos();
1750 
1751  bool isFSR = myHistory->showers->timesPtr->isTimelike(myHistory->state,
1752  rad, emt, rec, "");
1753  if (isFSR)
1754  stateVars = myHistory->showers->timesPtr->getStateVariables(
1755  myHistory->state,rad,emt,rec,"");
1756  else
1757  stateVars = myHistory->showers->spacePtr->getStateVariables(
1758  myHistory->state,rad,emt,rec,"");
1759 
1760  double z = stateVars["z"];
1761  double t = stateVars["t"];
1762 
1763  double m2dip = abs
1764  (-2.*myHistory->state[emt].p()*myHistory->state[rad].p()
1765  -2.*myHistory->state[emt].p()*myHistory->state[rec].p()
1766  +2.*myHistory->state[rad].p()*myHistory->state[rec].p());
1767  double kappa2 = t/m2dip;
1768  double xCS = (z*(1-z) - kappa2) / (1 -z);
1769 
1770  // For II dipoles, scale with 1/xCS.
1771  prob *= 1./xCS;
1772 
1773  // Multiply with ME correction.
1774  prob *= myHistory->MECnum/myHistory->MECden;
1775 
1776  // Attach point to list of shower subtractions.
1777  appendSubtraction( prob, psppoint);
1778 
1779  }
1780 
1781  // Restore stored hard process candidates
1782  mergingHooksPtr->storeHardProcessCandidates( myHistory->state );
1783 
1784  // Done
1785  return true;
1786 
1787 }
1788 
1789 //--------------------------------------------------------------------------
1790 
1791 // Function to calulate the weights used for UNLOPS merging.
1792 
1793 int DireMerging::calculateWeights( double RNpath, bool useAll ) {
1794 
1795  // Initialise which part of UNLOPS merging is applied.
1796  bool nloTilde = settingsPtr->flag("Merging:doUNLOPSTilde");
1797  bool doUNLOPSTree = settingsPtr->flag("Merging:doUNLOPSTree");
1798  bool doUNLOPSLoop = settingsPtr->flag("Merging:doUNLOPSLoop");
1799  bool doUNLOPSSubt = settingsPtr->flag("Merging:doUNLOPSSubt");
1800  bool doUNLOPSSubtNLO = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
1801  // Save number of looping steps
1802  mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
1803  int nRecluster = settingsPtr->mode("Merging:nRecluster");
1804 
1805  // Ensure that merging hooks to not remove emissions
1806  mergingHooksPtr->doIgnoreEmissions(true);
1807  mergingHooksPtr->setWeightCKKWL({1.});
1808  mergingHooksPtr->setWeightFIRST({0.});
1809 
1810  // Reset weight of the event.
1811  double wgt = 1.;
1812  // Reset the O(alphaS)-term of the UMEPS weight.
1813  double wgtFIRST = 0.;
1814  mergingHooksPtr->muMI(-1.);
1815 
1816  // Check if event passes the merging scale cut.
1817  double tmsval = mergingHooksPtr->tms();
1818 
1819  if (doMOPS) tmsval = 0.;
1820 
1821  // Get merging scale in current event.
1822  double tmsnow = mergingHooksPtr->tmsNow( myHistory->state );
1823  // Calculate number of clustering steps
1824  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( myHistory->state,
1825  true);
1826  int nRequested = mergingHooksPtr->nRequested();
1827 
1828  if (doMOPS && nSteps == 0) { return 1; }
1829 
1830  // Too few steps can be possible if a chain of resonance decays has been
1831  // removed. In this case, reject this event, since it will be handled in
1832  // lower-multiplicity samples.
1833  if (nSteps < nRequested) {
1834  string message="Warning in DireMerging::calculateWeights: "
1835  "Les Houches Event";
1836  message+=" after removing decay products does not contain enough partons.";
1837  infoPtr->errorMsg(message);
1838  if (allowReject) return -1;
1839  }
1840 
1841  // Reset the minimal tms value, if necessary.
1842  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
1843 
1844  // Do not apply cut if the configuration could not be projected onto an
1845  // underlying born configuration.
1846  bool applyCut = nSteps > 0 && myHistory->select(RNpath)->nClusterings() > 0;
1847 
1848  // Enfore merging scale cut if the event did not pass the merging scale
1849  // criterion.
1850  if ( enforceCutOnLHE && applyCut && nSteps == nRequested
1851  && tmsnow < tmsval && tmsval > 0.) {
1852  string message="Warning in DireMerging::calculateWeights: Les Houches";
1853  message+=" Event fails merging scale cut. Reject event.";
1854  infoPtr->errorMsg(message);
1855  if (allowReject) return -1;
1856  //return -1;
1857  }
1858 
1859  // Potentially recluster real emission jets for powheg input containing
1860  // "too many" jets, i.e. real-emission kinematics.
1861  bool containsRealKin = nSteps > nRequested && nSteps > 0;
1862  if ( containsRealKin ) nRecluster += nSteps - nRequested;
1863 
1864  // Remove real emission events without underlying Born configuration from
1865  // the loop sample, since such states will be taken care of by tree-level
1866  // samples.
1867  if ( doUNLOPSLoop && containsRealKin && !allowIncompleteReal
1868  && myHistory->select(RNpath)->nClusterings() == 0 ) {
1869  if (allowReject) return -1;
1870  }
1871 
1872  // Discard if the state could not be reclustered to any state above TMS.
1873  int nPerformed = 0;
1874  if ( nSteps > 0 && !allowIncompleteReal
1875  && ( doUNLOPSSubt || doUNLOPSSubtNLO || containsRealKin )
1876  && !myHistory->getFirstClusteredEventAboveTMS( RNpath, nRecluster,
1877  myHistory->state, nPerformed, false ) ) {
1878  if (allowReject) return -1;
1879  }
1880 
1881  // Check reclustering steps to correctly apply MPI.
1882  mergingHooksPtr->nMinMPI(nSteps - nPerformed);
1883 
1884  // Perform one reclustering for real emission kinematics, then apply
1885  // merging scale cut on underlying Born kinematics.
1886  if ( containsRealKin ) {
1887  Event dummy = Event();
1888  // Initialise temporary output of reclustering.
1889  dummy.clear();
1890  dummy.init( "(hard process-modified)", particleDataPtr );
1891  dummy.clear();
1892  // Recluster once.
1893  myHistory->getClusteredEvent( RNpath, nSteps, dummy );
1894  double tnowNew = mergingHooksPtr->tmsNow( dummy );
1895  // Veto if underlying Born kinematics do not pass merging scale cut.
1896  if (enforceCutOnLHE && nRequested > 0 && tnowNew < tmsval && tmsval > 0.) {
1897  string message="Warning in DireMerging::calculateWeights: Les Houches";
1898  message+=" Event fails merging scale cut. Reject event.";
1899  infoPtr->errorMsg(message);
1900  if (allowReject) return -1;
1901  //return -1;
1902  }
1903  }
1904 
1905  // Setup to choose shower starting conditions randomly.
1906  double sumAll(0.), sumFullAll(0.);
1907  for ( map<double, DireHistory*>::iterator it =
1908  myHistory->goodBranches.begin();
1909  it != myHistory->goodBranches.end(); ++it ) {
1910  sumAll += it->second->prodOfProbs;
1911  sumFullAll += it->second->prodOfProbsFull;
1912  }
1913 
1914  // Sum up signal and background probabilities.
1915  double sumSignalProb(0.), sumBkgrndProb(0.);
1916 
1917  // New UNLOPS strategy based on UN2LOPS.
1918  bool doUNLOPS2 = false;
1919  int depth = -1;
1920 
1921  if (!useAll) {
1922 
1923  // Calculate weights.
1924  if (doMOPS)
1925  wgt = myHistory->weightMOPS( trialPartonLevelPtr,
1926  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaEM_FSR(),
1927  RNpath);
1928  else if ( mergingHooksPtr->doCKKWLMerging() )
1929  wgt = myHistory->weightTREE( trialPartonLevelPtr,
1930  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1931  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1932  RNpath);
1933  else if ( mergingHooksPtr->doUMEPSTreeSave )
1934  wgt = myHistory->weight_UMEPS_TREE( trialPartonLevelPtr,
1935  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1936  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1937  RNpath);
1938  else if ( mergingHooksPtr->doUMEPSSubtSave )
1939  wgt = myHistory->weight_UMEPS_SUBT( trialPartonLevelPtr,
1940  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1941  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1942  RNpath);
1943  else if ( mergingHooksPtr->doUNLOPSTreeSave )
1944  wgt = myHistory->weight_UNLOPS_TREE( trialPartonLevelPtr,
1945  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1946  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1947  RNpath, depth);
1948  else if ( mergingHooksPtr->doUNLOPSLoopSave )
1949  wgt = myHistory->weight_UNLOPS_LOOP( trialPartonLevelPtr,
1950  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1951  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1952  RNpath, depth);
1953  else if ( mergingHooksPtr->doUNLOPSSubtNLOSave )
1954  wgt = myHistory->weight_UNLOPS_SUBTNLO( trialPartonLevelPtr,
1955  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1956  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1957  RNpath, depth);
1958  else if ( mergingHooksPtr->doUNLOPSSubtSave )
1959  wgt = myHistory->weight_UNLOPS_SUBT( trialPartonLevelPtr,
1960  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(),
1961  mergingHooksPtr->AlphaEM_FSR(), mergingHooksPtr->AlphaEM_ISR(),
1962  RNpath, depth);
1963 
1964  // For tree-level or subtractive sammples, rescale with k-Factor
1965  if ( doUNLOPSTree || doUNLOPSSubt ){
1966  // Find k-factor
1967  double kFactor = 1.;
1968  if ( nSteps > mergingHooksPtr->nMaxJetsNLO() )
1969  kFactor = mergingHooksPtr->kFactor( mergingHooksPtr->nMaxJetsNLO() );
1970  else kFactor = mergingHooksPtr->kFactor(nSteps);
1971  // For NLO merging, rescale CKKW-L weight with k-factor
1972  wgt *= (nRecluster == 2 && nloTilde) ? 1. : kFactor;
1973  }
1974 
1975  } else if (doMOPS) {
1976  // Calculate CKKWL reweighting for all paths.
1977  double wgtsum(0.);
1978  double lastp(0.);
1979 
1980  for ( map<double, DireHistory*>::iterator it =
1981  myHistory->goodBranches.begin();
1982  it != myHistory->goodBranches.end(); ++it ) {
1983 
1984  // Double to access path.
1985  double indexNow = (lastp + 0.5*(it->first - lastp))/sumAll;
1986  lastp = it->first;
1987 
1988  // Probability of path.
1989  double probPath = it->second->prodOfProbsFull/sumFullAll;
1990 
1991  myHistory->select(indexNow)->setSelectedChild();
1992 
1993  // Calculate CKKWL weight:
1994  double w = myHistory->weightMOPS( trialPartonLevelPtr,
1995  mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaEM_FSR(),
1996  indexNow);
1997 
1998  wgtsum += probPath*w;
1999 
2000  // Sum up signal and background probabilities separately.
2001  if (it->second->hasTag("higgs"))
2002  sumSignalProb += it->second->prodOfProbsFull*w;
2003  else
2004  sumBkgrndProb += it->second->prodOfProbsFull*w;
2005 
2006  }
2007  wgt = wgtsum;
2008 
2009  }
2010 
2011  mergingHooksPtr->setWeightCKKWL({wgt});
2012 
2013  // Check if we need to subtract the O(\alpha_s)-term. If the number
2014  // of additional partons is larger than the number of jets for
2015  // which loop matrix elements are available, do standard UMEPS.
2016  int nMaxNLO = mergingHooksPtr->nMaxJetsNLO();
2017  bool doOASTree = doUNLOPSTree && nSteps <= nMaxNLO;
2018  bool doOASSubt = doUNLOPSSubt && nSteps <= nMaxNLO+1 && nSteps > 0;
2019 
2020  // Now begin NLO part for tree-level events
2021  if ( doOASTree || doOASSubt ) {
2022 
2023  // Decide on which order to expand to.
2024  int order = ( nSteps > 0 && nSteps <= nMaxNLO) ? 1 : -1;
2025 
2026  // Exclusive inputs:
2027  // Subtract only the O(\alpha_s^{n+0})-term from the tree-level
2028  // subtraction, if we're at the highest NLO multiplicity (nMaxNLO).
2029  if ( nloTilde && doUNLOPSSubt && nRecluster == 1
2030  && nSteps == nMaxNLO+1 ) order = 0;
2031 
2032  // Exclusive inputs:
2033  // Do not remove the O(as)-term if the number of reclusterings
2034  // exceeds the number of NLO jets, or if more clusterings have
2035  // been performed.
2036  if (nloTilde && doUNLOPSSubt && ( nSteps > nMaxNLO+1
2037  || (nSteps == nMaxNLO+1 && nPerformed != nRecluster) ))
2038  order = -1;
2039 
2040  // Calculate terms in expansion of the CKKW-L weight.
2041  wgtFIRST = myHistory->weight_UNLOPS_CORRECTION( order,
2042  trialPartonLevelPtr, mergingHooksPtr->AlphaS_FSR(),
2043  mergingHooksPtr->AlphaS_ISR(), mergingHooksPtr->AlphaEM_FSR(),
2044  mergingHooksPtr->AlphaEM_ISR(), RNpath, rndmPtr );
2045 
2046  // Exclusive inputs:
2047  // Subtract the O(\alpha_s^{n+1})-term from the tree-level
2048  // subtraction, not the O(\alpha_s^{n+0})-terms.
2049  if ( nloTilde && doUNLOPSSubt && nRecluster == 1
2050  && nPerformed == nRecluster && nSteps <= nMaxNLO )
2051  wgtFIRST += 1.;
2052 
2053  // Subtract the O(\alpha_s)-term from the CKKW-L weight
2054  // If PDF contributions have not been included, subtract these later
2055  // New UNLOPS based on UN2LOPS.
2056  if (doUNLOPS2 && order > -1) wgt = -wgt*(wgtFIRST-1.);
2057  else if (order > -1) wgt = wgt - wgtFIRST;
2058 
2059  }
2060 
2061  // If no-emission probability is zero.
2062  if ( allowReject && wgt == 0. ) return 0;
2063  //if ( wgt == 0. ) return 0;
2064 
2065  // Done
2066  return 1;
2067 
2068 }
2069 
2070 //--------------------------------------------------------------------------
2071 
2072 // Function to perform UNLOPS merging on this event.
2073 
2074 int DireMerging::getStartingConditions( double RNpath, Event& process) {
2075 
2076  // Initialise which part of UNLOPS merging is applied.
2077  bool doUNLOPSSubt = settingsPtr->flag("Merging:doUNLOPSSubt");
2078  bool doUNLOPSSubtNLO = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
2079  // Save number of looping steps
2080  mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
2081  int nRecluster = settingsPtr->mode("Merging:nRecluster");
2082 
2083  // Calculate number of clustering steps
2084  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( myHistory->state,
2085  true);
2086  int nRequested = mergingHooksPtr->nRequested();
2087 
2088  // Potentially recluster real emission jets for powheg input containing
2089  // "too many" jets, i.e. real-emission kinematics.
2090  bool containsRealKin = nSteps > nRequested && nSteps > 0;
2091  if ( containsRealKin ) nRecluster += nSteps - nRequested;
2092 
2093  // Event with production scales set for further (trial) showering
2094  // and starting conditions for the shower.
2095  int nPerformed = 0;
2096  if (!doUNLOPSSubt && !doUNLOPSSubtNLO && !containsRealKin )
2097  myHistory->getStartingConditions(RNpath, process);
2098  // Do reclustering (looping) steps.
2099  else myHistory->getFirstClusteredEventAboveTMS( RNpath, nRecluster, process,
2100  nPerformed, true );
2101 
2102  // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
2103  // --> Set to minimal mT of partons.
2104  int nFinal = 0;
2105  double muf = process[0].e();
2106  for ( int i=0; i < process.size(); ++i )
2107  if ( process[i].isFinal()
2108  && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
2109  nFinal++;
2110  muf = min( muf, abs(process[i].mT()) );
2111  }
2112  // For pure QCD dijet events (only!), set the process scale to the
2113  // transverse momentum of the outgoing partons.
2114  if ( nSteps == 0 && nFinal == 2
2115  && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
2116  || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
2117  process.scale(muf);
2118 
2119  // Reset hard process candidates (changed after clustering a parton).
2120  mergingHooksPtr->storeHardProcessCandidates( process );
2121 
2122  // Check if resonance structure has been changed
2123  // (e.g. because of clustering W/Z/gluino)
2124  vector <int> oldResonance;
2125  for ( int i=0; i < myHistory->state.size(); ++i )
2126  if ( myHistory->state[i].status() == 22 )
2127  oldResonance.push_back(myHistory->state[i].id());
2128  vector <int> newResonance;
2129  for ( int i=0; i < process.size(); ++i )
2130  if ( process[i].status() == 22 )
2131  newResonance.push_back(process[i].id());
2132  // Compare old and new resonances
2133  for ( int i=0; i < int(oldResonance.size()); ++i )
2134  for ( int j=0; j < int(newResonance.size()); ++j )
2135  if ( newResonance[j] == oldResonance[i] ) {
2136  oldResonance[i] = 99;
2137  break;
2138  }
2139  bool hasNewResonances = (newResonance.size() != oldResonance.size());
2140  for ( int i=0; i < int(oldResonance.size()); ++i )
2141  hasNewResonances = (oldResonance[i] != 99);
2142 
2143  // If necessary, reattach resonance decay products.
2144  if (!hasNewResonances) mergingHooksPtr->reattachResonanceDecays(process);
2145 
2146  // Allow merging hooks to remove emissions from now on.
2147  mergingHooksPtr->doIgnoreEmissions(false);
2148 
2149  // If the resonance structure of the process has changed due to reclustering,
2150  // redo the resonance decays in Pythia::next()
2151  if (hasNewResonances) return 2;
2152 
2153  // Done
2154  return 1;
2155 
2156 }
2157 
2158 //--------------------------------------------------------------------------
2159 
2160 // Function to apply the merging scale cut on an input event.
2161 
2162 bool DireMerging::cutOnProcess( Event& process) {
2163 
2164  // Save number of looping steps
2165  mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
2166 
2167  // For now, prefer construction of ordered histories.
2168  mergingHooksPtr->orderHistories(true);
2169  // For pp > h, allow cut on state, so that underlying processes
2170  // can be clustered to gg > h
2171  if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
2172  mergingHooksPtr->allowCutOnRecState(true);
2173 
2174  // Reset any incoming spins for W+-.
2175  if (mergingHooksPtr->doWeakClustering())
2176  for (int i = 0;i < process.size();++i)
2177  process[i].pol(9);
2178 
2179  // Prepare process record for merging. If Pythia has already decayed
2180  // resonances used to define the hard process, remove resonance decay
2181  // products.
2182  Event newProcess( mergingHooksPtr->bareEvent( process, true) );
2183  // Store candidates for the splitting V -> qqbar'
2184  mergingHooksPtr->storeHardProcessCandidates( newProcess );
2185 
2186  // Check if event passes the merging scale cut.
2187  double tmsval = mergingHooksPtr->tms();
2188  // Get merging scale in current event.
2189  double tmsnow = mergingHooksPtr->tmsNow( newProcess );
2190  // Calculate number of clustering steps
2191  int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess, true);
2192 
2193  // Too few steps can be possible if a chain of resonance decays has been
2194  // removed. In this case, reject this event, since it will be handled in
2195  // lower-multiplicity samples.
2196  int nRequested = mergingHooksPtr->nRequested();
2197  if (nSteps < nRequested) return true;
2198 
2199  // Reset the minimal tms value, if necessary.
2200  tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
2201 
2202  // Potentially recluster real emission jets for powheg input containing
2203  // "too many" jets, i.e. real-emission kinematics.
2204  bool containsRealKin = nSteps > nRequested && nSteps > 0;
2205 
2206  // Get random number to choose a path.
2207  double RN = rndmPtr->flat();
2208  // Set dummy process scale.
2209  newProcess.scale(0.0);
2210  // Generate all histories
2211  DireHistory FullHistory( nSteps, 0.0, newProcess, DireClustering(),
2212  mergingHooksPtr,
2213  (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
2214  trialPartonLevelPtr, fsr, isr, psweights, coupSMPtr, true, true,
2215  1.0, 1.0, 1.0, 1.0, 0);
2216  // Project histories onto desired branches, e.g. only ordered paths.
2217  FullHistory.projectOntoDesiredHistories();
2218 
2219  // Remove real emission events without underlying Born configuration from
2220  // the loop sample, since such states will be taken care of by tree-level
2221  // samples.
2222  if ( containsRealKin && !allowIncompleteReal
2223  && FullHistory.select(RN)->nClusterings() == 0 )
2224  return true;
2225 
2226  // Cut if no history passes the cut on the lowest-multiplicity state.
2227  double dampWeight = mergingHooksPtr->dampenIfFailCuts(
2228  FullHistory.lowestMultProc(RN) );
2229  if ( dampWeight == 0. ) return true;
2230 
2231  // Do not apply cut if the configuration could not be projected onto an
2232  // underlying born configuration.
2233  if ( nSteps > 0 && FullHistory.select(RN)->nClusterings() == 0 )
2234  return false;
2235 
2236  // Now enfore merging scale cut if the event did not pass the merging scale
2237  // criterion.
2238  if ( nSteps > 0 && nSteps == nRequested && tmsnow < tmsval && tmsval > 0.) {
2239  string message="Warning in DireMerging::cutOnProcess: Les Houches Event";
2240  message+=" fails merging scale cut. Reject event.";
2241  infoPtr->errorMsg(message);
2242  return true;
2243  }
2244 
2245  // Check if more steps should be taken.
2246  int nFinalP = 0;
2247  int nFinalW = 0;
2248  Event coreProcess = Event();
2249  coreProcess.clear();
2250  coreProcess.init( "(hard process-modified)", particleDataPtr );
2251  coreProcess.clear();
2252  coreProcess = FullHistory.lowestMultProc(RN);
2253  for ( int i = 0; i < coreProcess.size(); ++i )
2254  if ( coreProcess[i].isFinal() ) {
2255  if ( coreProcess[i].colType() != 0 )
2256  nFinalP++;
2257  if ( coreProcess[i].idAbs() == 24 )
2258  nFinalW++;
2259  }
2260 
2261  bool complete = (FullHistory.select(RN)->nClusterings() == nSteps) ||
2262  ( mergingHooksPtr->doWeakClustering() && nFinalP == 2 && nFinalW == 0 );
2263 
2264  if ( !complete ) {
2265  string message="Warning in DireMerging::cutOnProcess: No clusterings";
2266  message+=" found. History incomplete.";
2267  infoPtr->errorMsg(message);
2268  }
2269 
2270  // Done if no real-emission jets are present.
2271  if ( !containsRealKin ) return false;
2272 
2273  // Now cut on events that contain an additional real-emission jet.
2274  // Perform one reclustering for real emission kinematics, then apply merging
2275  // scale cut on underlying Born kinematics.
2276  Event dummy = Event();
2277  // Initialise temporary output of reclustering.
2278  dummy.clear();
2279  dummy.init( "(hard process-modified)", particleDataPtr );
2280  dummy.clear();
2281  // Recluster once.
2282  FullHistory.getClusteredEvent( RN, nSteps, dummy );
2283  double tnowNew = mergingHooksPtr->tmsNow( dummy );
2284  // Veto if underlying Born kinematics do not pass merging scale cut.
2285  if ( nRequested > 0 && tnowNew < tmsval && tmsval > 0.) {
2286  string message="Warning in DireMerging::cutOnProcess: Les Houches Event";
2287  message+=" fails merging scale cut. Reject event.";
2288  infoPtr->errorMsg(message);
2289  return true;
2290  }
2291 
2292  // Done if only interested in cross section estimate after cuts.
2293  return false;
2294 
2295 }
2296 
2297 //==========================================================================
2298 
2299 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5