StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ColourReconnectionHooks.h
1 // ColourReconnectionHooks.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 // This file contains two UserHooks that, along with the internal models,
7 // implement all the models used for the top mass study in
8 // S. Argyropoulos and T. Sjostrand,
9 // arXiv:1407.6653 [hep-ph] (LU TP 14-23, DESY 14-134, MCnet-14-15)
10 
11 // MBReconUserHooks: can be used for all kinds of events, not only top ones.
12 // TopReconUserHooks: models intended specifically for top decay products,
13 // whereas the underlying event is handled by the default model.
14 
15 // Warning: some small modifications have been made when collecting
16 // the models, but nothing intended to change the behaviour.
17 // Note: the move model is also available with ColourReconnection:mode = 2,
18 // while the ColourReconnection:mode = 1 model has not been used here.
19 // Note: the new models tend to be slower than the default CR scenario,
20 // since they have to probe many more reconnection possibilities.
21 
22 #ifndef Pythia8_ColourReconnectionHooks_H
23 #define Pythia8_ColourReconnectionHooks_H
24 
25 // Includes
26 #include "Pythia8/Pythia.h"
27 namespace Pythia8 {
28 
29 //==========================================================================
30 
31 // Class for colour reconnection models of general validity.
32 
33 class MBReconUserHooks : public UserHooks {
34 
35 public:
36 
37  // Constructor and destructor.
38  // mode = 0: no reconnection (dummy option, does nothing);
39  // = 1: swap gluons to minimize lambda.
40  // = 2: move gluons to minimize lambda.
41  // flip = 0: no flip between quark-antiquark ends.
42  // = 1: flip between quark-antiquark ends, excluding junction systems.
43  // = 2: flip between quark-antiquark ends, including junction systems.
44  // dLamCut: smallest -delta-lambda value for which to swap/mode (positive).
45  // fracGluon: the fraction of gluons that will be studied for reconnection.
46  // m2Ref : squared reference mass scale for lambda measure calculation.
47  MBReconUserHooks(int modeIn = 0, int flipIn = 0, double dLamCutIn = 0.,
48  double fracGluonIn = 1.) : mode(modeIn), flip(flipIn), dLamCut(dLamCutIn),
49  fracGluon(fracGluonIn) { m2Ref = 1.; dLamCut = max(0., dLamCut); }
50  ~MBReconUserHooks() {}
51 
52  // Allow colour reconnection after resonance decays (early or late)...
53  virtual bool canReconnectResonanceSystems() {return true;}
54 
55  // ...which gives access to the event, for modification.
56  virtual bool doReconnectResonanceSystems( int, Event& event) {
57 
58  // Return without action for relevant mode numbers.
59  if (mode <= 0 || mode > 2) return true;
60 
61  // Double diffraction not yet implemented, so return without action.
62  // (But works for internal move implementation.)
63  if (infoPtr->isDiffractiveA() && infoPtr->isDiffractiveB()) return true;
64 
65  // Initial setup: relevant gluons and coloured partons.
66  if (!setupConfig( event)) return false;
67 
68  // Done if not enough gluons.
69  if ( (mode == 1 && nGlu < 2) || (mode == 2 && nGlu < 1) ) return true;
70 
71  // Colour reconnect. Return if failed.
72  bool hasRec = (mode == 1) ? doReconnectSwap( event)
73  : doReconnectMove( event);
74  if (!hasRec) return false;
75 
76  // Colour flip afterburner.
77  if (flip > 0) return doReconnectFlip( event);
78  return true;
79 
80  }
81 
82  // Return number of reconnections for current event.
83  //int numberReconnections() {return nRec;}
84  //double dLambdaReconnections() {return -dLamTot;}
85 
86 private:
87 
88  // Mode. Number of reconnections. lambda measure reference scale.
89  int mode, flip, nRec, nGlu, nAllCol, nCol;
90  double dLamCut, fracGluon, m2Ref, dLamTot;
91 
92  // Array of (indices of) final gluons.
93  vector<int> iGlu;
94 
95  // Array of (indices of) all final coloured particles.
96  vector<int> iToAllCol, iAllCol;
97 
98  // Maps telling where all colours and anticolours are stored.
99  map<int, int> colMap, acolMap;
100 
101  // Array of all lambda distances between coloured partons.
102  vector<double> lambdaij;
103 
104  // Function to return lambda value from array.
105  double lambda12( int i, int j) {
106  int iAC = iToAllCol[i]; int jAC = iToAllCol[j];
107  return lambdaij[nAllCol * min( iAC, jAC) + max( iAC, jAC)];
108  }
109 
110  // Function to return lambda(i,j) + lambda(i,k) - lambda(j,k).
111  double lambda123( int i, int j, int k) {
112  int iAC = iToAllCol[i]; int jAC = iToAllCol[j]; int kAC = iToAllCol[k];
113  return lambdaij[nAllCol * min( iAC, jAC) + max( iAC, jAC)]
114  + lambdaij[nAllCol * min( iAC, kAC) + max( iAC, kAC)]
115  - lambdaij[nAllCol * min( jAC, kAC) + max( jAC, kAC)];
116  }
117 
118  // Small nested class for the effect of a potential gluon swap or move.
119  class InfoSwapMove{
120  public:
121  InfoSwapMove(int i1in = 0, int i2in = 0) : i1(i1in), i2(i2in) {}
122  InfoSwapMove(int i1in, int i2in, int iCol1in, int iAcol1in, int iCol2in,
123  int iAcol2in, int dLamIn) : i1(i1in), i2(i2in), iCol1(iCol1in),
124  iAcol1(iAcol1in), iCol2(iCol2in), iAcol2(iAcol2in), dLam(dLamIn) {}
125  ~InfoSwapMove() {}
126  int i1, i2, col1, acol1, iCol1, iAcol1, col2, acol2, iCol2, iAcol2;
127  double lamNow, dLam;
128  };
129 
130  // Vector of (1) gluon-pair swap or (2) gluon move properties.
131  vector<InfoSwapMove> infoSM;
132 
133  //----------------------------------------------------------------------
134 
135  // Initial setup: relevant gluons and coloured partons.
136 
137  inline bool setupConfig(Event& event) {
138 
139  // Reset arrays to prepare for the new event analysis.
140  iGlu.clear();
141  iToAllCol.resize( event.size() );
142  iAllCol.clear();
143  colMap.clear();
144  acolMap.clear();
145  infoSM.clear();
146  nRec = 0;
147  dLamTot = 0.;
148 
149  // Loop over all final particles. Store (fraction of) gluons.
150  for (int i = 3; i < event.size(); ++i) if (event[i].isFinal()) {
151  if (event[i].id() == 21 && rndmPtr->flat() < fracGluon)
152  iGlu.push_back(i);
153 
154  // Store location of all colour and anticolour particles and indices.
155  if (event[i].col() > 0 || event[i].acol() > 0) {
156  iToAllCol[i] = iAllCol.size();
157  iAllCol.push_back(i);
158  if (event[i].col() > 0) colMap[event[i].col()] = i;
159  if (event[i].acol() > 0) acolMap[event[i].acol()] = i;
160  }
161  }
162 
163  // Erase (anti)colours for (anti)junctions and skip adjacent gluons.
164  for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
165  if (event.kindJunction(iJun) == 1) {
166  for (int j = 0; j < 3; ++j) {
167  int jCol = event.colJunction( iJun, j);
168  for (map<int, int>::iterator colM = colMap.begin();
169  colM != colMap.end(); ++colM)
170  if (colM->first == jCol) {
171  colMap.erase( colM);
172  break;
173  }
174  for (int iG = 0; iG < int(iGlu.size()); ++iG)
175  if (event[iGlu[iG]].col() == jCol) {
176  iGlu.erase(iGlu.begin() + iG);
177  break;
178  }
179  }
180  } else if (event.kindJunction(iJun) == 2) {
181  for (int j = 0; j < 3; ++j) {
182  int jCol = event.colJunction( iJun, j);
183  for (map<int, int>::iterator acolM = acolMap.begin();
184  acolM != acolMap.end(); ++acolM)
185  if (acolM->first == jCol) {
186  acolMap.erase( acolM);
187  break;
188  }
189  for (int iG = 0; iG < int(iGlu.size()); ++iG)
190  if (event[iGlu[iG]].acol() == jCol) {
191  iGlu.erase(iGlu.begin() + iG);
192  break;
193  }
194  }
195  }
196  }
197 
198  // Error checks.
199  nGlu = iGlu.size();
200  nCol = colMap.size();
201  if (int(acolMap.size()) != nCol) {
202  infoPtr->errorMsg("Error in MBReconUserHooks: map sizes do not match");
203  return false;
204  }
205  map<int, int>::iterator colM = colMap.begin();
206  map<int, int>::iterator acolM = acolMap.begin();
207  for (int iCol = 0; iCol < nCol; ++iCol) {
208  if (colM->first != acolM->first) {
209  infoPtr->errorMsg("Error in MBReconUserHooks: map elements"
210  " do not match");
211  return false;
212  }
213  ++colM;
214  ++acolM;
215  }
216 
217  // Calculate and tabulate lambda between any pair of coloured partons.
218  nAllCol = iAllCol.size();
219  lambdaij.resize( pow2(nAllCol) );
220  int i, j;
221  for (int iAC = 0; iAC < nAllCol - 1; ++iAC) {
222  i = iAllCol[iAC];
223  for (int jAC = iAC + 1; jAC < nAllCol; ++jAC) {
224  j = iAllCol[jAC];
225  lambdaij[nAllCol * iAC + jAC]
226  = log(1. + m2( event[i], event[j]) / m2Ref);
227  }
228  }
229 
230  // Done.
231  return true;
232 
233  }
234 
235  //----------------------------------------------------------------------
236 
237  // Swap gluons by lambda measure.
238 
239  inline bool doReconnectSwap(Event& event) {
240 
241  // Set up initial possible gluon swap pairs with lambda gains/losses.
242  for (int iG1 = 0; iG1 < nGlu - 1; ++iG1) {
243  int i1 = iGlu[iG1];
244  for (int iG2 = iG1 + 1; iG2 < nGlu; ++iG2) {
245  int i2 = iGlu[iG2];
246  InfoSwapMove tmpSM( i1, i2);
247  calcLamSwap( tmpSM, event);
248  infoSM.push_back( tmpSM);
249  }
250  }
251  int nPair = infoSM.size();
252 
253  // Keep on looping over swaps until no further negative dLambda.
254  for ( int iSwap = 0; iSwap < nGlu; ++iSwap) {
255  int iPairMin = -1;
256  double dLamMin = 1e4;
257 
258  // Find lowest dLambda.
259  for (int iPair = 0; iPair < nPair; ++iPair)
260  if (infoSM[iPair].dLam < dLamMin) {
261  iPairMin = iPair;
262  dLamMin = infoSM[iPair].dLam;
263  }
264 
265  // Break if no shift below upper limit found.
266  if (dLamMin > -dLamCut) break;
267  ++nRec;
268  dLamTot += dLamMin;
269  int i1min = infoSM[iPairMin].i1;
270  int i2min = infoSM[iPairMin].i2;
271 
272  // Swap the colours in the event record.
273  int col1 = event[i1min].col();
274  int acol1 = event[i1min].acol();
275  int col2 = event[i2min].col();
276  int acol2 = event[i2min].acol();
277  event[i1min].cols( col2, acol2);
278  event[i2min].cols( col1, acol1);
279 
280  // Swap the indices in the colour maps.
281  colMap[col1] = i2min;
282  acolMap[acol1] = i2min;
283  colMap[col2] = i1min;
284  acolMap[acol2] = i1min;
285 
286  // Remove already swapped pair from further consideration.
287  infoSM[iPairMin] = infoSM.back();
288  infoSM.pop_back();
289  --nPair;
290 
291  // Update all pairs that have been affected.
292  for (int iPair = 0; iPair < nPair; ++iPair) {
293  InfoSwapMove& tmpSM = infoSM[iPair];
294  if ( tmpSM.i1 == i1min || tmpSM.i1 == i2min
295  || tmpSM.i2 == i1min || tmpSM.i2 == i2min
296  || tmpSM.iCol1 == i1min || tmpSM.iCol1 == i2min
297  || tmpSM.iAcol1 == i1min || tmpSM.iAcol1 == i2min
298  || tmpSM.iCol2 == i1min || tmpSM.iCol2 == i2min
299  || tmpSM.iAcol2 == i1min || tmpSM.iAcol2 == i2min)
300  calcLamSwap( tmpSM, event);
301  }
302  }
303 
304  // Done.
305  return true;
306 
307  }
308 
309  //----------------------------------------------------------------------
310 
311  // Calculate pair swap properties.
312 
313  inline void calcLamSwap( InfoSwapMove& tmpSM, Event& event) {
314 
315  // Colour line tracing to neighbours.
316  tmpSM.col1 = event[tmpSM.i1].col();
317  tmpSM.acol1 = event[tmpSM.i1].acol();
318  tmpSM.iCol1 = acolMap[tmpSM.col1];
319  tmpSM.iAcol1 = colMap[tmpSM.acol1];
320  tmpSM.col2 = event[tmpSM.i2].col();
321  tmpSM.acol2 = event[tmpSM.i2].acol();
322  tmpSM.iCol2 = acolMap[tmpSM.col2];
323  tmpSM.iAcol2 = colMap[tmpSM.acol2];
324 
325  // Lambda swap properties.
326  double lam1c = lambda12( tmpSM.i1, tmpSM.iCol1);
327  double lam1a = lambda12( tmpSM.i1, tmpSM.iAcol1);
328  double lam2c = lambda12( tmpSM.i2, tmpSM.iCol2);
329  double lam2a = lambda12( tmpSM.i2, tmpSM.iAcol2);
330  double lam3c = lambda12( tmpSM.i1, tmpSM.iCol2);
331  double lam3a = lambda12( tmpSM.i1, tmpSM.iAcol2);
332  double lam4c = lambda12( tmpSM.i2, tmpSM.iCol1);
333  double lam4a = lambda12( tmpSM.i2, tmpSM.iAcol1);
334  if (tmpSM.col1 == tmpSM.acol2 && tmpSM.acol1 == tmpSM.col2)
335  tmpSM.dLam = 2e4;
336  else if (tmpSM.col1 == tmpSM.acol2)
337  tmpSM.dLam = (lam3c + lam4a) - (lam1a + lam2c);
338  else if (tmpSM.acol1 == tmpSM.col2)
339  tmpSM.dLam = (lam3a + lam4c) - (lam1c + lam2a);
340  else tmpSM.dLam = (lam3c + lam3a + lam4c + lam4a)
341  - (lam1c + lam1a + lam2c + lam2a);
342 
343  // Done.
344  }
345 
346  //----------------------------------------------------------------------
347 
348  // Move gluons by lambda measure.
349 
350  inline bool doReconnectMove(Event& event) {
351 
352  // Temporary variables.
353  int iNow, colNow, acolNow, iColNow, iAcolNow, col2Now;
354  double lamNow;
355 
356  // Set up initial possible gluon moves with lambda gains/losses.
357  for (int iG = 0; iG < nGlu; ++iG) {
358 
359  // Gluon and its neighbours.
360  iNow = iGlu[iG];
361  colNow = event[iNow].col();
362  acolNow = event[iNow].acol();
363  iColNow = acolMap[colNow];
364  iAcolNow = colMap[acolNow];
365 
366  // Addition to Lambda of gluon in current position.
367  lamNow = lambda123( iNow, iColNow, iAcolNow);
368 
369  // Loop over all colour lines where gluon could be inserted.
370  for (map<int, int>::iterator colM = colMap.begin();
371  colM != colMap.end(); ++colM) {
372  col2Now = colM->first;
373 
374  // New container for gluon and colour line information.
375  InfoSwapMove tmpSM( iNow);
376  tmpSM.col1 = colNow;
377  tmpSM.acol1 = acolNow;
378  tmpSM.iCol1 = iColNow;
379  tmpSM.iAcol1 = iAcolNow;
380  tmpSM.lamNow = lamNow;
381  tmpSM.col2 = col2Now;
382  tmpSM.iCol2 = colMap[col2Now];
383  tmpSM.iAcol2 = acolMap[col2Now];
384 
385  // Addition to Lambda if gluon inserted on line.
386  tmpSM.dLam = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1
387  || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4
388  : lambda123( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2) - tmpSM.lamNow;
389  infoSM.push_back( tmpSM);
390  }
391  }
392  int nPair = infoSM.size();
393 
394  // Keep on looping over moves until no further negative dLambda.
395  for ( int iMove = 0; iMove < nGlu; ++iMove) {
396  int iPairMin = -1;
397  double dLamMin = 1e4;
398 
399  // Find lowest dLambda.
400  for (int iPair = 0; iPair < nPair; ++iPair)
401  if (infoSM[iPair].dLam < dLamMin) {
402  iPairMin = iPair;
403  dLamMin = infoSM[iPair].dLam;
404  }
405 
406  // Break if no shift below upper limit found.
407  if (dLamMin > -dLamCut) break;
408  ++nRec;
409  dLamTot += dLamMin;
410 
411  // Partons and colours involved in move.
412  InfoSwapMove& selSM = infoSM[iPairMin];
413  int i1Sel = selSM.i1;
414  int iCol1Sel = selSM.iCol1;
415  int iAcol1Sel = selSM.iAcol1;
416  int iCol2Sel = selSM.iCol2;
417  int iAcol2Sel = selSM.iAcol2;
418  int iCol2Mod[3] = { iAcol1Sel , i1Sel , iCol2Sel };
419  int col2Mod[3] = { selSM.col1, selSM.col2, selSM.acol1};
420 
421  // Remove gluon from old colour line and insert on new colour line.
422  for (int i = 0; i < 3; ++i) {
423  event[ iCol2Mod[i] ].col( col2Mod[i] );
424  colMap[ col2Mod[i] ] = iCol2Mod[i];
425  }
426 
427  // Update info for partons with new colors.
428  int i1Now = 0;
429  bool doUpdate = false;
430  for (int iPair = 0; iPair < nPair; ++iPair) {
431  InfoSwapMove& tmpSM = infoSM[iPair];
432  if (tmpSM.i1 != i1Now) {
433  i1Now = tmpSM.i1;
434  doUpdate = false;
435  if (i1Now == i1Sel || i1Now == iCol1Sel || i1Now == iAcol1Sel
436  || i1Now == iCol2Sel || i1Now == iAcol2Sel) {
437  colNow = event[i1Now].col();
438  acolNow = event[i1Now].acol();
439  iColNow = acolMap[colNow];
440  iAcolNow = colMap[acolNow];
441  lamNow = lambda123( i1Now, iColNow, iAcolNow);
442  doUpdate = true;
443  }
444  }
445  if (doUpdate) {
446  tmpSM.col1 = colNow;
447  tmpSM.acol1 = acolNow;
448  tmpSM.iCol1 = iColNow;
449  tmpSM.iAcol1 = iAcolNow;
450  tmpSM.lamNow = lamNow;
451  }
452  }
453 
454  // Update info on dLambda for affected particles and colour lines.
455  for (int iPair = 0; iPair < nPair; ++iPair) {
456  InfoSwapMove& tmpSM = infoSM[iPair];
457  int iMod = -1;
458  for (int i = 0; i < 3; ++i) if (tmpSM.col2 == col2Mod[i]) iMod = i;
459  if (iMod > -1) tmpSM.iCol2 = iCol2Mod[iMod];
460  if (tmpSM.i1 == i1Sel || tmpSM.i1 == iCol1Sel || tmpSM.i1 == iAcol1Sel
461  || tmpSM.i1 == iCol2Sel || tmpSM.i1 == iAcol2Sel || iMod > -1)
462  tmpSM.dLam = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1
463  || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4
464  : lambda123( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2) - tmpSM.lamNow;
465  }
466 
467  // End of loop over gluon shifting.
468  }
469 
470  // Done.
471  return true;
472 
473  }
474 
475  //----------------------------------------------------------------------
476 
477  // Flip colour chains by lambda measure.
478 
479  inline bool doReconnectFlip(Event& event) {
480 
481  // Array with colour lines, and where each line begins and ends.
482  vector<int> iTmp, iVec, iBeg, iEnd;
483 
484  // Grab all colour ends.
485  for (int i = 3; i < event.size(); ++i)
486  if (event[i].isFinal() && event[i].col() > 0 && event[i].acol() == 0) {
487  iTmp.clear();
488  iTmp.push_back( i);
489 
490  // Step through colour neighbours to catch system.
491  int iNow = i;
492  map<int, int>::iterator acolM = acolMap.find( event[iNow].col() );
493  bool foundEnd = false;
494  while (acolM != acolMap.end()) {
495  iNow = acolM->second;
496  iTmp.push_back( iNow);
497  if (event[iNow].col() == 0) {
498  foundEnd = true;
499  break;
500  }
501  acolM = acolMap.find( event[iNow].col() );
502  }
503 
504  // Store acceptable system, optionally including junction legs.
505  if (foundEnd || flip == 2) {
506  iBeg.push_back( iVec.size());
507  for (int j = 0; j < int(iTmp.size()); ++j) iVec.push_back( iTmp[j]);
508  iEnd.push_back( iVec.size());
509  }
510  }
511 
512  // Optionally search for antijunction legs: grab all anticolour ends.
513  if (flip == 2) for (int i = 3; i < event.size(); ++i)
514  if (event[i].isFinal() && event[i].acol() > 0 && event[i].col() == 0) {
515  iTmp.clear();
516  iTmp.push_back( i);
517 
518  // Step through anticolour neighbours to catch system.
519  int iNow = i;
520  map<int, int>::iterator colM = colMap.find( event[iNow].acol() );
521  bool foundEnd = false;
522  while (colM != colMap.end()) {
523  iNow = colM->second;
524  iTmp.push_back( iNow);
525  if (event[iNow].acol() == 0) {
526  foundEnd = true;
527  break;
528  }
529  colM = colMap.find( event[iNow].acol() );
530  }
531 
532  // Store acceptable system, but do not doublecount q - (n g) - qbar.
533  if (!foundEnd) {
534  iBeg.push_back( iVec.size());
535  for (int j = 0; j < int(iTmp.size()); ++j) iVec.push_back( iTmp[j]);
536  iEnd.push_back( iVec.size());
537  }
538  }
539 
540 
541  // Variables for minimum search.
542  int nSys = iBeg.size();
543  int i1c, i1a, i2c, i2a, i1cMin, i1aMin, i2cMin, i2aMin, iSMin;
544  double dLam, dLamMin;
545  vector<InfoSwapMove> flipMin;
546 
547  // Loop through all system pairs.
548  for (int iSys1 = 0; iSys1 < nSys - 1; ++iSys1) if (iBeg[iSys1] >= 0)
549  for (int iSys2 = iSys1 + 1; iSys2 < nSys; ++iSys2) if (iBeg[iSys2] >= 0) {
550  i1cMin = 0;
551  i1aMin = 0;
552  i2cMin = 0;
553  i2aMin = 0;
554  dLamMin = 1e4;
555 
556  // Loop through all possible flip locations for a pair.
557  for (int j1 = iBeg[iSys1]; j1 < iEnd[iSys1] - 1; ++j1)
558  for (int j2 = iBeg[iSys2]; j2 < iEnd[iSys2] - 1; ++j2) {
559  i1c = iVec[j1];
560  i1a = iVec[j1 + 1];
561  i2c = iVec[j2];
562  i2a = iVec[j2 + 1];
563  dLam = lambda12( i1c, i2a) + lambda12( i2c, i1a)
564  - lambda12( i1c, i1a) - lambda12( i2c, i2a);
565  if (dLam < dLamMin) {
566  i1cMin = i1c;
567  i1aMin = i1a;
568  i2cMin = i2c;
569  i2aMin = i2a;
570  dLamMin = dLam;
571  }
572  }
573 
574  // Store possible flips if low enough dLamMin.
575  if (dLamMin < -dLamCut) flipMin.push_back( InfoSwapMove(
576  iSys1, iSys2, i1cMin, i1aMin, i2cMin, i2aMin, dLamMin) );
577  }
578  int flipSize = flipMin.size();
579 
580  // Search for lowest possible flip among unused systems.
581  for (int iFlip = 0; iFlip < min( nSys / 2, flipSize); ++iFlip) {
582  iSMin = -1;
583  dLamMin = 1e4;
584  for (int iSys12 = 0; iSys12 < flipSize; ++iSys12)
585  if (flipMin[iSys12].i1 >= 0 && flipMin[iSys12].dLam < dLamMin) {
586  iSMin = iSys12;
587  dLamMin = flipMin[iSys12].dLam;
588  }
589 
590  // Do flip. Mark flipped systems.
591  if (iSMin >= 0) {
592  InfoSwapMove& flipNow = flipMin[iSMin];
593  int iS1 = flipNow.i1;
594  int iS2 = flipNow.i2;
595  event[ flipNow.iAcol1 ].acol( event[flipNow.iCol2].col() );
596  event[ flipNow.iAcol2 ].acol( event[flipNow.iCol1].col() );
597  ++nRec;
598  dLamTot += flipNow.dLam;
599  for (int iSys12 = 0; iSys12 < flipSize; ++iSys12)
600  if ( flipMin[iSys12].i1 == iS1 || flipMin[iSys12].i1 == iS2
601  || flipMin[iSys12].i2 == iS1 || flipMin[iSys12].i2 == iS2)
602  flipMin[iSys12].i1 = -1;
603  }
604  else break;
605  }
606 
607  // Done.
608  return true;
609 
610  }
611 
612 };
613 
614 //==========================================================================
615 
616 
617 // Class for colour reconnection models specifically aimed at top decays.
618 
619 class TopReconUserHooks : public UserHooks {
620 
621 public:
622 
623  // Constructor and destructor.
624  // mode = 0: no reconnection of tops (dummy option, does nothing);
625  // = 1: reconnect with random background gluon;
626  // = 2: reconnect with nearest (smallest-mass) background gluon;
627  // = 3: reconnect with furthest (largest-mass) background gluon;
628  // = 4: reconnect with smallest (with sign) lambda measure shift;
629  // = 5: reconnect only if reduced lamda, and then to most reduction.
630  // strength: fraction of top gluons that is to be colour reconnected.
631  // nList: list first nList parton classifications.
632  // pTolerance: acceptable total momentum error (in reconstruction check).
633  // m2Ref: squared reference mass scale for lambda measure calculation.
634  // Possible variants for the future: swap with nearest in angle, not mass,
635  // and/or only allow a background gluon to swap colours once.
636 
637  TopReconUserHooks(int modeIn = 0, double strengthIn = 1.) : mode(modeIn),
638  strength(strengthIn) { iList = 0; nList = 0; pTolerance = 0.01;
639  m2Ref = 1.;}
640  ~TopReconUserHooks() {}
641 
642  // Allow colour reconnection after resonance decays (early or late)...
643  virtual bool canReconnectResonanceSystems() {return true;}
644 
645  // ...which gives access to the event, for modification.
646  virtual bool doReconnectResonanceSystems( int, Event& event) {
647 
648  // Return without action for relevant mode numbers.
649  if (mode <= 0 || mode > 5) return true;
650 
651  // Classify coloured final partons.
652  classifyFinalPartons(event);
653 
654  // Check that classification worked as expected.
655  if (!checkClassification(event)) return false;
656 
657  // List first few classifications, along with the event.
658  if (iList++ < nList) {
659  listClassification();
660  event.list();
661  }
662 
663  // Perform reconnection for t and tbar in random order.
664  bool tqrkFirst = (rndmPtr->flat() < 0.5);
665  doReconnect( tqrkFirst, event);
666  doReconnect(!tqrkFirst, event);
667 
668  // Done.
669  return true;
670  }
671 
672  // Return number of reconnections for current event.
673  //int numberReconnections() {return nRec;}
674 
675 private:
676 
677  // Mode. Counters for how many events to list. Allowed momentum error.
678  int mode, iList, nList, nRec;
679  double strength, pTolerance, m2Ref;
680 
681  // Arrays of (indices of) final partons from different sources.
682  // So far geared towards t -> b W decays only.
683  vector<int> iBqrk, iWpos, iTqrk, iBbar, iWneg, iTbar, iRest;
684 
685  // Maps telling where all colours and anticolours are stored.
686  map<int, int> colMap, acolMap;
687 
688  //----------------------------------------------------------------------
689 
690  // Classify all coloured partons at the end of showers by origin.
691  // Note: for now only t -> b W is fully classified.
692 
693  inline bool classifyFinalPartons(Event& event) {
694 
695  // Reset arrays to prepare for the new event analysis.
696  iBqrk.clear();
697  iWpos.clear();
698  iTqrk.clear();
699  iBbar.clear();
700  iWneg.clear();
701  iTbar.clear();
702  iRest.clear();
703  colMap.clear();
704  acolMap.clear();
705  nRec = 0;
706 
707  // Loop over all final particles. Tag coloured ones.
708  for (int i = 3; i < event.size(); ++i) if (event[i].isFinal()) {
709  bool hasCol = (event[i].colType() != 0);
710 
711  // Set up to find where each parton comes from.
712  bool fsrFromT = false;
713  bool fromTqrk = false;
714  bool fromBqrk = false;
715  bool fromWpos = false;
716  bool fromTbar = false;
717  bool fromBbar = false;
718  bool fromWneg = false;
719 
720  // Identify current particle.
721  int iNow = i;
722  int idOld = 0;
723  do {
724  int idNow = event[iNow].id();
725 
726  // Exclude FSR of gluons/photons from the t quark proper.
727  if (abs(idNow) == 6 && (idOld == 21 || idOld == 22)) fsrFromT = true;
728 
729  // Check if current particle matches any of the categories.
730  else if (idNow == 6) fromTqrk = true;
731  else if (idNow == 5) fromBqrk = true;
732  else if (idNow == 24) fromWpos = true;
733  else if (idNow == -6) fromTbar = true;
734  else if (idNow == -5) fromBbar = true;
735  else if (idNow == -24) fromWneg = true;
736 
737  // Step up through the history to the very top.
738  iNow = event[iNow].mother1();
739  idOld = idNow;
740  } while (iNow > 2 && !fsrFromT);
741 
742  // Bookkeep where the parton comes from. Note that b quarks also
743  // can appear in W decays, so order of checks is relevant.
744  if (fromTqrk && fromWpos && hasCol) iWpos.push_back(i);
745  else if (fromTqrk && fromBqrk && hasCol) iBqrk.push_back(i);
746  else if (fromTqrk) iTqrk.push_back(i);
747  else if (fromTbar && fromWneg && hasCol) iWneg.push_back(i);
748  else if (fromTbar && fromBbar && hasCol) iBbar.push_back(i);
749  else if (fromTbar) iTbar.push_back(i);
750  else if (hasCol) iRest.push_back(i);
751 
752  // Store location of all colour and anticolour indices.
753  if (hasCol && (mode == 4 || mode == 5)) {
754  if (event[i].col() > 0) colMap[event[i].col()] = i;
755  if (event[i].acol() > 0) acolMap[event[i].acol()] = i;
756  }
757  }
758 
759  // So far method always returns true.
760  return true;
761 
762  }
763 
764  //----------------------------------------------------------------------
765 
766  // Check that classification worked by summing up partons to t/tbar.
767 
768  inline bool checkClassification(Event& event) {
769 
770  // Find final copy of t and tbar quarks.
771  int iTqrkLoc = 0;
772  int iTbarLoc = 0;
773  for (int i = 3; i < event.size(); ++i) {
774  if(event[i].id() == 6) iTqrkLoc = i;
775  if(event[i].id() == -6) iTbarLoc = i;
776  }
777 
778  // Four-momentum of t minus all its decay products.
779  Vec4 tqrkDiff = event[iTqrkLoc].p();
780  for (int i = 0; i < int(iBqrk.size()); ++i)
781  tqrkDiff -= event[iBqrk[i]].p();
782  for (int i = 0; i < int(iWpos.size()); ++i)
783  tqrkDiff -= event[iWpos[i]].p();
784  for (int i = 0; i < int(iTqrk.size()); ++i)
785  tqrkDiff -= event[iTqrk[i]].p();
786 
787  // Four-momentum of tbar minus all its decay products.
788  Vec4 tbarDiff = event[iTbarLoc].p();
789  for (int i = 0; i < int(iBbar.size()); ++i)
790  tbarDiff -= event[iBbar[i]].p();
791  for (int i = 0; i < int(iWneg.size()); ++i)
792  tbarDiff -= event[iWneg[i]].p();
793  for (int i = 0; i < int(iTbar.size()); ++i)
794  tbarDiff -= event[iTbar[i]].p();
795 
796  // Print difference vectors and event if sum deviation is too big.
797  double totErr = abs(tqrkDiff.px()) + abs(tqrkDiff.py())
798  + abs(tqrkDiff.pz()) + abs(tqrkDiff.e()) + abs(tbarDiff.px())
799  + abs(tbarDiff.py()) + abs(tbarDiff.pz()) + abs(tqrkDiff.e());
800  if (totErr > pTolerance) {
801  infoPtr->errorMsg("Error in TopReconUserHooks::checkClassification");
802  cout << "\n Error in t/tbar daughter search: \n t difference "
803  << tqrkDiff << " tbar difference "<< tbarDiff;
804  listClassification();
805  event.list();
806  }
807 
808  // Done.
809  return (totErr < pTolerance);
810  }
811 
812  //----------------------------------------------------------------------
813 
814  // Print how final-state (mainly coloured) particles were classified.
815 
816  inline void listClassification() {
817 
818  cout << "\n Final-state coloured partons classified by source: ";
819  cout << "\n From Bqrk:";
820  for (int i = 0; i < int(iBqrk.size()); ++i) cout << " " << iBqrk[i];
821  cout << "\n From Wpos:";
822  for (int i = 0; i < int(iWpos.size()); ++i) cout << " " << iWpos[i];
823  cout << "\n From Tqrk:";
824  for (int i = 0; i < int(iTqrk.size()); ++i) cout << " " << iTqrk[i];
825  cout << "\n From Bbar:";
826  for (int i = 0; i < int(iBbar.size()); ++i) cout << " " << iBbar[i];
827  cout << "\n From Wneg:";
828  for (int i = 0; i < int(iWneg.size()); ++i) cout << " " << iWneg[i];
829  cout << "\n From Tbar:";
830  for (int i = 0; i < int(iTbar.size()); ++i) cout << " " << iTbar[i];
831  cout << "\n From Rest:";
832  for (int i = 0; i < int(iRest.size()); ++i) {
833  cout << " " << iRest[i];
834  if (i%20 == 19 && i + 1 != int(iRest.size())) cout << "\n ";
835  }
836  cout << endl;
837  }
838 
839  //----------------------------------------------------------------------
840 
841  // Reconnect gluons either from t or from tbar quark.
842 
843  inline bool doReconnect(bool doTqrk, Event& event) {
844 
845  // Gather coloured decay products either of t or of tbar.
846  vector<int> iTdec;
847  if (doTqrk) {
848  for (int i = 0; i < int(iBqrk.size()); ++i) iTdec.push_back(iBqrk[i]);
849  for (int i = 0; i < int(iWpos.size()); ++i) iTdec.push_back(iWpos[i]);
850  } else {
851  for (int i = 0; i < int(iBbar.size()); ++i) iTdec.push_back(iBbar[i]);
852  for (int i = 0; i < int(iWneg.size()); ++i) iTdec.push_back(iWneg[i]);
853  }
854 
855  // Extract the gluons from the t quark decay.
856  vector<int> iGT;
857  for (int i = 0; i < int(iTdec.size()); ++i) {
858  int colNow = event[iTdec[i]].col();
859  int acolNow = event[iTdec[i]].acol();
860  if (colNow > 0 && acolNow > 0) iGT.push_back(iTdec[i]);
861  }
862  int nGT = iGT.size();
863 
864  // Randomize their stored order.
865  if (nGT > 1) for (int i = 0; i < nGT; ++i) {
866  int j = min( int(nGT * rndmPtr->flat()), nGT - 1 );
867  swap( iGT[i], iGT[j]);
868  }
869 
870  // Also extract the rest of the gluons in the event.
871  vector<int> iGR;
872  for (int i = 0; i < int(iRest.size()); ++i) {
873  int colNow = event[iRest[i]].col();
874  int acolNow = event[iRest[i]].acol();
875  if (colNow > 0 && acolNow > 0) iGR.push_back(iRest[i]);
876  }
877  int nGR = iGR.size();
878  int iR, colT, acolT, iColT, iAcolT, colR, acolR, iColR, iAcolR;
879  double mTR2, mTR2now, dLam, lamT, lamNow, lamRec;
880 
881  // Loop through all top gluons; study fraction given by strength.
882  if (nGT > 0 && nGR > 0)
883  for (int iT = 0; iT < nGT; ++iT) {
884  if (strength < rndmPtr->flat()) continue;
885 
886  // Pick random gluon from rest of event.
887  if (mode == 1) iR = min( int(nGR * rndmPtr->flat()), nGR - 1 );
888 
889  // Find gluon from rest with lowest or highest invariant mass.
890  else if (mode < 4) {
891  iR = 0;
892  mTR2 = m2( event[iGT[iT]], event[iGR[iR]]);
893  for (int ii = 1; ii < nGR; ++ii) {
894  mTR2now = m2( event[iGT[iT]], event[iGR[ii]]);
895  if (mode == 2 && mTR2now < mTR2) {iR = ii; mTR2 = mTR2now;}
896  if (mode == 3 && mTR2now > mTR2) {iR = ii; mTR2 = mTR2now;}
897  }
898 
899  // Find gluon from rest with smallest lambda value shift.
900  } else {
901  iR = -1;
902  dLam = 1e10;
903  colT = event[iGT[iT]].col();
904  acolT = event[iGT[iT]].acol();
905  iColT = acolMap[colT];
906  iAcolT = colMap[acolT];
907  lamT = log(1. + m2( event[iGT[iT]], event[iColT]) / m2Ref)
908  + log(1. + m2( event[iGT[iT]], event[iAcolT]) / m2Ref);
909  for (int ii = 0; ii < nGR; ++ii) {
910  colR = event[iGR[ii]].col();
911  acolR = event[iGR[ii]].acol();
912  iColR = acolMap[colR];
913  iAcolR = colMap[acolR];
914  lamNow = lamT
915  + log(1. + m2( event[iGR[ii]], event[iColR]) / m2Ref)
916  + log(1. + m2( event[iGR[ii]], event[iAcolR]) / m2Ref);
917  lamRec = log(1. + m2( event[iGT[iT]], event[iColR]) / m2Ref)
918  + log(1. + m2( event[iGT[iT]], event[iAcolR]) / m2Ref)
919  + log(1. + m2( event[iGR[ii]], event[iColT]) / m2Ref)
920  + log(1. + m2( event[iGR[ii]], event[iAcolT]) / m2Ref);
921  if (lamRec - lamNow < dLam) {iR = ii; dLam = lamRec - lamNow;}
922  }
923  }
924  if (mode == 5 && dLam > 0.) continue;
925 
926  // Swap top and rest gluon colour and anticolour.
927  ++nRec;
928  swapCols( iGT[iT], iGR[iR], event);
929  }
930 
931  // Done.
932  return true;
933 
934  }
935 
936  //----------------------------------------------------------------------
937 
938  // Swap colours and/or anticolours in the event listing.
939 
940  inline void swapCols( int i, int j, Event& event) {
941 
942  // Swap the colours in the event record.
943  int coli = event[i].col();
944  int acoli = event[i].acol();
945  int colj = event[j].col();
946  int acolj = event[j].acol();
947  event[i].cols( colj, acolj);
948  event[j].cols( coli, acoli);
949 
950  // Swap the indices in the colour maps.
951  colMap[coli] = j;
952  acolMap[acoli] = j;
953  colMap[colj] = i;
954  acolMap[acolj] = i;
955 
956  }
957 
958 };
959 
960 //==========================================================================
961 
962 } // end namespace Pythia8
963 
964 #endif // end Pythia8_ColourReconnectionHooks_H