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