PYTHIA  8.311
ColourReconnection.h
1 // ColourReconnection.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 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 // Header file for the Colour reconnection handling.
7 // Reconnect the colours between the partons before hadronization.
8 // It Contains the following classes:
9 // ColourDipole, ColourParticle, ColourJunction, ColourReconnection.
10 
11 #ifndef Pythia8_ColourReconnection_H
12 #define Pythia8_ColourReconnection_H
13 
14 #include "Pythia8/Basics.h"
15 #include "Pythia8/BeamParticle.h"
16 #include "Pythia8/Event.h"
17 #include "Pythia8/FragmentationFlavZpT.h"
18 #include "Pythia8/Info.h"
19 #include "Pythia8/ParticleData.h"
20 #include "Pythia8/StringFragmentation.h"
21 #include "Pythia8/PartonDistributions.h"
22 #include "Pythia8/PartonSystems.h"
23 #include "Pythia8/PythiaStdlib.h"
24 #include "Pythia8/Settings.h"
25 #include "Pythia8/StringLength.h"
26 #include "Pythia8/StringInteractions.h"
27 
28 namespace Pythia8 {
29 
30 //==========================================================================
31 
32 // Contain a single colour chain. It always start from a quark and goes to
33 // an anti quark or from an anti-junction to at junction
34 // (or possible combinations).
35 
36 class ColourDipole {
37 
38 public:
39 
40  // Constructor.
41  ColourDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0,
42  int colReconnectionIn = 0, bool isJunIn = false, bool isAntiJunIn = false,
43  bool isActiveIn = true, bool isRealIn = false) : col(colIn), iCol(iColIn),
44  iAcol(iAcolIn), colReconnection(colReconnectionIn), isJun(isJunIn),
45  isAntiJun(isAntiJunIn),isActive(isActiveIn), isReal(isRealIn)
46  {iColLeg = 0; iAcolLeg = 0; printed = false; p1p2 = 0.;}
47 
48  double mDip(Event & event) {
49  if (isJun || isAntiJun) return 1E9;
50  else return m(event[iCol].p(),event[iAcol].p());
51  }
52 
53  // Members.
54  int col, iCol, iAcol, iColLeg, iAcolLeg, colReconnection;
55  bool isJun, isAntiJun, isActive, isReal, printed;
56  weak_ptr<ColourDipole> leftDip{}, rightDip{};
57  vector<weak_ptr<ColourDipole> > colDips, acolDips;
58  double p1p2;
59 
60  // Printing function, mainly intended for debugging.
61  void list();
62  long index{0};
63 
64 };
65 
66 // Comparison operator by index for two dipole pointers.
67 inline bool operator<(const ColourDipolePtr& d1, const ColourDipolePtr& d2) {
68  return ( d1 && d2 ? d1->index < d2->index : !d1 );}
69 
70 //==========================================================================
71 
72 // Junction class. In addition to the normal junction class, also contains a
73 // list of dipoles connected to it.
74 
75 class ColourJunction : public Junction {
76 
77 public:
78 
79  ColourJunction(const Junction& ju) : Junction(ju), dips(), dipsOrig() { }
80 
81  ColourJunction(const ColourJunction& ju) : Junction(Junction(ju)), dips(),
82  dipsOrig() { for(int i = 0;i < 3;++i) {
83  dips[i] = ju.dips[i]; dipsOrig[i] = ju.dipsOrig[i];}
84  }
85  ColourJunction& operator=( const ColourJunction& ju) {
86  this->Junction::operator=(ju);
87  for(int i = 0;i < 3;++i) {
88  dips[i] = ju.dips[i]; dipsOrig[i] = ju.dipsOrig[i];
89  }
90  return (*this);
91  }
92 
93  ColourDipolePtr getColDip(int i) {return dips[i];}
94  void setColDip(int i, ColourDipolePtr dip) {dips[i] = dip;}
95  ColourDipolePtr dips[3];
96  ColourDipolePtr dipsOrig[3];
97  void list();
98 
99 };
100 
101 //==========================================================================
102 
103 // TrialReconnection class.
104 
105 //--------------------------------------------------------------------------
106 
108 
109 public:
110 
111  TrialReconnection(ColourDipolePtr dip1In = 0, ColourDipolePtr dip2In = 0,
112  ColourDipolePtr dip3In = 0, ColourDipolePtr dip4In = 0, int modeIn = 0,
113  double lambdaDiffIn = 0) {
114  dips.push_back(dip1In); dips.push_back(dip2In);
115  dips.push_back(dip3In); dips.push_back(dip4In);
116  mode = modeIn; lambdaDiff = lambdaDiffIn;
117  }
118 
119  void list() {
120  cout << "mode: " << mode << " " << "lambdaDiff: " << lambdaDiff << endl;
121  for (int i = 0;i < int(dips.size()) && dips[i] != 0;++i) {
122  cout << " "; dips[i]->list(); }
123  }
124 
125  vector<ColourDipolePtr> dips;
126  int mode;
127  double lambdaDiff;
128 
129 };
130 
131 //==========================================================================
132 
133 // ColourParticle class.
134 
135 //--------------------------------------------------------------------------
136 
137 class ColourParticle : public Particle {
138 
139 public:
140 
141  ColourParticle(const Particle& ju) : Particle(ju), isJun(), junKind() {}
142 
143  vector<vector<ColourDipolePtr> > dips;
144  vector<bool> colEndIncluded, acolEndIncluded;
145  vector<ColourDipolePtr> activeDips;
146  bool isJun;
147  int junKind;
148 
149  // Printing functions, intended for debugging.
150  void listParticle();
151  void listActiveDips();
152  void listDips();
153 
154 };
155 
156 //==========================================================================
157 
158 // The ColourReconnection class handles the colour reconnection.
159 
160 //--------------------------------------------------------------------------
161 
163 
164 public:
165 
166  // Constructor
167  ColourReconnection() : allowJunctions(), sameNeighbourCol(),
168  singleReconOnly(), lowerLambdaOnly(), allowDiqJunCR(), nSys(),
169  nReconCols(), swap1(), swap2(), reconnectMode(), flipMode(),
170  timeDilationMode(), eCM(), sCM(), pT0(), pT20Rec(), pT0Ref(), ecmRef(),
171  ecmPow(), reconnectRange(), m0(), mPseudo(), m2Lambda(), fracGluon(),
172  dLambdaCut(), timeDilationPar(), timeDilationParGeV(), tfrag(), blowR(),
173  blowT(), rHadron(), kI(), dipMaxDist(), nColMove() {}
174 
175  // Initialization.
176  bool init();
177 
178  // New beams possible for handling of hard diffraction.
179  void reassignBeamPtrs( BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn)
180  {beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn;}
181 
182  // Do colour reconnection for current event.
183  bool next( Event & event, int oldSize);
184 
185 private:
186 
187  // Constants: could only be changed in the code itself.
188  static const double MINIMUMGAIN, MINIMUMGAINJUN, TINYP1P2;
189  static const int MAXRECONNECTIONS;
190 
191  // Variables needed.
192  bool allowJunctions, sameNeighbourCol, singleReconOnly, lowerLambdaOnly;
193  bool allowDiqJunCR;
194  int nSys, nReconCols, swap1, swap2, reconnectMode, flipMode,
195  timeDilationMode;
196  double eCM, sCM, pT0, pT20Rec, pT0Ref, ecmRef, ecmPow, reconnectRange,
197  m0, mPseudo, m2Lambda, fracGluon, dLambdaCut, timeDilationPar,
198  timeDilationParGeV, tfrag, blowR, blowT, rHadron, kI;
199  double dipMaxDist;
200 
201  // List of current dipoles.
202  vector<ColourDipolePtr> dipoles, usedDipoles;
203 
204  // Last used dipole index, used for sorting.
205  int dipoleIndex{0};
206 
207  // Add a dipole and increment index.
208  void addDipole(int colIn = 0, int iColIn = 0, int iAcolIn = 0,
209  int colReconnectionIn = 0, bool isJunIn = false, bool isAntiJunIn = false,
210  bool isActiveIn = true, bool isRealIn = false) {
211  dipoles.push_back(make_shared<ColourDipole>(colIn, iColIn, iAcolIn,
212  colReconnectionIn, isJunIn, isAntiJunIn, isActiveIn, isRealIn));
213  dipoles.back()->index = ++dipoleIndex;
214  }
215 
216  // Add a dipole and increment index.
217  void addDipole(const ColourDipole& dipole) {
218  dipoles.push_back(make_shared<ColourDipole>(dipole));
219  dipoles.back()->index = ++dipoleIndex;
220  }
221 
222  // Lists of particles, junctions and trials.
223  vector<ColourJunction> junctions;
224  vector<ColourParticle> particles;
225  vector<TrialReconnection> junTrials, dipTrials;
226  vector<vector<int> > iColJun;
227  map<int,double> formationTimes;
228 
229  // This is only to access the function call junctionRestFrame.
230  StringFragmentation stringFragmentation;
231 
232  // This class is used to calculate the string length.
233  StringLength stringLength;
234 
235  // Do colour reconnection for the event using the new model.
236  bool nextNew( Event & event, int oldSize);
237 
238  // Simple test swap between two dipoles.
239  void swapDipoles(ColourDipolePtr dip1, ColourDipolePtr dip2,
240  bool back = false);
241 
242  // Setup the dipoles.
243  void setupDipoles( Event& event, int iFirst = 0);
244 
245  // Form pseuparticle of a given dipole (or junction system).
246  void makePseudoParticle(ColourDipolePtr dip, int status,
247  bool setupDone = false);
248 
249  // Find the indices in the particle list of the junction and also their
250  // respectively leg numbers.
251  bool getJunctionIndices(ColourDipolePtr dip, int &iJun, int &i0, int &i1,
252  int &i2, int &junLeg0, int &junLeg1, int &junLeg2);
253 
254  // Form all possible pseudoparticles.
255  void makeAllPseudoParticles(Event & event, int iFirst = 0);
256 
257  // Update all colours in the event.
258  void updateEvent( Event& event, int iFirst = 0);
259 
260  double calculateStringLength( ColourDipolePtr dip,
261  vector<ColourDipolePtr> & dips);
262 
263  // Calculate the string length for two event indices.
264  double calculateStringLength( int i, int j);
265 
266  // Calculate the length of a single junction
267  // given the 3 entries in the particle list.
268  double calculateJunctionLength(int i, int j, int k);
269 
270  // Calculate the length of a double junction,
271  // given the 4 entries in the particle record.
272  // First two are expected to be the quarks and second two the anti quarks.
273  double calculateDoubleJunctionLength( int i, int j, int k, int l);
274 
275  // Find all the particles connected in the junction.
276  // If a single junction, the size of iParticles should be 3.
277  // For multiple junction structures, the size will increase.
278  bool findJunctionParticles( int iJun, vector<int>& iParticles,
279  vector<bool> &usedJuns, int &nJuns, vector<ColourDipolePtr> &dips);
280 
281  // Do a single trial reconnection.
282  void singleReconnection( ColourDipolePtr dip1, ColourDipolePtr dip2);
283 
284  // Do a single trial reconnection to form a junction.
285  void singleJunction(ColourDipolePtr dip1, ColourDipolePtr dip2);
286 
287  // Do a single trial reconnection to form a junction.
288  void singleJunction(ColourDipolePtr dip1, ColourDipolePtr dip2,
289  ColourDipolePtr dip3);
290 
291  // Print the chain containing the dipole.
292  void listChain(ColourDipolePtr dip);
293 
294  // Print all the chains.
295  void listAllChains();
296 
297  // Print dipoles, intended for debuggning purposes.
298  void listDipoles( bool onlyActive = false, bool onlyReal = false);
299 
300  // Print particles, intended for debugging purposes.
301  void listParticles();
302 
303  // Print junctions, intended for debugging purposes.
304  void listJunctions();
305 
306  // Check that the current dipole setup is consistent. Debug purpose only.
307  void checkDipoles();
308 
309  // Check that the current dipole setup is consistent. Debug purpose only.
310  void checkRealDipoles(Event& event, int iFirst);
311 
312  // Calculate the invariant mass of a dipole.
313  double mDip(ColourDipolePtr dip);
314 
315  // Find a vertex of the (anti)-colour side of the dipole. If
316  // connected to a junction, recurse using the other junction
317  // dipoles.
318  Vec4 getVProd(ColourDipolePtr dip, bool anti) const;
319 
320  // Find an average vertex of the (anti)-colour sides of the dipoles
321  // connected to the given junction (not incuding the given dipole).
322  Vec4 getVProd(int iJun, ColourDipolePtr dip, bool anti) const;
323 
324  // Check that the distance between the impact parameter centers of
325  // the dipoles are within the allowed range of dipMaxDist.
326  bool checkDist(ColourDipolePtr dip1, ColourDipolePtr dip2);
327 
328  // Find the neighbour to anti colour side. Return false if the dipole is
329  // connected to a junction or the new particle has a junction inside of it.
330  bool findAntiNeighbour(ColourDipolePtr& dip);
331 
332  // Find the neighbour to colour side. Return false if the dipole is
333  // connected to a junction or the new particle has a junction inside of it.
334  bool findColNeighbour(ColourDipolePtr& dip);
335 
336  // Check that trials do not contain junctions / unusable pseudoparticles.
337  bool checkJunctionTrials();
338 
339  // Store all dipoles connected to the ones used in the junction connection.
340  void storeUsedDips(TrialReconnection& trial);
341 
342  // Change colour structure to describe the reconnection in juncTrial.
343  bool doJunctionTrial(Event& event, TrialReconnection& juncTrial);
344 
345  // Change colour structure to describe the reconnection in juncTrial.
346  void doDipoleTrial(TrialReconnection& trial);
347 
348  // Change colour structure if it is three dipoles forming a junction system.
349  bool doTripleJunctionTrial(Event& event, TrialReconnection& juncTrial);
350 
351  // Calculate the difference between the old and new lambda.
352  double getLambdaDiff(ColourDipolePtr dip1, ColourDipolePtr dip2,
353  ColourDipolePtr dip3, ColourDipolePtr dip4, int mode);
354 
355  // Calculate the difference between the old and new lambda (dipole swap).
356  double getLambdaDiff(ColourDipolePtr dip1, ColourDipolePtr dip2);
357 
358  // Update the list of dipole trial swaps to account for latest swap.
359  void updateDipoleTrials();
360 
361  // Update the list of dipole trial swaps to account for latest swap.
362  void updateJunctionTrials();
363 
364  // Check whether up to four dipoles are 'causally' connected.
365  bool checkTimeDilation(ColourDipolePtr dip1 = 0, ColourDipolePtr dip2 = 0,
366  ColourDipolePtr dip3 = 0, ColourDipolePtr dip4 = 0);
367 
368  // Check whether two four momenta are casually connected.
369  bool checkTimeDilation(Vec4 p1, Vec4 p2, double t1, double t2);
370 
371  // Find the momentum of the dipole.
372  Vec4 getDipoleMomentum(ColourDipolePtr dip);
373 
374  // Find all particles connected to a junction system (particle list).
375  void addJunctionIndices(int iSinglePar, vector<int> &iPar,
376  vector<int> &usedJuncs);
377 
378  // Find all the formation times.
379  void setupFormationTimes( Event & event);
380 
381  // Get the mass of all partons connected to a junction system (event list).
382  double getJunctionMass(Event & event, int col);
383 
384  // Find all particles connected to a junction system (event list).
385  void addJunctionIndices(Event & event, int iSinglePar,
386  vector<int> &iPar, vector<int> &usedJuncs);
387 
388  // The old MPI-based scheme.
389  bool reconnectMPIs( Event& event, int oldSize);
390 
391  // Vectors and methods needed for the new gluon-move model.
392 
393  // Array of (indices of) all final coloured particles.
394  vector<int> iReduceCol, iExpandCol;
395 
396  // Array of all lambda distances between coloured partons.
397  int nColMove;
398  vector<double> lambdaijMove;
399 
400  // Function to return lambda value from array.
401  double lambda12Move( int i, int j) {
402  int iAC = iReduceCol[i]; int jAC = iReduceCol[j];
403  return lambdaijMove[nColMove * min( iAC, jAC) + max( iAC, jAC)];
404  }
405 
406  // Function to return lambda(i,j) + lambda(i,k) - lambda(j,k).
407  double lambda123Move( int i, int j, int k) {
408  int iAC = iReduceCol[i]; int jAC = iReduceCol[j]; int kAC = iReduceCol[k];
409  return lambdaijMove[nColMove * min( iAC, jAC) + max( iAC, jAC)]
410  + lambdaijMove[nColMove * min( iAC, kAC) + max( iAC, kAC)]
411  - lambdaijMove[nColMove * min( jAC, kAC) + max( jAC, kAC)];
412  }
413 
414  // The new gluon-move scheme.
415  bool reconnectMove( Event& event, int oldSize);
416 
417  // The common part for both Type I and II reconnections in e+e..
418  bool reconnectTypeCommon( Event& event, int oldSize);
419 
420  // The e+e- type I CR model.
421  map<double,pair<int,int> > reconnectTypeI( Event& event,
422  vector<vector<ColourDipole> > &dips, Vec4 decays[2]);
423  // bool reconnectTypeI( Event& event, int oldSize);
424 
425  // The e+e- type II CR model.
426  map<double,pair<int,int> > reconnectTypeII( Event& event,
427  vector<vector<ColourDipole> > &dips, Vec4 decays[2]);
428  // bool reconnectTypeII( Event& event, int oldSize);
429 
430  // calculate the determinant of 3 * 3 matrix.
431  double determinant3(vector<vector< double> >& vec);
432 };
433 
434 //==========================================================================
435 
436 } // end namespace Pythia8
437 
438 #endif // Pythia8_ColourReconnection_H
Definition: ColourReconnection.h:36
TrialReconnection class.
Definition: ColourReconnection.h:107
The Event class holds all info on the generated event.
Definition: Event.h:453
Definition: BeamParticle.h:133
void reassignBeamPtrs(BeamParticle *beamAPtrIn, BeamParticle *beamBPtrIn)
New beams possible for handling of hard diffraction.
Definition: ColourReconnection.h:179
Definition: Event.h:338
Definition: StringFragmentation.h:119
int col
Members.
Definition: ColourReconnection.h:54
Definition: StringInteractions.h:78
ColourDipole(int colIn=0, int iColIn=0, int iAcolIn=0, int colReconnectionIn=0, bool isJunIn=false, bool isAntiJunIn=false, bool isActiveIn=true, bool isRealIn=false)
Constructor.
Definition: ColourReconnection.h:41
StringLength class. It is used to calculate the lambda measure.
Definition: StringLength.h:23
ColourParticle class.
Definition: ColourReconnection.h:137
Definition: Event.h:32
Definition: ColourReconnection.h:75
The ColourReconnection class handles the colour reconnection.
Definition: ColourReconnection.h:162
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
double m(const Vec4 &v1)
Invariant mass and its square.
Definition: Basics.cc:576
void list()
Printing function, mainly intended for debugging.
Definition: ColourReconnection.cc:38
Definition: Basics.h:32
bool operator<(const ColourDipolePtr &d1, const ColourDipolePtr &d2)
Comparison operator by index for two dipole pointers.
Definition: ColourReconnection.h:67
ColourReconnection()
Constructor.
Definition: ColourReconnection.h:167