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 Merging class.
8 :
9 : #include "Pythia8/Merging.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // The Merging class.
16 :
17 : //--------------------------------------------------------------------------
18 :
19 : // Factor by which the maximal value of the merging scale can deviate before
20 : // a warning is printed.
21 : const double Merging::TMSMISMATCH = 1.5;
22 :
23 : //--------------------------------------------------------------------------
24 :
25 : // Initialise Merging class
26 :
27 : void Merging::init( Settings* settingsPtrIn, Info* infoPtrIn,
28 : ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
29 : BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
30 : MergingHooks* mergingHooksPtrIn, PartonLevel* trialPartonLevelPtrIn ){
31 :
32 : // Save pointers.
33 0 : settingsPtr = settingsPtrIn;
34 0 : infoPtr = infoPtrIn;
35 0 : particleDataPtr = particleDataPtrIn;
36 0 : rndmPtr = rndmPtrIn;
37 0 : mergingHooksPtr = mergingHooksPtrIn;
38 0 : trialPartonLevelPtr = trialPartonLevelPtrIn;
39 0 : beamAPtr = beamAPtrIn;
40 0 : beamBPtr = beamBPtrIn;
41 : // Reset minimal tms value.
42 0 : tmsNowMin = infoPtr->eCM();
43 :
44 0 : }
45 :
46 : //--------------------------------------------------------------------------
47 :
48 : // Function to print information.
49 : void Merging::statistics( ostream& os ) {
50 :
51 : // Recall switch to enfore merging scale cut.
52 0 : bool enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
53 : // Recall merging scale value.
54 0 : double tmsval = mergingHooksPtr->tms();
55 0 : bool printBanner = enforceCutOnLHE && tmsNowMin > TMSMISMATCH*tmsval;
56 : // Reset minimal tms value.
57 0 : tmsNowMin = infoPtr->eCM();
58 :
59 0 : if (!printBanner) return;
60 :
61 : // Header.
62 0 : os << "\n *------- PYTHIA Matrix Element Merging Information ------"
63 0 : << "-------------------------------------------------------*\n"
64 0 : << " | "
65 0 : << " |\n";
66 : // Print warning if the minimal tms value of any event was significantly
67 : // above the desired merging scale value.
68 0 : os << " | Warning in Merging::statistics: All Les Houches events"
69 0 : << " significantly above Merging:TMS cut. Please check. |\n";
70 :
71 : // Listing finished.
72 0 : os << " | "
73 0 : << " |\n"
74 0 : << " *------- End PYTHIA Matrix Element Merging Information -----"
75 0 : << "-----------------------------------------------------*" << endl;
76 0 : }
77 :
78 : //--------------------------------------------------------------------------
79 :
80 : // Function to steer different merging prescriptions.
81 :
82 : int Merging::mergeProcess(Event& process){
83 :
84 : int vetoCode = 1;
85 :
86 : // Reinitialise hard process.
87 0 : mergingHooksPtr->hardProcess.clear();
88 0 : mergingHooksPtr->processSave = settingsPtr->word("Merging:Process");
89 0 : mergingHooksPtr->hardProcess.initOnProcess(
90 0 : settingsPtr->word("Merging:Process"), particleDataPtr);
91 :
92 0 : mergingHooksPtr->doUserMergingSave
93 0 : = settingsPtr->flag("Merging:doUserMerging");
94 0 : mergingHooksPtr->doMGMergingSave
95 0 : = settingsPtr->flag("Merging:doMGMerging");
96 0 : mergingHooksPtr->doKTMergingSave
97 0 : = settingsPtr->flag("Merging:doKTMerging");
98 0 : mergingHooksPtr->doPTLundMergingSave
99 0 : = settingsPtr->flag("Merging:doPTLundMerging");
100 0 : mergingHooksPtr->doCutBasedMergingSave
101 0 : = settingsPtr->flag("Merging:doCutBasedMerging");
102 0 : mergingHooksPtr->doNL3TreeSave
103 0 : = settingsPtr->flag("Merging:doNL3Tree");
104 0 : mergingHooksPtr->doNL3LoopSave
105 0 : = settingsPtr->flag("Merging:doNL3Loop");
106 0 : mergingHooksPtr->doNL3SubtSave
107 0 : = settingsPtr->flag("Merging:doNL3Subt");
108 0 : mergingHooksPtr->doUNLOPSTreeSave
109 0 : = settingsPtr->flag("Merging:doUNLOPSTree");
110 0 : mergingHooksPtr->doUNLOPSLoopSave
111 0 : = settingsPtr->flag("Merging:doUNLOPSLoop");
112 0 : mergingHooksPtr->doUNLOPSSubtSave
113 0 : = settingsPtr->flag("Merging:doUNLOPSSubt");
114 0 : mergingHooksPtr->doUNLOPSSubtNLOSave
115 0 : = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
116 0 : mergingHooksPtr->doUMEPSTreeSave
117 0 : = settingsPtr->flag("Merging:doUMEPSTree");
118 0 : mergingHooksPtr->doUMEPSSubtSave
119 0 : = settingsPtr->flag("Merging:doUMEPSSubt");
120 0 : mergingHooksPtr->nReclusterSave
121 0 : = settingsPtr->mode("Merging:nRecluster");
122 :
123 : // Possibility to apply merging scale to an input event.
124 0 : bool applyTMSCut = settingsPtr->flag("Merging:doXSectionEstimate");
125 0 : if ( applyTMSCut && cutOnProcess(process) ) return -1;
126 : // Done if only a cut should be applied.
127 0 : if ( applyTMSCut ) return 1;
128 :
129 : // Possibility to perform CKKW-L merging on this event.
130 0 : if ( mergingHooksPtr->doCKKWLMerging() )
131 0 : vetoCode = mergeProcessCKKWL(process);
132 :
133 : // Possibility to perform UMEPS merging on this event.
134 0 : if ( mergingHooksPtr->doUMEPSMerging() )
135 0 : vetoCode = mergeProcessUMEPS(process);
136 :
137 : // Possibility to perform NL3 NLO merging on this event.
138 0 : if ( mergingHooksPtr->doNL3Merging() )
139 0 : vetoCode = mergeProcessNL3(process);
140 :
141 : // Possibility to perform UNLOPS merging on this event.
142 0 : if ( mergingHooksPtr->doUNLOPSMerging() )
143 0 : vetoCode = mergeProcessUNLOPS(process);
144 :
145 0 : return vetoCode;
146 :
147 0 : }
148 :
149 : //--------------------------------------------------------------------------
150 :
151 : // Function to perform CKKW-L merging on this event.
152 :
153 : int Merging::mergeProcessCKKWL( Event& process) {
154 :
155 : // Ensure that merging hooks to not veto events in the trial showers.
156 0 : mergingHooksPtr->doIgnoreStep(true);
157 : // For pp > h, allow cut on state, so that underlying processes
158 : // can be clustered to gg > h
159 0 : if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0 )
160 0 : mergingHooksPtr->allowCutOnRecState(true);
161 : // For now, prefer construction of ordered histories.
162 0 : mergingHooksPtr->orderHistories(true);
163 :
164 : // Reset weight of the event.
165 : double wgt = 1.0;
166 0 : mergingHooksPtr->setWeightCKKWL(1.);
167 0 : mergingHooksPtr->muMI(-1.);
168 :
169 : // Prepare process record for merging. If Pythia has already decayed
170 : // resonances used to define the hard process, remove resonance decay
171 : // products.
172 0 : Event newProcess( mergingHooksPtr->bareEvent( process, true) );
173 : // Store candidates for the splitting V -> qqbar'.
174 0 : mergingHooksPtr->storeHardProcessCandidates( newProcess);
175 :
176 : // Check if event passes the merging scale cut.
177 0 : double tmsval = mergingHooksPtr->tms();
178 : // Get merging scale in current event.
179 0 : double tmsnow = mergingHooksPtr->tmsNow( newProcess );
180 : // Calculate number of clustering steps.
181 0 : int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess);
182 :
183 : // Too few steps can be possible if a chain of resonance decays has been
184 : // removed. In this case, reject this event, since it will be handled in
185 : // lower-multiplicity samples.
186 0 : int nRequested = settingsPtr->mode("Merging:nRequested");
187 0 : if (nSteps < nRequested) {
188 0 : mergingHooksPtr->setWeightCKKWL(0.);
189 0 : return -1;
190 : }
191 :
192 : // Reset the minimal tms value, if necessary.
193 0 : tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
194 :
195 : // Get random number to choose a path.
196 0 : double RN = rndmPtr->flat();
197 : // Set dummy process scale.
198 0 : newProcess.scale(0.0);
199 : // Generate all histories.
200 0 : History FullHistory( nSteps, 0.0, newProcess, Clustering(), mergingHooksPtr,
201 0 : (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr, true, true,
202 : true, true, 1.0, 0);
203 : // Project histories onto desired branches, e.g. only ordered paths.
204 0 : FullHistory.projectOntoDesiredHistories();
205 :
206 : // Do not apply cut if the configuration could not be projected onto an
207 : // underlying born configuration.
208 0 : bool applyCut = nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
209 :
210 : // Enfore merging scale cut if the event did not pass the merging scale
211 : // criterion.
212 0 : bool enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
213 0 : if ( enforceCutOnLHE && applyCut && tmsnow < tmsval ) {
214 0 : string message="Warning in Merging::mergeProcessCKKWL: Les Houches Event";
215 0 : message+=" fails merging scale cut. Reject event.";
216 0 : infoPtr->errorMsg(message);
217 0 : mergingHooksPtr->setWeightCKKWL(0.);
218 : return -1;
219 0 : }
220 :
221 0 : if ( FullHistory.select(RN)->nClusterings() < nSteps) {
222 0 : string message="Warning in Merging::mergeProcessCKKWL: No clusterings";
223 0 : message+=" found. History incomplete.";
224 0 : infoPtr->errorMsg(message);
225 0 : }
226 :
227 : // Calculate CKKWL weight:
228 : // Perform reweighting with Sudakov factors, save alpha_s ratios and
229 : // PDF ratio weights.
230 0 : wgt = FullHistory.weightTREE( trialPartonLevelPtr,
231 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
232 :
233 : // Event with production scales set for further (trial) showering
234 : // and starting conditions for the shower.
235 0 : FullHistory.getStartingConditions( RN, process );
236 : // If necessary, reattach resonance decay products.
237 0 : mergingHooksPtr->reattachResonanceDecays(process);
238 :
239 : // Allow to dampen histories in which the lowest multiplicity reclustered
240 : // state does not pass the lowest multiplicity cut of the matrix element.
241 0 : double dampWeight = mergingHooksPtr->dampenIfFailCuts(
242 0 : FullHistory.lowestMultProc(RN) );
243 : // Save the weight of the event for histogramming. Only change the
244 : // event weight after trial shower on the matrix element
245 : // multiplicity event (= in doVetoStep).
246 0 : wgt *= dampWeight;
247 :
248 : // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
249 : // --> Set to minimal mT of partons.
250 : int nFinal = 0;
251 0 : double muf = process[0].e();
252 0 : for ( int i=0; i < process.size(); ++i )
253 0 : if ( process[i].isFinal()
254 0 : && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
255 0 : nFinal++;
256 0 : muf = min( muf, abs(process[i].mT()) );
257 0 : }
258 : // For pure QCD dijet events (only!), set the process scale to the
259 : // transverse momentum of the outgoing partons.
260 0 : if ( nSteps == 0 && nFinal == 2
261 0 : && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
262 0 : || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
263 0 : process.scale(muf);
264 :
265 : // Save the weight of the event for histogramming.
266 0 : mergingHooksPtr->setWeightCKKWL(wgt);
267 :
268 : // Allow merging hooks to veto events from now on.
269 0 : mergingHooksPtr->doIgnoreStep(false);
270 :
271 : // If no-emission probability is zero.
272 0 : if ( wgt == 0. ) return 0;
273 :
274 : // Done
275 0 : return 1;
276 :
277 0 : }
278 :
279 : //--------------------------------------------------------------------------
280 :
281 : // Function to perform UMEPS merging on this event.
282 :
283 : int Merging::mergeProcessUMEPS( Event& process) {
284 :
285 : // Initialise which part of UMEPS merging is applied.
286 0 : bool doUMEPSTree = settingsPtr->flag("Merging:doUMEPSTree");
287 0 : bool doUMEPSSubt = settingsPtr->flag("Merging:doUMEPSSubt");
288 : // Save number of looping steps
289 0 : mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
290 0 : int nRecluster = settingsPtr->mode("Merging:nRecluster");
291 :
292 : // Ensure that merging hooks does not remove emissions.
293 0 : mergingHooksPtr->doIgnoreEmissions(true);
294 : // For pp > h, allow cut on state, so that underlying processes
295 : // can be clustered to gg > h
296 0 : if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0 )
297 0 : mergingHooksPtr->allowCutOnRecState(true);
298 : // For now, prefer construction of ordered histories.
299 0 : mergingHooksPtr->orderHistories(true);
300 :
301 : // Reset weights of the event.
302 : double wgt = 1.;
303 0 : mergingHooksPtr->setWeightCKKWL(1.);
304 0 : mergingHooksPtr->muMI(-1.);
305 :
306 : // Prepare process record for merging. If Pythia has already decayed
307 : // resonances used to define the hard process, remove resonance decay
308 : // products.
309 0 : Event newProcess( mergingHooksPtr->bareEvent( process, true) );
310 : // Store candidates for the splitting V -> qqbar'.
311 0 : mergingHooksPtr->storeHardProcessCandidates( newProcess );
312 :
313 : // Check if event passes the merging scale cut.
314 0 : double tmsval = mergingHooksPtr->tms();
315 : // Get merging scale in current event.
316 0 : double tmsnow = mergingHooksPtr->tmsNow( newProcess );
317 : // Calculate number of clustering steps.
318 0 : int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess );
319 :
320 : // Too few steps can be possible if a chain of resonance decays has been
321 : // removed. In this case, reject this event, since it will be handled in
322 : // lower-multiplicity samples.
323 0 : int nRequested = settingsPtr->mode("Merging:nRequested");
324 0 : if (nSteps < nRequested) {
325 0 : mergingHooksPtr->setWeightCKKWL(0.);
326 0 : return -1;
327 : }
328 :
329 : // Reset the minimal tms value, if necessary.
330 0 : tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
331 :
332 : // Get random number to choose a path.
333 0 : double RN = rndmPtr->flat();
334 : // Set dummy process scale.
335 0 : newProcess.scale(0.0);
336 : // Generate all histories.
337 0 : History FullHistory( nSteps, 0.0, newProcess, Clustering(), mergingHooksPtr,
338 0 : (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr, true, true,
339 : true, true, 1.0, 0);
340 : // Project histories onto desired branches, e.g. only ordered paths.
341 0 : FullHistory.projectOntoDesiredHistories();
342 :
343 : // Do not apply cut if the configuration could not be projected onto an
344 : // underlying born configuration.
345 0 : bool applyCut = nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
346 :
347 : // Enfore merging scale cut if the event did not pass the merging scale
348 : // criterion.
349 0 : bool enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
350 0 : if ( enforceCutOnLHE && applyCut && tmsnow < tmsval ) {
351 0 : string message="Warning in Merging::mergeProcessUMEPS: Les Houches Event";
352 0 : message+=" fails merging scale cut. Reject event.";
353 0 : infoPtr->errorMsg(message);
354 0 : mergingHooksPtr->setWeightCKKWL(0.);
355 : return -1;
356 0 : }
357 :
358 : // Check reclustering steps to correctly apply MPI.
359 0 : int nPerformed = 0;
360 0 : if ( nSteps > 0 && doUMEPSSubt
361 0 : && !FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster,
362 : newProcess, nPerformed, false ) ) {
363 : // Discard if the state could not be reclustered to a state above TMS.
364 0 : mergingHooksPtr->setWeightCKKWL(0.);
365 0 : return -1;
366 : }
367 :
368 0 : mergingHooksPtr->nMinMPI(nSteps - nPerformed);
369 :
370 : // Calculate CKKWL weight:
371 : // Perform reweighting with Sudakov factors, save alpha_s ratios and
372 : // PDF ratio weights.
373 0 : if ( doUMEPSTree ) {
374 0 : wgt = FullHistory.weight_UMEPS_TREE( trialPartonLevelPtr,
375 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN );
376 0 : } else {
377 0 : wgt = FullHistory.weight_UMEPS_SUBT( trialPartonLevelPtr,
378 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN );
379 : }
380 :
381 : // Event with production scales set for further (trial) showering
382 : // and starting conditions for the shower.
383 0 : if ( doUMEPSTree ) FullHistory.getStartingConditions( RN, process );
384 : // Do reclustering (looping) steps.
385 0 : else FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster, process,
386 : nPerformed, true );
387 :
388 : // Allow to dampen histories in which the lowest multiplicity reclustered
389 : // state does not pass the lowest multiplicity cut of the matrix element
390 0 : double dampWeight = mergingHooksPtr->dampenIfFailCuts(
391 0 : FullHistory.lowestMultProc(RN) );
392 : // Save the weight of the event for histogramming. Only change the
393 : // event weight after trial shower on the matrix element
394 : // multiplicity event (= in doVetoStep)
395 0 : wgt *= dampWeight;
396 :
397 : // Save the weight of the event for histogramming.
398 0 : mergingHooksPtr->setWeightCKKWL(wgt);
399 :
400 : // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
401 : // --> Set to minimal mT of partons.
402 : int nFinal = 0;
403 0 : double muf = process[0].e();
404 0 : for ( int i=0; i < process.size(); ++i )
405 0 : if ( process[i].isFinal()
406 0 : && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
407 0 : nFinal++;
408 0 : muf = min( muf, abs(process[i].mT()) );
409 0 : }
410 :
411 : // For pure QCD dijet events (only!), set the process scale to the
412 : // transverse momentum of the outgoing partons.
413 : // Calculate number of clustering steps.
414 0 : int nStepsNew = mergingHooksPtr->getNumberOfClusteringSteps( process );
415 0 : if ( nStepsNew == 0
416 0 : && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
417 0 : || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
418 0 : process.scale(muf);
419 :
420 : // Reset hard process candidates (changed after clustering a parton).
421 0 : mergingHooksPtr->storeHardProcessCandidates( process );
422 : // If necessary, reattach resonance decay products.
423 0 : mergingHooksPtr->reattachResonanceDecays(process);
424 :
425 : // Allow merging hooks to remove emissions from now on.
426 0 : mergingHooksPtr->doIgnoreEmissions(false);
427 :
428 : // If no-emission probability is zero.
429 0 : if ( wgt == 0. ) return 0;
430 :
431 : // Done
432 0 : return 1;
433 :
434 0 : }
435 :
436 : //--------------------------------------------------------------------------
437 :
438 : // Function to perform NL3 NLO merging on this event.
439 :
440 : int Merging::mergeProcessNL3( Event& process) {
441 :
442 : // Initialise which part of NL3 merging is applied.
443 0 : bool doNL3Tree = settingsPtr->flag("Merging:doNL3Tree");
444 0 : bool doNL3Loop = settingsPtr->flag("Merging:doNL3Loop");
445 0 : bool doNL3Subt = settingsPtr->flag("Merging:doNL3Subt");
446 : // Save number of looping steps.
447 0 : int nRequested = settingsPtr->mode("Merging:nRequested");
448 :
449 : // Ensure that hooks (NL3 part) to not remove emissions.
450 0 : mergingHooksPtr->doIgnoreEmissions(true);
451 : // Ensure that hooks (CKKWL part) to not veto events in trial showers.
452 0 : mergingHooksPtr->doIgnoreStep(true);
453 : // For pp > h, allow cut on state, so that underlying processes
454 : // can be clustered to gg > h
455 0 : if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
456 0 : mergingHooksPtr->allowCutOnRecState(true);
457 : // For now, prefer construction of ordered histories.
458 0 : mergingHooksPtr->orderHistories(true);
459 :
460 : // Reset weight of the event
461 : double wgt = 1.;
462 0 : mergingHooksPtr->setWeightCKKWL(1.);
463 : // Reset the O(alphaS)-term of the CKKW-L weight.
464 : double wgtFIRST = 0.;
465 0 : mergingHooksPtr->setWeightFIRST(0.);
466 0 : mergingHooksPtr->muMI(-1.);
467 :
468 : // Prepare process record for merging. If Pythia has already decayed
469 : // resonances used to define the hard process, remove resonance decay
470 : // products.
471 0 : Event newProcess( mergingHooksPtr->bareEvent( process, true) );
472 : // Store candidates for the splitting V -> qqbar'
473 0 : mergingHooksPtr->storeHardProcessCandidates( newProcess);
474 :
475 : // Check if event passes the merging scale cut.
476 0 : double tmsval = mergingHooksPtr->tms();
477 : // Get merging scale in current event.
478 0 : double tmsnow = mergingHooksPtr->tmsNow( newProcess );
479 : // Calculate number of clustering steps
480 0 : int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess);
481 :
482 : // Too few steps can be possible if a chain of resonance decays has been
483 : // removed. In this case, reject this event, since it will be handled in
484 : // lower-multiplicity samples.
485 0 : if (nSteps < nRequested) {
486 0 : mergingHooksPtr->setWeightCKKWL(0.);
487 0 : mergingHooksPtr->setWeightFIRST(0.);
488 0 : return -1;
489 : }
490 :
491 : // Reset the minimal tms value, if necessary.
492 0 : tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
493 :
494 : // Enfore merging scale cut if the event did not pass the merging scale
495 : // criterion.
496 0 : bool enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
497 0 : if ( enforceCutOnLHE && nSteps > 0 && nSteps == nRequested
498 0 : && tmsnow < tmsval ) {
499 0 : string message="Warning in Merging::mergeProcessNL3: Les Houches Event";
500 0 : message+=" fails merging scale cut. Reject event.";
501 0 : infoPtr->errorMsg(message);
502 0 : mergingHooksPtr->setWeightCKKWL(0.);
503 0 : mergingHooksPtr->setWeightFIRST(0.);
504 : return -1;
505 0 : }
506 :
507 : // Get random number to choose a path.
508 0 : double RN = rndmPtr->flat();
509 : // Set dummy process scale.
510 0 : newProcess.scale(0.0);
511 : // Generate all histories
512 0 : History FullHistory( nSteps, 0.0, newProcess, Clustering(), mergingHooksPtr,
513 0 : (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr, true, true,
514 : true, true, 1.0, 0);
515 : // Project histories onto desired branches, e.g. only ordered paths.
516 0 : FullHistory.projectOntoDesiredHistories();
517 :
518 : // Discard states that cannot be projected unto a state with one less jet.
519 0 : if ( nSteps > 0 && doNL3Subt
520 0 : && FullHistory.select(RN)->nClusterings() == 0 ){
521 0 : mergingHooksPtr->setWeightCKKWL(0.);
522 0 : mergingHooksPtr->setWeightFIRST(0.);
523 0 : return -1;
524 : }
525 :
526 : // Potentially recluster real emission jets for powheg input containing
527 : // "too many" jets, i.e. real-emission kinematics.
528 0 : bool containsRealKin = nSteps > nRequested && nSteps > 0;
529 :
530 : // Perform one reclustering for real emission kinematics, then apply merging
531 : // scale cut on underlying Born kinematics.
532 0 : if ( containsRealKin ) {
533 0 : Event dummy = Event();
534 : // Initialise temporary output of reclustering.
535 0 : dummy.clear();
536 0 : dummy.init( "(hard process-modified)", particleDataPtr );
537 0 : dummy.clear();
538 : // Recluster once.
539 0 : if ( !FullHistory.getClusteredEvent( RN, nSteps, dummy )) {
540 0 : mergingHooksPtr->setWeightCKKWL(0.);
541 0 : mergingHooksPtr->setWeightFIRST(0.);
542 0 : return -1;
543 : }
544 0 : double tnowNew = mergingHooksPtr->tmsNow( dummy );
545 : // Veto if underlying Born kinematics do not pass merging scale cut.
546 0 : if ( enforceCutOnLHE && nSteps > 0 && nRequested > 0
547 0 : && tnowNew < tmsval ) {
548 0 : mergingHooksPtr->setWeightCKKWL(0.);
549 0 : mergingHooksPtr->setWeightFIRST(0.);
550 0 : return -1;
551 : }
552 0 : }
553 :
554 : // Remember number of jets, to include correct MPI no-emission probabilities.
555 0 : if ( doNL3Subt || containsRealKin ) mergingHooksPtr->nMinMPI(nSteps - 1);
556 0 : else mergingHooksPtr->nMinMPI(nSteps);
557 :
558 : // Calculate weight
559 : // Do LO or first part of NLO tree-level reweighting
560 0 : if( doNL3Tree ) {
561 : // Perform reweighting with Sudakov factors, save as ratios and
562 : // PDF ratio weights
563 0 : wgt = FullHistory.weightTREE( trialPartonLevelPtr,
564 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
565 0 : } else if( doNL3Loop || doNL3Subt ) {
566 : // No reweighting, just set event scales properly and incorporate MPI
567 : // no-emission probabilities.
568 0 : wgt = FullHistory.weightLOOP( trialPartonLevelPtr, RN);
569 0 : }
570 :
571 : // Event with production scales set for further (trial) showering
572 : // and starting conditions for the shower
573 0 : if ( !doNL3Subt && !containsRealKin )
574 0 : FullHistory.getStartingConditions(RN, process);
575 : // For sutraction of nSteps-additional resolved partons from
576 : // the nSteps-1 parton phase space, recluster the last parton
577 : // in nSteps-parton events, and sutract later
578 : else {
579 : // Function to return the reclustered event
580 0 : if ( !FullHistory.getClusteredEvent( RN, nSteps, process )) {
581 0 : mergingHooksPtr->setWeightCKKWL(0.);
582 0 : mergingHooksPtr->setWeightFIRST(0.);
583 0 : return -1;
584 : }
585 : }
586 :
587 : // Allow to dampen histories in which the lowest multiplicity reclustered
588 : // state does not pass the lowest multiplicity cut of the matrix element
589 0 : double dampWeight = mergingHooksPtr->dampenIfFailCuts(
590 0 : FullHistory.lowestMultProc(RN) );
591 : // Save the weight of the event for histogramming. Only change the
592 : // event weight after trial shower on the matrix element
593 : // multiplicity event (= in doVetoStep)
594 0 : wgt *= dampWeight;
595 :
596 : // For tree level samples in NL3, rescale with k-Factor
597 0 : if (doNL3Tree ){
598 : // Find k-factor
599 : double kFactor = 1.;
600 0 : if( nSteps > mergingHooksPtr->nMaxJetsNLO() )
601 0 : kFactor = mergingHooksPtr->kFactor( mergingHooksPtr->nMaxJetsNLO() );
602 0 : else kFactor = mergingHooksPtr->kFactor(nSteps);
603 : // For NLO merging, rescale CKKW-L weight with k-factor
604 0 : wgt *= kFactor;
605 0 : }
606 :
607 : // Save the weight of the event for histogramming
608 0 : mergingHooksPtr->setWeightCKKWL(wgt);
609 :
610 : // Check if we need to subtract the O(\alpha_s)-term. If the number
611 : // of additional partons is larger than the number of jets for
612 : // which loop matrix elements are available, do standard CKKW-L
613 0 : bool doOASTree = doNL3Tree && nSteps <= mergingHooksPtr->nMaxJetsNLO();
614 :
615 : // Now begin NLO part for tree-level events
616 0 : if ( doOASTree ) {
617 : // Calculate the O(\alpha_s)-term of the CKKWL weight
618 0 : wgtFIRST = FullHistory.weightFIRST( trialPartonLevelPtr,
619 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN,
620 0 : rndmPtr );
621 : // If necessary, also dampen the O(\alpha_s)-term
622 0 : wgtFIRST *= dampWeight;
623 : // Set the subtractive weight to the value calculated so far
624 0 : mergingHooksPtr->setWeightFIRST(wgtFIRST);
625 : // Subtract the O(\alpha_s)-term from the CKKW-L weight
626 : // If PDF contributions have not been included, subtract these later
627 : wgt = wgt - wgtFIRST;
628 0 : }
629 :
630 : // Set qcd 2->2 starting scale different from arbirtrary scale in LHEF!
631 : // --> Set to pT of partons
632 : double pT = 0.;
633 0 : for( int i=0; i < process.size(); ++i)
634 0 : if(process[i].isFinal() && process[i].colType() != 0) {
635 0 : pT = sqrt(pow(process[i].px(),2) + pow(process[i].py(),2));
636 0 : break;
637 : }
638 : // For pure QCD dijet events (only!), set the process scale to the
639 : // transverse momentum of the outgoing partons.
640 0 : if ( nSteps == 0
641 0 : && mergingHooksPtr->getProcessString().compare("pp>jj") == 0)
642 0 : process.scale(pT);
643 :
644 : // Reset hard process candidates (changed after clustering a parton).
645 0 : mergingHooksPtr->storeHardProcessCandidates( process );
646 : // If necessary, reattach resonance decay products.
647 0 : mergingHooksPtr->reattachResonanceDecays(process);
648 :
649 : // Allow merging hooks (NL3 part) to remove emissions from now on.
650 0 : mergingHooksPtr->doIgnoreEmissions(false);
651 : // Allow merging hooks (CKKWL part) to veto events from now on.
652 0 : mergingHooksPtr->doIgnoreStep(false);
653 :
654 : // Done
655 : return 1;
656 :
657 0 : }
658 :
659 : //--------------------------------------------------------------------------
660 :
661 : // Function to perform UNLOPS merging on this event.
662 :
663 : int Merging::mergeProcessUNLOPS( Event& process) {
664 :
665 : // Initialise which part of UNLOPS merging is applied.
666 0 : bool nloTilde = settingsPtr->flag("Merging:doUNLOPSTilde");
667 0 : bool doUNLOPSTree = settingsPtr->flag("Merging:doUNLOPSTree");
668 0 : bool doUNLOPSLoop = settingsPtr->flag("Merging:doUNLOPSLoop");
669 0 : bool doUNLOPSSubt = settingsPtr->flag("Merging:doUNLOPSSubt");
670 0 : bool doUNLOPSSubtNLO = settingsPtr->flag("Merging:doUNLOPSSubtNLO");
671 : // Save number of looping steps
672 0 : mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
673 0 : int nRecluster = settingsPtr->mode("Merging:nRecluster");
674 0 : int nRequested = settingsPtr->mode("Merging:nRequested");
675 :
676 : // Ensure that merging hooks to not remove emissions
677 0 : mergingHooksPtr->doIgnoreEmissions(true);
678 : // For now, prefer construction of ordered histories.
679 0 : mergingHooksPtr->orderHistories(true);
680 : // For pp > h, allow cut on state, so that underlying processes
681 : // can be clustered to gg > h
682 0 : if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
683 0 : mergingHooksPtr->allowCutOnRecState(true);
684 :
685 : // Reset weight of the event.
686 : double wgt = 1.;
687 0 : mergingHooksPtr->setWeightCKKWL(1.);
688 : // Reset the O(alphaS)-term of the UMEPS weight.
689 : double wgtFIRST = 0.;
690 0 : mergingHooksPtr->setWeightFIRST(0.);
691 0 : mergingHooksPtr->muMI(-1.);
692 :
693 : // Prepare process record for merging. If Pythia has already decayed
694 : // resonances used to define the hard process, remove resonance decay
695 : // products.
696 0 : Event newProcess( mergingHooksPtr->bareEvent( process, true) );
697 : // Store candidates for the splitting V -> qqbar'
698 0 : mergingHooksPtr->storeHardProcessCandidates( newProcess );
699 :
700 : // Check if event passes the merging scale cut.
701 0 : double tmsval = mergingHooksPtr->tms();
702 : // Get merging scale in current event.
703 0 : double tmsnow = mergingHooksPtr->tmsNow( newProcess );
704 : // Calculate number of clustering steps
705 0 : int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess);
706 :
707 : // Too few steps can be possible if a chain of resonance decays has been
708 : // removed. In this case, reject this event, since it will be handled in
709 : // lower-multiplicity samples.
710 0 : if (nSteps < nRequested) {
711 0 : string message="Warning in Merging::mergeProcessUNLOPS: Les Houches Event";
712 0 : message+=" after removing decay products does not contain enough partons.";
713 0 : infoPtr->errorMsg(message);
714 0 : mergingHooksPtr->setWeightCKKWL(0.);
715 0 : mergingHooksPtr->setWeightFIRST(0.);
716 : return -1;
717 0 : }
718 :
719 : // Reset the minimal tms value, if necessary.
720 0 : tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
721 :
722 : // Get random number to choose a path.
723 0 : double RN = rndmPtr->flat();
724 : // Set dummy process scale.
725 0 : newProcess.scale(0.0);
726 : // Generate all histories
727 0 : History FullHistory( nSteps, 0.0, newProcess, Clustering(), mergingHooksPtr,
728 0 : (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr, true, true,
729 : true, true, 1.0, 0);
730 : // Project histories onto desired branches, e.g. only ordered paths.
731 0 : FullHistory.projectOntoDesiredHistories();
732 :
733 : // Do not apply cut if the configuration could not be projected onto an
734 : // underlying born configuration.
735 0 : bool applyCut = nSteps > 0 && FullHistory.select(RN)->nClusterings() > 0;
736 :
737 : // Enfore merging scale cut if the event did not pass the merging scale
738 : // criterion.
739 0 : bool enforceCutOnLHE = settingsPtr->flag("Merging:enforceCutOnLHE");
740 0 : if ( enforceCutOnLHE && applyCut && nSteps == nRequested
741 0 : && tmsnow < tmsval ) {
742 0 : string message="Warning in Merging::mergeProcessUNLOPS: Les Houches Event";
743 0 : message+=" fails merging scale cut. Reject event.";
744 0 : infoPtr->errorMsg(message);
745 0 : mergingHooksPtr->setWeightCKKWL(0.);
746 0 : mergingHooksPtr->setWeightFIRST(0.);
747 : return -1;
748 0 : }
749 :
750 : // Potentially recluster real emission jets for powheg input containing
751 : // "too many" jets, i.e. real-emission kinematics.
752 0 : bool containsRealKin = nSteps > nRequested && nSteps > 0;
753 0 : if ( containsRealKin ) nRecluster += nSteps - nRequested;
754 :
755 : // Remove real emission events without underlying Born configuration from
756 : // the loop sample, since such states will be taken care of by tree-level
757 : // samples.
758 0 : bool allowIncompleteReal =
759 0 : settingsPtr->flag("Merging:allowIncompleteHistoriesInReal");
760 0 : if ( doUNLOPSLoop && containsRealKin && !allowIncompleteReal
761 0 : && FullHistory.select(RN)->nClusterings() == 0 ) {
762 0 : mergingHooksPtr->setWeightCKKWL(0.);
763 0 : mergingHooksPtr->setWeightFIRST(0.);
764 0 : return -1;
765 : }
766 :
767 : // Discard if the state could not be reclustered to any state above TMS.
768 0 : int nPerformed = 0;
769 0 : if ( nSteps > 0 && !allowIncompleteReal
770 0 : && ( doUNLOPSSubt || doUNLOPSSubtNLO || containsRealKin )
771 0 : && !FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster,
772 : newProcess, nPerformed, false ) ) {
773 0 : mergingHooksPtr->setWeightCKKWL(0.);
774 0 : mergingHooksPtr->setWeightFIRST(0.);
775 0 : return -1;
776 : }
777 : // Check reclustering steps to correctly apply MPI.
778 0 : mergingHooksPtr->nMinMPI(nSteps - nPerformed);
779 :
780 : // Perform one reclustering for real emission kinematics, then apply
781 : // merging scale cut on underlying Born kinematics.
782 0 : if ( containsRealKin ) {
783 0 : Event dummy = Event();
784 : // Initialise temporary output of reclustering.
785 0 : dummy.clear();
786 0 : dummy.init( "(hard process-modified)", particleDataPtr );
787 0 : dummy.clear();
788 : // Recluster once.
789 0 : FullHistory.getClusteredEvent( RN, nSteps, dummy );
790 0 : double tnowNew = mergingHooksPtr->tmsNow( dummy );
791 : // Veto if underlying Born kinematics do not pass merging scale cut.
792 0 : if ( enforceCutOnLHE && nSteps > 0 && nRequested > 0
793 0 : && tnowNew < tmsval ) {
794 0 : mergingHooksPtr->setWeightCKKWL(0.);
795 0 : mergingHooksPtr->setWeightFIRST(0.);
796 0 : return -1;
797 : }
798 0 : }
799 :
800 : // Calculate weights.
801 : // Do LO or first part of NLO tree-level reweighting
802 0 : if( doUNLOPSTree ) {
803 : // Perform reweighting with Sudakov factors, save as ratios and
804 : // PDF ratio weights
805 0 : wgt = FullHistory.weight_UNLOPS_TREE( trialPartonLevelPtr,
806 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
807 0 : } else if( doUNLOPSLoop ) {
808 : // No reweighting, just set event scales properly
809 0 : wgt = FullHistory.weight_UNLOPS_LOOP( trialPartonLevelPtr, RN);
810 0 : } else if( doUNLOPSSubtNLO ) {
811 : // The standard prescripition contains no real-emission parts
812 : // No reweighting, just set event scales properly
813 0 : wgt = FullHistory.weight_UNLOPS_SUBTNLO( trialPartonLevelPtr, RN);
814 0 : } else if( doUNLOPSSubt ) {
815 : // The standard prescripition contains no subtraction parts
816 : // No reweighting, just set event scales properly
817 0 : wgt = FullHistory.weight_UNLOPS_SUBT( trialPartonLevelPtr,
818 0 : mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
819 0 : }
820 :
821 : // Event with production scales set for further (trial) showering
822 : // and starting conditions for the shower.
823 0 : if (!doUNLOPSSubt && !doUNLOPSSubtNLO && !containsRealKin )
824 0 : FullHistory.getStartingConditions(RN, process);
825 : // Do reclustering (looping) steps.
826 0 : else FullHistory.getFirstClusteredEventAboveTMS( RN, nRecluster, process,
827 : nPerformed, true );
828 :
829 : // Allow to dampen histories in which the lowest multiplicity reclustered
830 : // state does not pass the lowest multiplicity cut of the matrix element
831 0 : double dampWeight = mergingHooksPtr->dampenIfFailCuts(
832 0 : FullHistory.lowestMultProc(RN) );
833 : // Save the weight of the event for histogramming. Only change the
834 : // event weight after trial shower on the matrix element
835 : // multiplicity event (= in doVetoStep)
836 0 : wgt *= dampWeight;
837 :
838 : // For tree-level or subtractive sammples, rescale with k-Factor
839 0 : if ( doUNLOPSTree || doUNLOPSSubt ){
840 : // Find k-factor
841 : double kFactor = 1.;
842 0 : if ( nSteps > mergingHooksPtr->nMaxJetsNLO() )
843 0 : kFactor = mergingHooksPtr->kFactor( mergingHooksPtr->nMaxJetsNLO() );
844 0 : else kFactor = mergingHooksPtr->kFactor(nSteps);
845 : // For NLO merging, rescale CKKW-L weight with k-factor
846 0 : wgt *= (nRecluster == 2 && nloTilde) ? 1. : kFactor;
847 0 : }
848 :
849 : // Save the weight of the event for histogramming
850 0 : mergingHooksPtr->setWeightCKKWL(wgt);
851 :
852 : // Check if we need to subtract the O(\alpha_s)-term. If the number
853 : // of additional partons is larger than the number of jets for
854 : // which loop matrix elements are available, do standard UMEPS.
855 0 : int nMaxNLO = mergingHooksPtr->nMaxJetsNLO();
856 0 : bool doOASTree = doUNLOPSTree && nSteps <= nMaxNLO;
857 0 : bool doOASSubt = doUNLOPSSubt && nSteps <= nMaxNLO+1 && nSteps > 0;
858 :
859 : // Now begin NLO part for tree-level events
860 0 : if ( doOASTree || doOASSubt ) {
861 :
862 : // Decide on which order to expand to.
863 0 : int order = ( nSteps > 0 && nSteps <= nMaxNLO) ? 1 : -1;
864 :
865 : // Exclusive inputs:
866 : // Subtract only the O(\alpha_s^{n+0})-term from the tree-level
867 : // subtraction, if we're at the highest NLO multiplicity (nMaxNLO).
868 0 : if ( nloTilde && doUNLOPSSubt && nRecluster == 1
869 0 : && nSteps == nMaxNLO+1 ) order = 0;
870 :
871 : // Exclusive inputs:
872 : // Do not remove the O(as)-term if the number of reclusterings
873 : // exceeds the number of NLO jets, or if more clusterings have
874 : // been performed.
875 0 : if (nloTilde && doUNLOPSSubt && ( nSteps > nMaxNLO+1
876 0 : || (nSteps == nMaxNLO+1 && nPerformed != nRecluster) ))
877 0 : order = -1;
878 :
879 : // Calculate terms in expansion of the CKKW-L weight.
880 0 : wgtFIRST = FullHistory.weight_UNLOPS_CORRECTION( order,
881 0 : trialPartonLevelPtr, mergingHooksPtr->AlphaS_FSR(),
882 0 : mergingHooksPtr->AlphaS_ISR(), RN, rndmPtr );
883 :
884 : // Exclusive inputs:
885 : // Subtract the O(\alpha_s^{n+1})-term from the tree-level
886 : // subtraction, not the O(\alpha_s^{n+0})-terms.
887 0 : if ( nloTilde && doUNLOPSSubt && nRecluster == 1
888 0 : && nPerformed == nRecluster && nSteps <= nMaxNLO )
889 0 : wgtFIRST += 1.;
890 :
891 : // If necessary, also dampen the O(\alpha_s)-term
892 0 : wgtFIRST *= dampWeight;
893 : // Set the subtractive weight to the value calculated so far
894 0 : mergingHooksPtr->setWeightFIRST(wgtFIRST);
895 : // Subtract the O(\alpha_s)-term from the CKKW-L weight
896 : // If PDF contributions have not been included, subtract these later
897 0 : wgt = wgt - wgtFIRST;
898 :
899 0 : }
900 :
901 : // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
902 : // --> Set to minimal mT of partons.
903 : int nFinal = 0;
904 0 : double muf = process[0].e();
905 0 : for ( int i=0; i < process.size(); ++i )
906 0 : if ( process[i].isFinal()
907 0 : && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
908 0 : nFinal++;
909 0 : muf = min( muf, abs(process[i].mT()) );
910 0 : }
911 : // For pure QCD dijet events (only!), set the process scale to the
912 : // transverse momentum of the outgoing partons.
913 0 : if ( nSteps == 0 && nFinal == 2
914 0 : && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
915 0 : || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
916 0 : process.scale(muf);
917 :
918 : // Reset hard process candidates (changed after clustering a parton).
919 0 : mergingHooksPtr->storeHardProcessCandidates( process );
920 :
921 : // Check if resonance structure has been changed
922 : // (e.g. because of clustering W/Z/gluino)
923 0 : vector <int> oldResonance;
924 0 : for ( int i=0; i < newProcess.size(); ++i )
925 0 : if ( newProcess[i].status() == 22 )
926 0 : oldResonance.push_back(newProcess[i].id());
927 0 : vector <int> newResonance;
928 0 : for ( int i=0; i < process.size(); ++i )
929 0 : if ( process[i].status() == 22 )
930 0 : newResonance.push_back(process[i].id());
931 : // Compare old and new resonances
932 0 : for ( int i=0; i < int(oldResonance.size()); ++i )
933 0 : for ( int j=0; j < int(newResonance.size()); ++j )
934 0 : if ( newResonance[j] == oldResonance[i] ) {
935 0 : oldResonance[i] = 99;
936 0 : break;
937 : }
938 0 : bool hasNewResonances = (newResonance.size() != oldResonance.size());
939 0 : for ( int i=0; i < int(oldResonance.size()); ++i )
940 0 : hasNewResonances = (oldResonance[i] != 99);
941 :
942 : // If necessary, reattach resonance decay products.
943 0 : if (!hasNewResonances) mergingHooksPtr->reattachResonanceDecays(process);
944 :
945 : // Allow merging hooks to remove emissions from now on.
946 0 : mergingHooksPtr->doIgnoreEmissions(false);
947 :
948 : // If no-emission probability is zero.
949 0 : if ( wgt == 0. ) return 0;
950 :
951 : // If the resonance structure of the process has changed due to reclustering,
952 : // redo the resonance decays in Pythia::next()
953 0 : if (hasNewResonances) return 2;
954 :
955 : // Done
956 0 : return 1;
957 :
958 0 : }
959 :
960 : //--------------------------------------------------------------------------
961 :
962 : // Function to apply the merging scale cut on an input event.
963 :
964 : bool Merging::cutOnProcess( Event& process) {
965 :
966 : // Save number of looping steps
967 0 : mergingHooksPtr->nReclusterSave = settingsPtr->mode("Merging:nRecluster");
968 0 : int nRequested = settingsPtr->mode("Merging:nRequested");
969 :
970 : // For now, prefer construction of ordered histories.
971 0 : mergingHooksPtr->orderHistories(true);
972 : // For pp > h, allow cut on state, so that underlying processes
973 : // can be clustered to gg > h
974 0 : if ( mergingHooksPtr->getProcessString().compare("pp>h") == 0)
975 0 : mergingHooksPtr->allowCutOnRecState(true);
976 :
977 : // Prepare process record for merging. If Pythia has already decayed
978 : // resonances used to define the hard process, remove resonance decay
979 : // products.
980 0 : Event newProcess( mergingHooksPtr->bareEvent( process, true) );
981 : // Store candidates for the splitting V -> qqbar'
982 0 : mergingHooksPtr->storeHardProcessCandidates( newProcess );
983 :
984 : // Check if event passes the merging scale cut.
985 0 : double tmsval = mergingHooksPtr->tms();
986 : // Get merging scale in current event.
987 0 : double tmsnow = mergingHooksPtr->tmsNow( newProcess );
988 : // Calculate number of clustering steps
989 0 : int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( newProcess);
990 :
991 : // Too few steps can be possible if a chain of resonance decays has been
992 : // removed. In this case, reject this event, since it will be handled in
993 : // lower-multiplicity samples.
994 0 : if (nSteps < nRequested) return true;
995 :
996 : // Reset the minimal tms value, if necessary.
997 0 : tmsNowMin = (nSteps == 0) ? 0. : min(tmsNowMin, tmsnow);
998 :
999 : // Potentially recluster real emission jets for powheg input containing
1000 : // "too many" jets, i.e. real-emission kinematics.
1001 0 : bool containsRealKin = nSteps > nRequested && nSteps > 0;
1002 :
1003 : // Get random number to choose a path.
1004 0 : double RN = rndmPtr->flat();
1005 : // Set dummy process scale.
1006 0 : newProcess.scale(0.0);
1007 : // Generate all histories
1008 0 : History FullHistory( nSteps, 0.0, newProcess, Clustering(), mergingHooksPtr,
1009 0 : (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr, true, true,
1010 : true, true, 1.0, 0);
1011 : // Project histories onto desired branches, e.g. only ordered paths.
1012 0 : FullHistory.projectOntoDesiredHistories();
1013 :
1014 : // Remove real emission events without underlying Born configuration from
1015 : // the loop sample, since such states will be taken care of by tree-level
1016 : // samples.
1017 0 : bool allowIncompleteReal =
1018 0 : settingsPtr->flag("Merging:allowIncompleteHistoriesInReal");
1019 0 : if ( containsRealKin && !allowIncompleteReal
1020 0 : && FullHistory.select(RN)->nClusterings() == 0 )
1021 0 : return true;
1022 :
1023 : // Cut if no history passes the cut on the lowest-multiplicity state.
1024 0 : double dampWeight = mergingHooksPtr->dampenIfFailCuts(
1025 0 : FullHistory.lowestMultProc(RN) );
1026 0 : if ( dampWeight == 0. ) return true;
1027 :
1028 : // Do not apply cut if the configuration could not be projected onto an
1029 : // underlying born configuration.
1030 0 : if ( nSteps > 0 && FullHistory.select(RN)->nClusterings() == 0 )
1031 0 : return false;
1032 :
1033 : // Now enfore merging scale cut if the event did not pass the merging scale
1034 : // criterion.
1035 0 : if ( nSteps > 0 && nSteps == nRequested && tmsnow < tmsval ) {
1036 0 : string message="Warning in Merging::cutOnProcess: Les Houches Event";
1037 0 : message+=" fails merging scale cut. Reject event.";
1038 0 : infoPtr->errorMsg(message);
1039 : return true;
1040 0 : }
1041 :
1042 0 : if ( FullHistory.select(RN)->nClusterings() < nSteps) {
1043 0 : string message="Warning in Merging::cutOnProcess: No clusterings";
1044 0 : message+=" found. History incomplete.";
1045 0 : infoPtr->errorMsg(message);
1046 0 : }
1047 :
1048 : // Done if no real-emission jets are present.
1049 0 : if ( !containsRealKin ) return false;
1050 :
1051 : // Now cut on events that contain an additional real-emission jet.
1052 : // Perform one reclustering for real emission kinematics, then apply merging
1053 : // scale cut on underlying Born kinematics.
1054 0 : if ( containsRealKin ) {
1055 0 : Event dummy = Event();
1056 : // Initialise temporary output of reclustering.
1057 0 : dummy.clear();
1058 0 : dummy.init( "(hard process-modified)", particleDataPtr );
1059 0 : dummy.clear();
1060 : // Recluster once.
1061 0 : FullHistory.getClusteredEvent( RN, nSteps, dummy );
1062 0 : double tnowNew = mergingHooksPtr->tmsNow( dummy );
1063 : // Veto if underlying Born kinematics do not pass merging scale cut.
1064 0 : if ( nSteps > 0 && nRequested > 0 && tnowNew < tmsval ) return true;
1065 0 : }
1066 :
1067 : // Done if only interested in cross section estimate after cuts.
1068 0 : return false;
1069 :
1070 0 : }
1071 :
1072 : //==========================================================================
1073 :
1074 : } // end namespace Pythia8
|
