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