StRoot  1
StringLength.cc
1 // StringLength.cc 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
7 // StringLength class.
8
9 // Calculate the lambda measure for string and junctions.
10
11 #include "Pythia8/StringLength.h"
12
13 namespace Pythia8 {
14
15 //==========================================================================
16
17 // The StringLength class.
18
19 //--------------------------------------------------------------------------
20
21 // Constants: could be changed here if desired, but normally should not.
22 // These are of technical nature, as described for each.
23
24 // Minimum energy of a parton and minimum angle between two partons.
25 // This is to avoid problems with infinities.
26 const double StringLength::TINY = 1e-20;
27 const double StringLength::MINANGLE = 1e-7;
28
29 //--------------------------------------------------------------------------
30
31 void StringLength::init(Info* infoPtrIn, Settings& settings) {
32
33  // Save pointers.
34  infoPtr = infoPtrIn;
35
36  // Store variables.
37  m0 = settings.parm("ColourReconnection:m0");
38  m0sqr = pow2(m0);
39  juncCorr = settings.parm("ColourReconnection:junctionCorrection");
40  sqrt2 = sqrt(2.);
41  lambdaForm = settings.mode("ColourReconnection:lambdaForm");
42 }
43
44 //--------------------------------------------------------------------------
45
46 // Calculate string length for two indices in the event record.
47
48 double StringLength::getStringLength( Event& event, int i, int j) {
49
50  // Find rest frame of particles.
51  Vec4 p1 = event[i].p();
52  Vec4 p2 = event[j].p();
53
54  return getStringLength(p1, p2);
55 }
56
57 //--------------------------------------------------------------------------
58
59 // Calculate string length for two particles given their four-momenta.
60
61 double StringLength::getStringLength( Vec4 p1, Vec4 p2) {
62
63  // Check for very small energies and angles.
64  if (p1.e() < TINY || p2.e() < TINY || theta(p1,p2) < MINANGLE) return 1e9;
65
66  // Boost to restframe.
67  Vec4 pSum = p1 + p2;
68  p1.bstback(pSum);
69  p2.bstback(pSum);
70
71  // Calculate string length.
72  Vec4 p0( 0., 0., 0., 1.);
73
74  return getLength(p1,p0) + getLength(p2,p0);
75 }
76
77 //--------------------------------------------------------------------------
78
79 // Calculate string length of a single particle.
80 // The first vector is the 4 vector of the particle.
81 // The second vector represents (1,0,0,0) in dipole restframe.
82
83 double StringLength::getLength(Vec4 p, Vec4 v, bool isJunc) {
84
85  double m = m0;
86  if (isJunc) m *= juncCorr;
87
88  if (lambdaForm == 0) return log (1. + sqrt2 * v * p/ m );
89  else if (lambdaForm == 1) return log (1. + 2 * v * p/ m );
90  else if (lambdaForm == 2) return log (2. * v * p / m);
91  else return 1e9;
92 }
93
94 //--------------------------------------------------------------------------
95
96 // Calculate the length of a single junction given the 3 entries in the event.
97
98 double StringLength::getJuncLength( Event& event, int i, int j, int k) {
99
100  if (i == j || i == k || j == k) return 1e9;
101
102  Vec4 p1 = event[i].p();
103  Vec4 p2 = event[j].p();
104  Vec4 p3 = event[k].p();
105
106  return getJuncLength(p1, p2, p3);
107 }
108 //--------------------------------------------------------------------------
109
110 // Calculate the length of a single junction given the 3 four-momenta.
111
112 double StringLength::getJuncLength(Vec4 p1, Vec4 p2, Vec4 p3) {
113
114  // Check for very small energies and angles.
115  if (p1.e() < TINY || p2.e() < TINY || p3.e() < TINY
116  || theta(p1,p2) < MINANGLE || theta(p1,p3) < MINANGLE
117  || theta(p2,p3) < MINANGLE) return 1e9;
118
119  // Find the junction rest frame.
120  RotBstMatrix MfromJRF1 = stringFragmentation.junctionRestFrame(p1,p2,p3);
121  MfromJRF1.invert();
122  Vec4 v1( 0., 0., 0., 1.);
123  v1.rotbst(MfromJRF1);
124
125  // Possible problem when the right system rest frame system is not found.
126  if ( pow2(p1*v1) - p1*p1 < 0. || pow2(p2*v1) - p2*p2 < 0.
127  || pow2(p3*v1) - p3*p3 < 0.) return 1e9;
128
129  // Calcualte the junction length.
130  return getLength(p1, v1, true) + getLength(p2, v1, true)
131  + getLength(p3, v1, true);
132 }
133
134 //--------------------------------------------------------------------------
135
136 // Calculate the length of a double junction given the 4 entries in the event.
137 // The first two are expected to be quarks, the second two to be anti quarks.
138
139 double StringLength::getJuncLength( Event& event, int i, int j, int k, int l) {
140
141  if (i == j || i == k || i == l || j == k || j == l || k == l) return 1e9;
142
143  // Simple minimum check of lengths.
144  double origLength = getStringLength(event, i, k)
145  + getStringLength(event, j, l);
146  double minLength = getStringLength(event, i, j)
147  + getStringLength(event, k, l);
148
149  if (origLength < minLength) return minLength;
150
151  Vec4 p1 = event[i].p();
152  Vec4 p2 = event[j].p();
153  Vec4 p3 = event[k].p();
154  Vec4 p4 = event[l].p();
155
156  return getJuncLength(p1, p2, p3, p4);
157 }
158
159 //--------------------------------------------------------------------------
160
161 // Calculate the length of a double junction given the 4 four-momenta.
162 // The first two are expected to be quarks, the second two to be anti quarks.
163
164 double StringLength::getJuncLength(Vec4 p1, Vec4 p2, Vec4 p3, Vec4 p4) {
165
166  // Check for very small energies and angles.
167  if ( p1.e() < TINY || p2.e() < TINY || p3.e() < TINY || p4.e() < TINY
168  || theta(p1,p2) < MINANGLE || theta(p1,p3) < MINANGLE
169  || theta(p1,p4) < MINANGLE || theta(p2,p3) < MINANGLE
170  || theta(p2,p4) < MINANGLE || theta(p3,p4) < MINANGLE) return 1e9;
171
172  // Calculate velocity of first junction.
173  Vec4 pSum1 = p3 +p4;
174
175  RotBstMatrix MfromJRF1 = stringFragmentation.junctionRestFrame(p1,p2,pSum1);
176  MfromJRF1.invert();
177  Vec4 v1( 0., 0., 0., 1.);
178  v1.rotbst(MfromJRF1);
179
180  // Calculate velocity of second junction.
181  Vec4 pSum2 = p1 + p2;
182
183  RotBstMatrix MfromJRF2 = stringFragmentation.junctionRestFrame(p3,p4,pSum2);
184  MfromJRF2.invert();
185  Vec4 v2( 0., 0., 0., 1.);
186  v2.rotbst(MfromJRF2);
187
188  // This only happens if it is not possible to find the correct rest frame.
189  if ( pow2(p1*v1) - p1*p1 < 0. || pow2(p2*v1) - p2*p2 < 0.
190  || pow2(p3*v2) - p3*p3 < 0. || pow2(p4*v2) - p4*p4 < 0.) return 1e9;
191
192  return getLength(p1, v1, true) + getLength(p2, v1, true)
193  + getLength(p3, v2, true) + getLength(p4, v2, true)
194  + log(v1*v2 + sqrt(pow2(v1*v2)-1));
195 }
196
197 //==========================================================================
198
199 } // end namespace Pythia8
Definition: AgUStep.h:26