Line data Source code
1 : // MiniStringFragmentation.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 : // MiniStringFragmentation class
8 :
9 : #include "Pythia8/MiniStringFragmentation.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // The MiniStringFragmentation class.
16 :
17 : //--------------------------------------------------------------------------
18 :
19 : // Constants: could be changed here if desired, but normally should not.
20 : // These are of technical nature, as described for each.
21 :
22 : // Since diffractive by definition is > 1 particle, try hard.
23 : const int MiniStringFragmentation::NTRYDIFFRACTIVE = 200;
24 :
25 : // After one-body fragmentation failed, try two-body once more.
26 : const int MiniStringFragmentation::NTRYLASTRESORT = 100;
27 :
28 : // Loop try to combine available endquarks to valid hadron.
29 : const int MiniStringFragmentation::NTRYFLAV = 10;
30 :
31 : //--------------------------------------------------------------------------
32 :
33 : // Initialize and save pointers.
34 :
35 : void MiniStringFragmentation::init(Info* infoPtrIn, Settings& settings,
36 : ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
37 : StringFlav* flavSelPtrIn, StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {
38 :
39 : // Save pointers.
40 0 : infoPtr = infoPtrIn;
41 0 : particleDataPtr = particleDataPtrIn;
42 0 : rndmPtr = rndmPtrIn;
43 0 : flavSelPtr = flavSelPtrIn;
44 0 : pTSelPtr = pTSelPtrIn;
45 0 : zSelPtr = zSelPtrIn;
46 :
47 : // Initialize the MiniStringFragmentation class proper.
48 0 : nTryMass = settings.mode("MiniStringFragmentation:nTry");
49 :
50 : // Initialize the b parameter of the z spectrum, used when joining jets.
51 0 : bLund = zSelPtr->bAreaLund();
52 :
53 0 : }
54 :
55 : //--------------------------------------------------------------------------
56 :
57 : // Do the fragmentation: driver routine.
58 :
59 : bool MiniStringFragmentation::fragment(int iSub, ColConfig& colConfig,
60 : Event& event, bool isDiff) {
61 :
62 : // Junction topologies not yet considered - is very rare.
63 0 : iParton = colConfig[iSub].iParton;
64 0 : if (iParton.front() < 0) {
65 0 : infoPtr->errorMsg("Error in MiniStringFragmentation::fragment: "
66 : "very low-mass junction topologies not yet handled");
67 0 : return false;
68 : }
69 :
70 : // Read in info on system to be treated.
71 0 : flav1 = FlavContainer( event[ iParton.front() ].id() );
72 0 : flav2 = FlavContainer( event[ iParton.back() ].id() );
73 0 : pSum = colConfig[iSub].pSum;
74 0 : mSum = colConfig[iSub].mass;
75 0 : m2Sum = mSum*mSum;
76 0 : isClosed = colConfig[iSub].isClosed;
77 :
78 : // Do not want diffractive systems to easily collapse to one particle.
79 0 : int nTryFirst = (isDiff) ? NTRYDIFFRACTIVE : nTryMass;
80 :
81 : // First try to produce two particles from the system.
82 0 : if (ministring2two( nTryFirst, event)) return true;
83 :
84 : // If this fails, then form one hadron and shuffle momentum.
85 0 : if (ministring2one( iSub, colConfig, event)) return true;
86 :
87 : // If also this fails, then try harder to produce two particles.
88 0 : if (ministring2two( NTRYLASTRESORT, event)) return true;
89 :
90 : // Else complete failure.
91 0 : infoPtr->errorMsg("Error in MiniStringFragmentation::fragment: "
92 : "no 1- or 2-body state found above mass threshold");
93 0 : return false;
94 :
95 0 : }
96 :
97 : //--------------------------------------------------------------------------
98 :
99 : // Attempt to produce two particles from the ministring.
100 :
101 : bool MiniStringFragmentation::ministring2two( int nTry, Event& event) {
102 :
103 : // Properties of the produced hadrons.
104 : int idHad1 = 0;
105 : int idHad2 = 0;
106 : double mHad1 = 0.;
107 : double mHad2 = 0.;
108 : double mHadSum = 0.;
109 :
110 : // Allow a few attempts to find a particle pair with low enough masses.
111 0 : for (int iTry = 0; iTry < nTry; ++iTry) {
112 :
113 : // For closed gluon loop need to pick an initial flavour.
114 0 : if (isClosed) do {
115 0 : int idStart = flavSelPtr->pickLightQ();
116 0 : FlavContainer flavStart(idStart, 1);
117 0 : flavStart = flavSelPtr->pick( flavStart);
118 0 : flav1 = flavSelPtr->pick( flavStart);
119 0 : flav2.anti(flav1);
120 0 : } while (flav1.id == 0 || flav1.nPop > 0);
121 :
122 : // Create a new q qbar flavour to form two hadrons.
123 : // Start from a diquark, if any.
124 : do {
125 0 : FlavContainer flav3 =
126 0 : (flav1.isDiquark() || (!flav2.isDiquark() && rndmPtr->flat() < 0.5) )
127 0 : ? flavSelPtr->pick( flav1) : flavSelPtr->pick( flav2).anti();
128 0 : idHad1 = flavSelPtr->combine( flav1, flav3);
129 0 : idHad2 = flavSelPtr->combine( flav2, flav3.anti());
130 0 : } while (idHad1 == 0 || idHad2 == 0);
131 :
132 : // Check whether the mass sum fits inside the available phase space.
133 0 : mHad1 = particleDataPtr->mSel(idHad1);
134 0 : mHad2 = particleDataPtr->mSel(idHad2);
135 0 : mHadSum = mHad1 + mHad2;
136 0 : if (mHadSum < mSum) break;
137 : }
138 0 : if (mHadSum >= mSum) return false;
139 :
140 : // Define an effective two-parton string, by splitting intermediate
141 : // gluon momenta in proportion to their closeness to either endpoint.
142 0 : Vec4 pSum1 = event[ iParton.front() ].p();
143 0 : Vec4 pSum2 = event[ iParton.back() ].p();
144 0 : if (iParton.size() > 2) {
145 0 : Vec4 pEnd1 = pSum1;
146 0 : Vec4 pEnd2 = pSum2;
147 0 : Vec4 pEndSum = pEnd1 + pEnd2;
148 0 : for (int i = 1; i < int(iParton.size()) - 1 ; ++i) {
149 0 : Vec4 pNow = event[ iParton[i] ].p();
150 0 : double ratio = (pEnd2 * pNow) / (pEndSum * pNow);
151 0 : pSum1 += ratio * pNow;
152 0 : pSum2 += (1. - ratio) * pNow;
153 0 : }
154 0 : }
155 :
156 : // If split did not provide an axis then pick random axis to break tie.
157 : // (Needed for low-mass q-g-qbar with q-qbar perfectly parallel.)
158 0 : if (pSum1.mCalc() + pSum2.mCalc() > 0.999999 * mSum) {
159 0 : double cthe = 2. * rndmPtr->flat() - 1.;
160 0 : double sthe = sqrtpos(1. - cthe * cthe);
161 0 : double phi = 2. * M_PI * rndmPtr->flat();
162 0 : Vec4 delta = 0.5 * min( pSum1.e(), pSum2.e())
163 0 : * Vec4( sthe * sin(phi), sthe * cos(phi), cthe, 0.);
164 0 : pSum1 += delta;
165 0 : pSum2 -= delta;
166 0 : infoPtr->errorMsg("Warning in MiniStringFragmentation::ministring2two: "
167 : "random axis needed to break tie");
168 0 : }
169 :
170 : // Set up a string region based on the two effective endpoints.
171 0 : StringRegion region;
172 0 : region.setUp( pSum1, pSum2);
173 :
174 : // Generate an isotropic decay in the ministring rest frame,
175 : // suppressed at large pT by a fragmentation pT Gaussian.
176 0 : double pAbs2 = 0.25 * ( pow2(m2Sum - mHad1*mHad1 - mHad2*mHad2)
177 0 : - pow2(2. * mHad1 * mHad2) ) / m2Sum;
178 : double pT2 = 0.;
179 0 : do {
180 0 : double cosTheta = rndmPtr->flat();
181 0 : pT2 = (1. - pow2(cosTheta)) * pAbs2;
182 0 : } while (pTSelPtr->suppressPT2(pT2) < rndmPtr->flat() );
183 :
184 : // Construct the forward-backward asymmetry of the two particles.
185 0 : double mT21 = mHad1*mHad1 + pT2;
186 0 : double mT22 = mHad2*mHad2 + pT2;
187 0 : double lambda = sqrtpos( pow2(m2Sum - mT21 - mT22) - 4. * mT21 * mT22 );
188 0 : double probReverse = 1. / (1. + exp( min( 50., bLund * lambda) ) );
189 :
190 : // Construct kinematics, as viewed in the transverse rest frame.
191 0 : double xpz1 = 0.5 * lambda/ m2Sum;
192 0 : if (probReverse > rndmPtr->flat()) xpz1 = -xpz1;
193 0 : double xmDiff = (mT21 - mT22) / m2Sum;
194 0 : double xe1 = 0.5 * (1. + xmDiff);
195 0 : double xe2 = 0.5 * (1. - xmDiff );
196 :
197 : // Distribute pT isotropically in angle.
198 0 : double phi = 2. * M_PI * rndmPtr->flat();
199 0 : double pT = sqrt(pT2);
200 0 : double px = pT * cos(phi);
201 0 : double py = pT * sin(phi);
202 :
203 : // Translate this into kinematics in the string frame.
204 0 : Vec4 pHad1 = region.pHad( xe1 + xpz1, xe1 - xpz1, px, py);
205 0 : Vec4 pHad2 = region.pHad( xe2 - xpz1, xe2 + xpz1, -px, -py);
206 :
207 : // Mark hadrons from junction fragmentation with different status.
208 : int statusHadPos = 82, statusHadNeg = 82;
209 0 : if (abs(idHad1) > 1000 && abs(idHad1) < 10000 &&
210 0 : abs(idHad2) > 1000 && abs(idHad2) < 10000) {
211 0 : if (event[ iParton.front() ].statusAbs() == 74) statusHadPos = 89;
212 0 : if (event[ iParton.back() ].statusAbs() == 74) statusHadNeg = 89;
213 : }
214 0 : else if (abs(idHad1) > 1000 && abs(idHad1) < 10000) {
215 0 : if (event[ iParton.front() ].statusAbs() == 74 ||
216 0 : event[ iParton.back() ].statusAbs() == 74) statusHadPos = 89;
217 : }
218 0 : else if (abs(idHad2) > 1000 && abs(idHad2) < 10000) {
219 0 : if (event[ iParton.front() ].statusAbs() == 74 ||
220 0 : event[ iParton.back() ].statusAbs() == 74) statusHadNeg = 89;
221 : }
222 : // Add produced particles to the event record.
223 0 : int iFirst = event.append( idHad1, statusHadPos, iParton.front(),
224 0 : iParton.back(), 0, 0, 0, 0, pHad1, mHad1);
225 0 : int iLast = event.append( idHad2, statusHadNeg, iParton.front(),
226 0 : iParton.back(), 0, 0, 0, 0, pHad2, mHad2);
227 :
228 : // Set decay vertex when this is displaced.
229 0 : if (event[iParton.front()].hasVertex()) {
230 0 : Vec4 vDec = event[iParton.front()].vDec();
231 0 : event[iFirst].vProd( vDec );
232 0 : event[iLast].vProd( vDec );
233 0 : }
234 :
235 : // Set lifetime of hadrons.
236 0 : event[iFirst].tau( event[iFirst].tau0() * rndmPtr->exp() );
237 0 : event[iLast].tau( event[iLast].tau0() * rndmPtr->exp() );
238 :
239 : // Mark original partons as hadronized and set their daughter range.
240 0 : for (int i = 0; i < int(iParton.size()); ++i) {
241 0 : event[ iParton[i] ].statusNeg();
242 0 : event[ iParton[i] ].daughters(iFirst, iLast);
243 : }
244 :
245 : // Successfully done.
246 : return true;
247 :
248 0 : }
249 :
250 : //--------------------------------------------------------------------------
251 :
252 : // Attempt to produce one particle from a ministring.
253 : // Current algorithm: find the system with largest invariant mass
254 : // relative to the existing one, and boost that system appropriately.
255 : // Try more sophisticated alternatives later?? (Z0 mass shifted??)
256 : // Also, if problems, attempt several times to obtain closer mass match??
257 :
258 : bool MiniStringFragmentation::ministring2one( int iSub,
259 : ColConfig& colConfig, Event& event) {
260 :
261 : // Cannot handle qq + qbarqbar system.
262 0 : if (abs(flav1.id) > 100 && abs(flav2.id) > 100) return false;
263 :
264 : // For closed gluon loop need to pick an initial flavour.
265 0 : if (isClosed) do {
266 0 : int idStart = flavSelPtr->pickLightQ();
267 0 : FlavContainer flavStart(idStart, 1);
268 0 : flav1 = flavSelPtr->pick( flavStart);
269 0 : flav2 = flav1.anti();
270 0 : } while (abs(flav1.id) > 100);
271 :
272 : // Select hadron flavour from available quark flavours.
273 : int idHad = 0;
274 0 : for (int iTryFlav = 0; iTryFlav < NTRYFLAV; ++iTryFlav) {
275 0 : idHad = flavSelPtr->combine( flav1, flav2);
276 0 : if (idHad != 0) break;
277 : }
278 0 : if (idHad == 0) return false;
279 :
280 : // Find mass.
281 0 : double mHad = particleDataPtr->mSel(idHad);
282 :
283 : // Find the untreated parton system which combines to the largest
284 : // squared mass above mimimum required.
285 : int iMax = -1;
286 0 : double deltaM2 = mHad*mHad - mSum*mSum;
287 : double delta2Max = 0.;
288 0 : for (int iRec = iSub + 1; iRec < colConfig.size(); ++iRec) {
289 0 : double delta2Rec = 2. * (pSum * colConfig[iRec].pSum) - deltaM2
290 0 : - 2. * mHad * colConfig[iRec].mass;
291 0 : if (delta2Rec > delta2Max) { iMax = iRec; delta2Max = delta2Rec;}
292 : }
293 0 : if (iMax == -1) return false;
294 :
295 : // Construct kinematics of the hadron and recoiling system.
296 0 : Vec4& pRec = colConfig[iMax].pSum;
297 0 : double mRec = colConfig[iMax].mass;
298 0 : double vecProd = pSum * pRec;
299 0 : double coefOld = mSum*mSum + vecProd;
300 0 : double coefNew = mHad*mHad + vecProd;
301 0 : double coefRec = mRec*mRec + vecProd;
302 0 : double coefSum = coefOld + coefNew;
303 0 : double sHat = coefOld + coefRec;
304 0 : double root = sqrtpos( (pow2(coefSum) - 4. * sHat * mHad*mHad)
305 0 : / (pow2(vecProd) - pow2(mSum * mRec)) );
306 0 : double k2 = 0.5 * (coefOld * root - coefSum) / sHat;
307 0 : double k1 = (coefRec * k2 + 0.5 * deltaM2) / coefOld;
308 0 : Vec4 pHad = (1. + k1) * pSum - k2 * pRec;
309 0 : Vec4 pRecNew = (1. + k2) * pRec - k1 * pSum;
310 :
311 : // Mark hadrons from junction split off with status 89.
312 : int statusHad = 81;
313 0 : if (abs(idHad) > 1000 && abs(idHad) < 10000 &&
314 0 : (event[ iParton.front() ].statusAbs() == 74 ||
315 0 : event[ iParton.back() ].statusAbs() == 74)) statusHad = 89;
316 :
317 : // Add the produced particle to the event record.
318 0 : int iHad = event.append( idHad, statusHad, iParton.front(), iParton.back(),
319 0 : 0, 0, 0, 0, pHad, mHad);
320 :
321 : // Set decay vertex when this is displaced.
322 0 : if (event[iParton.front()].hasVertex()) {
323 0 : Vec4 vDec = event[iParton.front()].vDec();
324 0 : event[iHad].vProd( vDec );
325 0 : }
326 :
327 : // Set lifetime of hadron.
328 0 : event[iHad].tau( event[iHad].tau0() * rndmPtr->exp() );
329 :
330 : // Mark original partons as hadronized and set their daughter range.
331 0 : for (int i = 0; i < int(iParton.size()); ++i) {
332 0 : event[ iParton[i] ].statusNeg();
333 0 : event[ iParton[i] ].daughters(iHad, iHad);
334 : }
335 :
336 : // Copy down recoiling system, with boosted momentum. Update current partons.
337 0 : RotBstMatrix M;
338 0 : M.bst(pRec, pRecNew);
339 0 : for (int i = 0; i < colConfig[iMax].size(); ++i) {
340 0 : int iOld = colConfig[iMax].iParton[i];
341 : // Do not touch negative iOld = beginning of new junction leg.
342 0 : if (iOld >= 0) {
343 : int iNew;
344 : // Keep track of 74 throughout the event.
345 0 : if (event[iOld].status() == 74) iNew = event.copy(iOld, 74);
346 0 : else iNew = event.copy(iOld, 72);
347 0 : event[iNew].rotbst(M);
348 0 : colConfig[iMax].iParton[i] = iNew;
349 0 : }
350 : }
351 0 : colConfig[iMax].pSum = pRecNew;
352 0 : colConfig[iMax].isCollected = true;
353 :
354 : // Successfully done.
355 : return true;
356 :
357 0 : }
358 :
359 : //==========================================================================
360 :
361 : } // end namespace Pythia8
|
