StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
6 // Function definitions (not found in the header) for the
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