1 | // RHadrons.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 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 RHadrons class. |
---|
7 | |
---|
8 | #include "RHadrons.h" |
---|
9 | |
---|
10 | namespace Pythia8 { |
---|
11 | |
---|
12 | //========================================================================== |
---|
13 | |
---|
14 | // The RHadrons class. |
---|
15 | |
---|
16 | //-------------------------------------------------------------------------- |
---|
17 | |
---|
18 | // Constants: could be changed here if desired, but normally should not. |
---|
19 | // These are of technical nature, as described for each. |
---|
20 | |
---|
21 | const int RHadrons::IDRHADSB[14] = { 1000512, 1000522, 1000532, |
---|
22 | 1000542, 1000552, 1005113, 1005211, 1005213, 1005223, 1005311, |
---|
23 | 1005313, 1005321, 1005323, 1005333 }; |
---|
24 | |
---|
25 | const int RHadrons::IDRHADST[14] = { 1000612, 1000622, 1000632, |
---|
26 | 1000642, 1000652, 1006113, 1006211, 1006213, 1006223, 1006311, |
---|
27 | 1006313, 1006321, 1006323, 1006333 }; |
---|
28 | |
---|
29 | // Gluino code and list of gluino R-hadron codes. |
---|
30 | const int RHadrons::IDRHADGO[38] = { 1000993, 1009113, 1009213, |
---|
31 | 1009223, 1009313, 1009323, 1009333, 1009413, 1009423, 1009433, |
---|
32 | 1009443, 1009513, 1009523, 1009533, 1009543, 1009553, 1091114, |
---|
33 | 1092114, 1092214, 1092224, 1093114, 1093214, 1093224, 1093314, |
---|
34 | 1093324, 1093334, 1094114, 1094214, 1094224, 1094314, 1094324, |
---|
35 | 1094334, 1095114, 1095214, 1095224, 1095314, 1095324, 1095334 }; |
---|
36 | |
---|
37 | // Allow a few tries for flavour selection since failure is allowed. |
---|
38 | const int RHadrons::NTRYMAX = 10; |
---|
39 | |
---|
40 | // Safety margin (in GeV) when constructing kinematics of system. |
---|
41 | const double RHadrons::MSAFETY = 0.1; |
---|
42 | |
---|
43 | // Maximal energy to borrow for gluon to insert on leg in to junction. |
---|
44 | const double RHadrons::EGBORROWMAX = 4.; |
---|
45 | |
---|
46 | //-------------------------------------------------------------------------- |
---|
47 | |
---|
48 | // Main routine to initialize R-hadron handling. |
---|
49 | |
---|
50 | bool RHadrons::init( Info* infoPtrIn, Settings& settings, |
---|
51 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn) { |
---|
52 | |
---|
53 | // Store input pointers for future use. |
---|
54 | infoPtr = infoPtrIn; |
---|
55 | particleDataPtr = particleDataPtrIn; |
---|
56 | rndmPtr = rndmPtrIn; |
---|
57 | |
---|
58 | // Flags and parameters related to R-hadron formation and decay. |
---|
59 | allowRH = settings.flag("RHadrons:allow"); |
---|
60 | maxWidthRH = settings.parm("RHadrons:maxWidth"); |
---|
61 | idRSb = settings.mode("RHadrons:idSbottom"); |
---|
62 | idRSt = settings.mode("RHadrons:idStop"); |
---|
63 | idRGo = settings.mode("RHadrons:idGluino"); |
---|
64 | setMassesRH = settings.flag("RHadrons:setMasses"); |
---|
65 | probGluinoballRH = settings.parm("RHadrons:probGluinoball"); |
---|
66 | mOffsetCloudRH = settings.parm("RHadrons:mOffsetCloud"); |
---|
67 | mCollapseRH = settings.parm("RHadrons:mCollapse"); |
---|
68 | diquarkSpin1RH = settings.parm("RHadrons:diquarkSpin1"); |
---|
69 | |
---|
70 | // Check whether sbottom, stop or gluino should form R-hadrons. |
---|
71 | allowRSb = allowRH && idRSb > 0 |
---|
72 | && (particleDataPtr->mWidth(idRSb) < maxWidthRH); |
---|
73 | allowRSt = allowRH && idRSt > 0 |
---|
74 | && (particleDataPtr->mWidth(idRSt) < maxWidthRH); |
---|
75 | allowRGo = allowRH && idRGo > 0 |
---|
76 | && (particleDataPtr->mWidth(idRGo) < maxWidthRH); |
---|
77 | allowSomeR = allowRSb || allowRSt || allowRGo; |
---|
78 | |
---|
79 | // Set nominal masses of sbottom R-mesons and R-baryons. |
---|
80 | if (allowRSb) { |
---|
81 | m0Sb = particleDataPtr->m0(idRSb); |
---|
82 | if (setMassesRH) { |
---|
83 | for (int i = 0; i < 14; ++i) { |
---|
84 | int idR = IDRHADSB[i]; |
---|
85 | double m0RHad = m0Sb + mOffsetCloudRH; |
---|
86 | m0RHad += particleDataPtr->constituentMass( (idR%100)/10); |
---|
87 | if (i > 4) |
---|
88 | m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 ); |
---|
89 | particleDataPtr->m0( idR, m0RHad); |
---|
90 | } |
---|
91 | } |
---|
92 | |
---|
93 | // Set widths and lifetimes of sbottom R-hadrons. |
---|
94 | double mWidthRHad = particleDataPtr->mWidth(idRSb); |
---|
95 | double tau0RHad = particleDataPtr->tau0( idRSb); |
---|
96 | for (int i = 0; i < 14; ++i) { |
---|
97 | particleDataPtr->mWidth( IDRHADSB[i], mWidthRHad); |
---|
98 | particleDataPtr->tau0( IDRHADSB[i], tau0RHad); |
---|
99 | } |
---|
100 | } |
---|
101 | |
---|
102 | // Set nominal masses of stop R-mesons and R-baryons. |
---|
103 | if (allowRSt) { |
---|
104 | m0St = particleDataPtr->m0(idRSt); |
---|
105 | if (setMassesRH) { |
---|
106 | for (int i = 0; i < 14; ++i) { |
---|
107 | int idR = IDRHADST[i]; |
---|
108 | double m0RHad = m0St + mOffsetCloudRH; |
---|
109 | m0RHad += particleDataPtr->constituentMass( (idR%100)/10); |
---|
110 | if (i > 4) |
---|
111 | m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 ); |
---|
112 | particleDataPtr->m0( idR, m0RHad); |
---|
113 | } |
---|
114 | } |
---|
115 | |
---|
116 | // Set widths and lifetimes of stop R-hadrons. |
---|
117 | double mWidthRHad = particleDataPtr->mWidth(idRSt); |
---|
118 | double tau0RHad = particleDataPtr->tau0( idRSt); |
---|
119 | for (int i = 0; i < 14; ++i) { |
---|
120 | particleDataPtr->mWidth( IDRHADST[i], mWidthRHad); |
---|
121 | particleDataPtr->tau0( IDRHADST[i], tau0RHad); |
---|
122 | } |
---|
123 | } |
---|
124 | |
---|
125 | // Set nominal masses of gluino R-glueballs, R-mesons and R-baryons. |
---|
126 | if (allowRGo) { |
---|
127 | m0Go = particleDataPtr->m0(idRGo); |
---|
128 | if (setMassesRH) { |
---|
129 | particleDataPtr->m0( IDRHADGO[0], m0Go + 2. * mOffsetCloudRH |
---|
130 | + particleDataPtr->constituentMass(21) ); |
---|
131 | for (int i = 1; i < 38; ++i) { |
---|
132 | int idR = IDRHADGO[i]; |
---|
133 | double m0RHad = m0Go + 2. * mOffsetCloudRH; |
---|
134 | m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 ); |
---|
135 | m0RHad += particleDataPtr->constituentMass( (idR%100)/10); |
---|
136 | if (i > 15) |
---|
137 | m0RHad += particleDataPtr->constituentMass( (idR%10000)/1000 ); |
---|
138 | particleDataPtr->m0( idR, m0RHad); |
---|
139 | } |
---|
140 | } |
---|
141 | |
---|
142 | // Set widths and lifetimes of gluino R-hadrons. |
---|
143 | double mWidthRHad = particleDataPtr->mWidth(idRGo); |
---|
144 | double tau0RHad = particleDataPtr->tau0( idRGo); |
---|
145 | for (int i = 0; i < 38; ++i) { |
---|
146 | particleDataPtr->mWidth( IDRHADGO[i], mWidthRHad); |
---|
147 | particleDataPtr->tau0( IDRHADGO[i], tau0RHad); |
---|
148 | } |
---|
149 | } |
---|
150 | |
---|
151 | // Done. |
---|
152 | return true; |
---|
153 | |
---|
154 | } |
---|
155 | //-------------------------------------------------------------------------- |
---|
156 | |
---|
157 | // Check if a given particle can produce R-hadrons. |
---|
158 | |
---|
159 | bool RHadrons::givesRHadron( int id) { |
---|
160 | if (allowRSb && abs(id) == idRSb) return true; |
---|
161 | if (allowRSt && abs(id) == idRSt) return true; |
---|
162 | if (allowRGo && id == idRGo) return true; |
---|
163 | return false; |
---|
164 | |
---|
165 | } |
---|
166 | |
---|
167 | //-------------------------------------------------------------------------- |
---|
168 | |
---|
169 | // Produce R-hadrons by fragmenting them off from existing strings. |
---|
170 | |
---|
171 | bool RHadrons::produce( ColConfig& colConfig, Event& event) { |
---|
172 | |
---|
173 | // Check whether some sparticles are allowed to hadronize. |
---|
174 | if (!allowSomeR) return true; |
---|
175 | |
---|
176 | // Reset arrays and values. |
---|
177 | iBefRHad.resize(0); |
---|
178 | iCreRHad.resize(0); |
---|
179 | iRHadron.resize(0); |
---|
180 | iAftRHad.resize(0); |
---|
181 | isTriplet.resize(0); |
---|
182 | nRHad = 0; |
---|
183 | |
---|
184 | // Loop over event and identify hadronizing sparticles. |
---|
185 | for (int i = 0; i < event.size(); ++i) |
---|
186 | if (event[i].isFinal() && givesRHadron(event[i].id())) { |
---|
187 | iBefRHad.push_back(i); |
---|
188 | iCreRHad.push_back(i); |
---|
189 | iRHadron.push_back(0); |
---|
190 | iAftRHad.push_back(0); |
---|
191 | isTriplet.push_back(true); |
---|
192 | } |
---|
193 | nRHad = iRHadron.size(); |
---|
194 | |
---|
195 | // Done if no hadronizing sparticles. |
---|
196 | if (nRHad == 0) return true; |
---|
197 | |
---|
198 | // Max two R-hadrons. Randomize order of processing. |
---|
199 | if (nRHad > 2) { |
---|
200 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
201 | "cannot handle more than two R-hadrons"); |
---|
202 | return false; |
---|
203 | } |
---|
204 | if (nRHad > 1 && rndmPtr->flat() > 0.5) swap( iBefRHad[0], iBefRHad[1]); |
---|
205 | |
---|
206 | // Split a system with both a sparticle and a junction. |
---|
207 | iBef = iBefRHad[0]; |
---|
208 | iSys = colConfig.findSinglet( iBef); |
---|
209 | systemPtr = &colConfig[iSys]; |
---|
210 | if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) { |
---|
211 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
212 | "cannot handle system with junction"); |
---|
213 | return false; |
---|
214 | } |
---|
215 | if (nRHad == 2) { |
---|
216 | iBef = iBefRHad[1]; |
---|
217 | iSys = colConfig.findSinglet( iBefRHad[1]); |
---|
218 | systemPtr = &colConfig[iSys]; |
---|
219 | if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) { |
---|
220 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
221 | "cannot handle system with junction"); |
---|
222 | return false; |
---|
223 | } |
---|
224 | } |
---|
225 | |
---|
226 | // Open up a closed gluon/gluino loop. |
---|
227 | iBef = iBefRHad[0]; |
---|
228 | iSys = colConfig.findSinglet( iBef); |
---|
229 | systemPtr = &colConfig[iSys]; |
---|
230 | if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) { |
---|
231 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
232 | "cannot open up closed gluon/gluino loop"); |
---|
233 | return false; |
---|
234 | } |
---|
235 | if (nRHad == 2) { |
---|
236 | iBef = iBefRHad[1]; |
---|
237 | iSys = colConfig.findSinglet( iBefRHad[1]); |
---|
238 | systemPtr = &colConfig[iSys]; |
---|
239 | if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) { |
---|
240 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
241 | "cannot open up closed gluon/gluino loop"); |
---|
242 | return false; |
---|
243 | } |
---|
244 | } |
---|
245 | |
---|
246 | // Split up a colour singlet system that should give two R-hadrons. |
---|
247 | if (nRHad == 2) { |
---|
248 | int iSys1 = colConfig.findSinglet( iBefRHad[0]); |
---|
249 | int iSys2 = colConfig.findSinglet( iBefRHad[1]); |
---|
250 | if (iSys2 == iSys1) { |
---|
251 | iSys = iSys1; |
---|
252 | systemPtr = &colConfig[iSys]; |
---|
253 | if ( !splitSystem( colConfig, event) ) { |
---|
254 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
255 | "failed to handle two sparticles in same system"); |
---|
256 | return false; |
---|
257 | } |
---|
258 | } |
---|
259 | } |
---|
260 | |
---|
261 | // Loop over R-hadrons to be formed. Find its colour singlet system. |
---|
262 | for (iRHad = 0; iRHad < nRHad; ++iRHad) { |
---|
263 | iBef = iBefRHad[iRHad]; |
---|
264 | iSys = colConfig.findSinglet( iBef); |
---|
265 | if (iSys < 0) { |
---|
266 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
267 | "sparticle not in any colour singlet"); |
---|
268 | return false; |
---|
269 | } |
---|
270 | systemPtr = &colConfig[iSys]; |
---|
271 | |
---|
272 | // For now don't handle systems involving junctions or loops. |
---|
273 | if (systemPtr->hasJunction) { |
---|
274 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
275 | "cannot handle system with junction"); |
---|
276 | return false; |
---|
277 | } |
---|
278 | if (systemPtr->isClosed) { |
---|
279 | infoPtr->errorMsg("Error in RHadrons::produce: " |
---|
280 | "cannot handle closed colour loop"); |
---|
281 | return false; |
---|
282 | } |
---|
283 | |
---|
284 | // Handle formation of R-hadron separately for gluino and squark. |
---|
285 | if (event[iBef].id() == idRGo) isTriplet[iRHad] = false; |
---|
286 | bool formed = (isTriplet[iRHad]) ? produceSquark( colConfig, event) |
---|
287 | : produceGluino( colConfig, event); |
---|
288 | if (!formed) return false; |
---|
289 | |
---|
290 | // End of loop over R-hadrons. Done. |
---|
291 | } |
---|
292 | return true; |
---|
293 | |
---|
294 | } |
---|
295 | |
---|
296 | //-------------------------------------------------------------------------- |
---|
297 | |
---|
298 | // Decay R-hadrons by resolving them into string systems and letting the |
---|
299 | // heavy unstable particle decay as normal. |
---|
300 | |
---|
301 | bool RHadrons::decay( Event& event) { |
---|
302 | |
---|
303 | // Loop over R-hadrons to decay. |
---|
304 | for (iRHad = 0; iRHad < nRHad; ++iRHad) { |
---|
305 | int iRNow = iRHadron[iRHad]; |
---|
306 | int iRBef = iBefRHad[iRHad]; |
---|
307 | int idRHad = event[iRNow].id(); |
---|
308 | double mRHad = event[iRNow].m(); |
---|
309 | double mRBef = event[iRBef].m(); |
---|
310 | int iR0 = 0; |
---|
311 | int iR2 = 0; |
---|
312 | |
---|
313 | // Find flavour content of squark or gluino R-hadron. |
---|
314 | pair<int,int> idPair = (isTriplet[iRHad]) |
---|
315 | ? fromIdWithSquark( idRHad) : fromIdWithGluino( idRHad); |
---|
316 | int id1 = idPair.first; |
---|
317 | int id2 = idPair.second; |
---|
318 | |
---|
319 | // Sharing of momentum: the squark/gluino should be restored |
---|
320 | // to original mass, but error if negative-mass spectators. |
---|
321 | double fracR = mRBef / mRHad; |
---|
322 | if (fracR >= 1.) { |
---|
323 | infoPtr->errorMsg("Error in RHadrons::decay: " |
---|
324 | "too low R-hadron mass for decay"); |
---|
325 | return false; |
---|
326 | } |
---|
327 | |
---|
328 | // Squark: new colour needed in the breakup. |
---|
329 | if (isTriplet[iRHad]) { |
---|
330 | int colNew = event.nextColTag(); |
---|
331 | int col = (event[iRBef].col() != 0) ? colNew : 0; |
---|
332 | int acol = (col == 0) ? colNew : 0; |
---|
333 | |
---|
334 | // Store the constituents of a squark R-hadron. |
---|
335 | iR0 = event.append( id1, 106, iRNow, 0, 0, 0, col, acol, |
---|
336 | fracR * event[iRNow].p(), fracR * mRHad, 0.); |
---|
337 | iR2 = event.append( id2, 106, iRNow, 0, 0, 0, acol, col, |
---|
338 | (1. - fracR) * event[iRNow].p(), (1. - fracR) * mRHad, 0.); |
---|
339 | |
---|
340 | // Gluino: set mass sharing between two spectators. |
---|
341 | } else { |
---|
342 | double m1Eff = particleDataPtr->constituentMass(id1) + mOffsetCloudRH; |
---|
343 | double m2Eff = particleDataPtr->constituentMass(id2) + mOffsetCloudRH; |
---|
344 | double frac1 = (1. - fracR) * m1Eff / ( m1Eff + m2Eff); |
---|
345 | double frac2 = (1. - fracR) * m2Eff / ( m1Eff + m2Eff); |
---|
346 | |
---|
347 | // Two new colours needed in the breakups. |
---|
348 | int col1 = event.nextColTag(); |
---|
349 | int col2 = event.nextColTag(); |
---|
350 | |
---|
351 | // Store the constituents of a gluino R-hadron. |
---|
352 | iR0 = event.append( idRGo, 106, iRNow, 0, 0, 0, col2, col1, |
---|
353 | fracR * event[iRNow].p(), fracR * mRHad, 0.); |
---|
354 | event.append( id1, 106, iRNow, 0, 0, 0, col1, 0, |
---|
355 | frac1 * event[iRNow].p(), frac1 * mRHad, 0.); |
---|
356 | iR2 = event.append( id2, 106, iRNow, 0, 0, 0, 0, col2, |
---|
357 | frac2 * event[iRNow].p(), frac2 * mRHad, 0.); |
---|
358 | } |
---|
359 | |
---|
360 | // Mark R-hadron as decayed and update history. |
---|
361 | event[iRNow].statusNeg(); |
---|
362 | event[iRNow].daughters( iR0, iR2); |
---|
363 | iAftRHad[iRHad] = iR0; |
---|
364 | |
---|
365 | // Set secondary vertex for decay products, but no lifetime. |
---|
366 | Vec4 vDec = event[iRNow].vProd() + event[iRNow].tau() |
---|
367 | * event[iR0].p() / event[iR0].m(); |
---|
368 | for (int iRd = iR0; iRd <= iR2; ++iRd) event[iRd].vProd( vDec); |
---|
369 | |
---|
370 | // End loop over R-hadron decays, based on velocity of squark. |
---|
371 | |
---|
372 | } |
---|
373 | |
---|
374 | // Done. |
---|
375 | return true; |
---|
376 | |
---|
377 | } |
---|
378 | |
---|
379 | //-------------------------------------------------------------------------- |
---|
380 | |
---|
381 | // Split a system that contains both a sparticle and a junction. |
---|
382 | |
---|
383 | bool RHadrons::splitOffJunction( ColConfig& colConfig, Event& event) { |
---|
384 | |
---|
385 | // Classify system into three legs, and find sparticle location. |
---|
386 | vector<int> leg1, leg2, leg3; |
---|
387 | int iLegSP = 0; |
---|
388 | int iPosSP = 0; |
---|
389 | int iLeg = 0; |
---|
390 | int iPos = 0; |
---|
391 | for (int i = 0; i < systemPtr->size(); ++i) { |
---|
392 | ++iPos; |
---|
393 | int iP = systemPtr->iParton[i]; |
---|
394 | if (iP < 0) { |
---|
395 | ++iLeg; |
---|
396 | iPos = -1; |
---|
397 | } else if (iLeg == 1) leg1.push_back( iP); |
---|
398 | else if (iLeg == 2) leg2.push_back( iP); |
---|
399 | else if (iLeg == 3) leg3.push_back( iP); |
---|
400 | if (iP == iBef) { |
---|
401 | iLegSP = iLeg; |
---|
402 | iPosSP = iPos; |
---|
403 | } |
---|
404 | } |
---|
405 | if (iLegSP == 0) return false; |
---|
406 | |
---|
407 | // Swap so leg 1 contains sparticle. If not innermost sparticle then |
---|
408 | // skip for now, and wait for this other (gluino!) to be split off. |
---|
409 | if (iLegSP == 2) swap(leg2, leg1); |
---|
410 | else if (iLegSP == 3) swap(leg3, leg1); |
---|
411 | for (int i = 0; i < iPosSP; ++i) |
---|
412 | if (event[leg1[i]].id() != 21) return true; |
---|
413 | |
---|
414 | // No gluon between sparticle and junction: find kinetic energy of system. |
---|
415 | if (iPosSP == 0) { |
---|
416 | Vec4 pSP = event[iBef].p(); |
---|
417 | Vec4 pRec = 0.; |
---|
418 | for (int i = 0; i < int(leg2.size()); ++i) pRec += event[leg2[i]].p(); |
---|
419 | for (int i = 0; i < int(leg3.size()); ++i) pRec += event[leg3[i]].p(); |
---|
420 | double mSys = (pSP + pRec).mCalc(); |
---|
421 | double mSP = pSP.mCalc(); |
---|
422 | double mRec = pRec.mCalc(); |
---|
423 | double eKin = mSys - mSP - mRec; |
---|
424 | |
---|
425 | // Insert new gluon that borrows part of kinetic energy. |
---|
426 | double mNewG = EGBORROWMAX * eKin / (EGBORROWMAX + eKin) ; |
---|
427 | Vec4 pNewG = (mNewG / mSys) * (pSP + pRec); |
---|
428 | int newCol = event.nextColTag(); |
---|
429 | bool isCol = (event[leg1.back()].col() > 0); |
---|
430 | int colNG = (isCol)? newCol : event[iBef].acol(); |
---|
431 | int acolNG = (isCol) ? event[iBef].col() : newCol; |
---|
432 | int iNewG = event.append( 21, 101, iBef, 0, 0, 0, colNG, acolNG, |
---|
433 | pNewG, mNewG, 0.); |
---|
434 | leg1.insert( leg1.begin(), iNewG); |
---|
435 | ++iPosSP; |
---|
436 | |
---|
437 | // Boosts for remainder systems that gave up energy. |
---|
438 | double mNewSys = mSys - mNewG; |
---|
439 | double pAbsOld = 0.5 * sqrtpos( pow2(mSys*mSys - mSP*mSP - mRec*mRec) |
---|
440 | - pow2(2. * mSP * mRec) ) / mSys; |
---|
441 | double pAbsNew = 0.5 * sqrtpos( pow2(mNewSys*mNewSys - mSP*mSP - mRec*mRec) |
---|
442 | - pow2(2. * mSP * mRec) ) / mNewSys; |
---|
443 | RotBstMatrix MtoCM, MfromCM, MSP, MRec; |
---|
444 | MtoCM.toCMframe( pSP, pRec); |
---|
445 | MfromCM = MtoCM; |
---|
446 | MfromCM.invert(); |
---|
447 | MSP = MtoCM; |
---|
448 | MSP.bst( 0., 0., -pAbsOld / sqrt(mSP * mSP + pAbsOld * pAbsOld)); |
---|
449 | MSP.bst( 0., 0., pAbsNew / sqrt(mSP * mSP + pAbsNew * pAbsNew)); |
---|
450 | MSP.rotbst( MfromCM); |
---|
451 | MRec = MtoCM; |
---|
452 | MRec.bst( 0., 0., pAbsOld / sqrt(mRec * mRec + pAbsOld * pAbsOld)); |
---|
453 | MRec.bst( 0., 0., -pAbsNew / sqrt(mRec * mRec + pAbsNew * pAbsNew)); |
---|
454 | MRec.rotbst( MfromCM); |
---|
455 | |
---|
456 | // Copy down recoling partons and boost their momenta. |
---|
457 | int iNewSP = event.copy( iBef, 101); |
---|
458 | event[iNewSP].mother2(0); |
---|
459 | event[iBef].daughter1(iNewG); |
---|
460 | event[iNewSP].rotbst( MSP); |
---|
461 | leg1[iPosSP] = iNewSP; |
---|
462 | if (iBefRHad[0] == iBef) iBefRHad[0] = iNewSP; |
---|
463 | else if (nRHad > 1 && iBefRHad[1] == iBef) iBefRHad[1] = iNewSP; |
---|
464 | iBef = iNewSP; |
---|
465 | for (int i = 0; i < int(leg2.size()); ++i) { |
---|
466 | int iCopy = event.copy( leg2[i], 101); |
---|
467 | event[iCopy].rotbst( MRec); |
---|
468 | if (iBefRHad[0] == leg2[i]) iBefRHad[0] = iCopy; |
---|
469 | else if (nRHad > 1 && iBefRHad[1] == leg2[i]) iBefRHad[1] = iCopy; |
---|
470 | leg2[i] = iCopy; |
---|
471 | } |
---|
472 | for (int i = 0; i < int(leg3.size()); ++i) { |
---|
473 | int iCopy = event.copy( leg3[i], 101); |
---|
474 | event[iCopy].rotbst( MRec); |
---|
475 | if (iBefRHad[0] == leg3[i]) iBefRHad[0] = iCopy; |
---|
476 | else if (nRHad > 1 && iBefRHad[1] == leg3[i]) iBefRHad[1] = iCopy; |
---|
477 | leg3[i] = iCopy; |
---|
478 | } |
---|
479 | |
---|
480 | // Now always at least one gluon between sparticle and junction. |
---|
481 | } |
---|
482 | |
---|
483 | // Find gluon with largest energy in sparticle rest frame. |
---|
484 | int iGspl = 0; |
---|
485 | double eGspl = event[leg1[0]].p() * event[iBef].p(); |
---|
486 | for (int i = 1; i < iPosSP; ++i) { |
---|
487 | double eGnow = event[leg1[i]].p() * event[iBef].p(); |
---|
488 | if (eGnow > eGspl) { |
---|
489 | iGspl = i; |
---|
490 | eGspl = eGnow; |
---|
491 | } |
---|
492 | } |
---|
493 | int iG = leg1[iGspl]; |
---|
494 | |
---|
495 | // Split this gluon into a collinear quark.antiquark pair. |
---|
496 | int idNewQ = flavSelPtr->pickLightQ(); |
---|
497 | int iNewQ = event.append( idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0, |
---|
498 | 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.); |
---|
499 | int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(), |
---|
500 | 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.); |
---|
501 | event[iG].statusNeg(); |
---|
502 | event[iG].daughters( iNewQ, iNewQb); |
---|
503 | if (event[leg1.back()].col() == 0) swap( iNewQ, iNewQb); |
---|
504 | |
---|
505 | // Collect two new systems after split. |
---|
506 | vector<int> iNewSys1, iNewSys2; |
---|
507 | iNewSys1.push_back( iNewQb); |
---|
508 | for (int i = iGspl + 1; i < int(leg1.size()); ++i) |
---|
509 | iNewSys1.push_back( leg1[i]); |
---|
510 | iNewSys2.push_back( -10); |
---|
511 | for (int i = 0; i < iGspl; ++i) iNewSys2.push_back( leg1[i]); |
---|
512 | iNewSys2.push_back( iNewQ); |
---|
513 | iNewSys2.push_back( -11); |
---|
514 | for (int i = 0; i < int(leg2.size()); ++i) iNewSys2.push_back( leg2[i]); |
---|
515 | iNewSys2.push_back( -12); |
---|
516 | for (int i = 0; i < int(leg3.size()); ++i) iNewSys2.push_back( leg3[i]); |
---|
517 | |
---|
518 | // Remove old system and insert two new ones. |
---|
519 | colConfig.erase(iSys); |
---|
520 | colConfig.insert( iNewSys1, event); |
---|
521 | colConfig.insert( iNewSys2, event); |
---|
522 | |
---|
523 | // Done. |
---|
524 | return true; |
---|
525 | |
---|
526 | } |
---|
527 | |
---|
528 | //-------------------------------------------------------------------------- |
---|
529 | |
---|
530 | // Open up a closed gluon/gluino loop. |
---|
531 | |
---|
532 | bool RHadrons::openClosedLoop( ColConfig& colConfig, Event& event) { |
---|
533 | |
---|
534 | // Find gluon with largest energy in gluino rest frame. |
---|
535 | int iGspl = -1; |
---|
536 | double eGspl = 0.; |
---|
537 | for (int i = 0; i < systemPtr->size(); ++i) { |
---|
538 | int iGNow = systemPtr->iParton[i]; |
---|
539 | if (event[iGNow].id() == 21) { |
---|
540 | double eGnow = event[iGNow].p() * event[iBef].p(); |
---|
541 | if (eGnow > eGspl) { |
---|
542 | iGspl = i; |
---|
543 | eGspl = eGnow; |
---|
544 | } |
---|
545 | } |
---|
546 | } |
---|
547 | if (iGspl == -1) return false; |
---|
548 | |
---|
549 | // Split this gluon into a collinear quark.antiquark pair. |
---|
550 | int iG = systemPtr->iParton[iGspl]; |
---|
551 | int idNewQ = flavSelPtr->pickLightQ(); |
---|
552 | int iNewQ = event.append( idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0, |
---|
553 | 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.); |
---|
554 | int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(), |
---|
555 | 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.); |
---|
556 | event[iG].statusNeg(); |
---|
557 | event[iG].daughters( iNewQ, iNewQb); |
---|
558 | |
---|
559 | // Order partons in new system. Note order of colour flow. |
---|
560 | int iNext = iGspl + 1; |
---|
561 | if (iNext == systemPtr->size()) iNext = 0; |
---|
562 | if (event[ systemPtr->iParton[iNext]].acol() != event[iNewQ].col()) |
---|
563 | swap( iNewQ, iNewQb); |
---|
564 | vector<int> iNewSys; |
---|
565 | iNewSys.push_back( iNewQ); |
---|
566 | for (int i = iGspl + 1; i < systemPtr->size(); ++i) |
---|
567 | iNewSys.push_back( systemPtr->iParton[i]); |
---|
568 | for (int i = 0; i < iGspl; ++i) |
---|
569 | iNewSys.push_back( systemPtr->iParton[i]); |
---|
570 | iNewSys.push_back( iNewQb); |
---|
571 | |
---|
572 | // Erase the old system and insert the new one instead. |
---|
573 | colConfig.erase(iSys); |
---|
574 | colConfig.insert( iNewSys, event); |
---|
575 | |
---|
576 | // Done. |
---|
577 | return true; |
---|
578 | |
---|
579 | } |
---|
580 | |
---|
581 | //-------------------------------------------------------------------------- |
---|
582 | |
---|
583 | // Split a single colour singlet that contains two sparticles. |
---|
584 | // To fix : if nLeg >= 3 && mMin large handle as in nLeg == 1?? |
---|
585 | |
---|
586 | bool RHadrons::splitSystem( ColConfig& colConfig, Event& event) { |
---|
587 | |
---|
588 | // First and second R-hadron mother. Number of legs in between. |
---|
589 | int iFirst = -1; |
---|
590 | int iSecond = -1; |
---|
591 | for (int i = 0; i < int(systemPtr->size()); ++i) { |
---|
592 | int iTmp = systemPtr->iParton[i]; |
---|
593 | if ( givesRHadron( event[iTmp].id()) ) { |
---|
594 | if (iFirst == -1) iFirst = i; |
---|
595 | else iSecond = i; |
---|
596 | } |
---|
597 | } |
---|
598 | int nLeg = iSecond - iFirst; |
---|
599 | |
---|
600 | // New flavour pair for breaking the string, and its mass. |
---|
601 | int idNewQ = flavSelPtr->pickLightQ(); |
---|
602 | double mNewQ = particleDataPtr->constituentMass( idNewQ); |
---|
603 | vector<int> iNewSys1, iNewSys2; |
---|
604 | |
---|
605 | // If sparticles are neighbours then need new q-qbar pair in between. |
---|
606 | if (nLeg == 1) { |
---|
607 | |
---|
608 | // The location of the two sparticles. |
---|
609 | int i1Old = systemPtr->iParton[iFirst]; |
---|
610 | int i2Old = systemPtr->iParton[iSecond]; |
---|
611 | |
---|
612 | // Take a fraction of momentum to give breakup quark pair. |
---|
613 | double mHat = (event[i1Old].p() + event[i2Old].p()).mCalc(); |
---|
614 | double mMax = mHat - event[i1Old].m() - event[i2Old].m(); |
---|
615 | if (mMax < 2. * (mNewQ + MSAFETY)) return false; |
---|
616 | double mEff = min( 2. * (mNewQ + mOffsetCloudRH), mMax - 2. * MSAFETY); |
---|
617 | double frac = mEff / mHat; |
---|
618 | Vec4 pEff = frac * (event[i1Old].p() + event[i2Old].p()); |
---|
619 | |
---|
620 | // New kinematics by (1) go to same mHat with bigger masses, and |
---|
621 | // (2) rescale back down to original masses and reduced mHat. |
---|
622 | Vec4 p1New, p2New; |
---|
623 | if ( !newKin( event[i1Old].p(), event[i2Old].p(), |
---|
624 | event[i1Old].m() / (1. - frac), event[i2Old].m() / (1. - frac), |
---|
625 | p1New, p2New) ) return false; |
---|
626 | p1New *= 1. - frac; |
---|
627 | p2New *= 1. - frac; |
---|
628 | |
---|
629 | // Fill in new partons after branching, with correct colour/flavour sign. |
---|
630 | int i1New, i2New, i3New, i4New; |
---|
631 | int newCol = event.nextColTag(); |
---|
632 | i1New = event.copy( i1Old, 101); |
---|
633 | if (event[i2Old].acol() == event[i1Old].col()) { |
---|
634 | i3New = event.append( -idNewQ, 101, i1Old, 0, 0, 0, |
---|
635 | 0, event[i2Old].acol(), 0.5 * pEff, 0.5 * mEff, 0.); |
---|
636 | i2New = event.copy( i2Old, 101); |
---|
637 | event[i2New].acol( newCol); |
---|
638 | i4New = event.append( idNewQ, 101, i2Old, 0, 0, 0, |
---|
639 | newCol, 0, 0.5 * pEff, 0.5 * mEff, 0.); |
---|
640 | } else { |
---|
641 | i3New = event.append( idNewQ, 101, i1Old, 0, 0, 0, |
---|
642 | event[i2Old].col(), 0, 0.5 * pEff, 0.5 * mEff, 0.); |
---|
643 | i2New = event.copy( i2Old, 101); |
---|
644 | event[i2New].col( newCol); |
---|
645 | i4New = event.append( -idNewQ, 101, i2Old, 0, 0, 0, |
---|
646 | 0, newCol, 0.5 * pEff, 0.5 * mEff, 0.); |
---|
647 | } |
---|
648 | |
---|
649 | // Modify momenta and history. For iBotCopyId tracing asymmetric |
---|
650 | // bookkeeping where one sparticle is radiator and the other recoiler. |
---|
651 | event[i1New].p( p1New); |
---|
652 | event[i2New].p( p2New); |
---|
653 | event[i1Old].daughters( i1New, i3New); |
---|
654 | event[i1New].mother2( 0); |
---|
655 | event[i2Old].daughters( i2New, i4New); |
---|
656 | event[i2New].mother2( 0); |
---|
657 | iBefRHad[0] = i1New; |
---|
658 | iBefRHad[1] = i2New; |
---|
659 | |
---|
660 | // Partons in the two new systems. |
---|
661 | for (int i = 0; i < iFirst; ++i) |
---|
662 | iNewSys1.push_back( systemPtr->iParton[i]); |
---|
663 | iNewSys1.push_back( i1New); |
---|
664 | iNewSys1.push_back( i3New); |
---|
665 | iNewSys2.push_back( i4New); |
---|
666 | iNewSys2.push_back( i2New); |
---|
667 | for (int i = iSecond + 1; i < int(systemPtr->size()); ++i) |
---|
668 | iNewSys2.push_back( systemPtr->iParton[i]); |
---|
669 | |
---|
670 | // If one gluon between sparticles then split it into a |
---|
671 | // collinear q-qbar pair. |
---|
672 | } else if (nLeg == 2) { |
---|
673 | |
---|
674 | // Gluon to be split and its two daughters. |
---|
675 | int iOld = systemPtr->iParton[iFirst + 1]; |
---|
676 | int i1New = event.append( idNewQ, 101, iOld, 0, 0, 0, |
---|
677 | event[iOld].col(), 0, 0.5 * event[iOld].p(), |
---|
678 | 0.5 * event[iOld].m(), 0.); |
---|
679 | int i2New = event.append( -idNewQ, 101, iOld, 0, 0, 0, |
---|
680 | 0, event[iOld].acol(), 0.5 * event[iOld].p(), |
---|
681 | 0.5 * event[iOld].m(), 0.); |
---|
682 | event[iOld].statusNeg(); |
---|
683 | event[iOld].daughters( i1New, i2New); |
---|
684 | |
---|
685 | // Partons in the two new systems. |
---|
686 | if (event[systemPtr->iParton[iFirst]].col() == event[i2New].acol()) |
---|
687 | swap( i1New, i2New); |
---|
688 | for (int i = 0; i <= iFirst; ++i) |
---|
689 | iNewSys1.push_back( systemPtr->iParton[i]); |
---|
690 | iNewSys1.push_back( i1New); |
---|
691 | iNewSys2.push_back( i2New); |
---|
692 | for (int i = iSecond; i < int(systemPtr->size()); ++i) |
---|
693 | iNewSys2.push_back( systemPtr->iParton[i]); |
---|
694 | |
---|
695 | // If several gluons between sparticles then find lowest-mass gluon pair |
---|
696 | // and replace it by a q-qbar pair. |
---|
697 | } else { |
---|
698 | |
---|
699 | // Find lowest-mass gluon pair and adjust effective quark mass. |
---|
700 | int iMin = 0; |
---|
701 | int i1Old = 0; |
---|
702 | int i2Old = 0; |
---|
703 | double mMin = 1e20; |
---|
704 | for (int i = iFirst + 1; i < iSecond - 1; ++i) { |
---|
705 | int i1Tmp = systemPtr->iParton[i]; |
---|
706 | int i2Tmp = systemPtr->iParton[i + 1]; |
---|
707 | double mTmp = (event[i1Tmp].p() + event[i2Tmp].p()).mCalc(); |
---|
708 | if (mTmp < mMin) { |
---|
709 | iMin = i; |
---|
710 | i1Old = i1Tmp; |
---|
711 | i2Old = i2Tmp; |
---|
712 | mMin = mTmp; |
---|
713 | } |
---|
714 | } |
---|
715 | double mEff = min( mNewQ + mOffsetCloudRH, 0.4 * mMin); |
---|
716 | |
---|
717 | // New kinematics by sharing mHat appropriately. |
---|
718 | Vec4 p1New, p2New; |
---|
719 | if ( !newKin( event[i1Old].p(), event[i2Old].p(), |
---|
720 | mEff, mEff, p1New, p2New, false) ) return false; |
---|
721 | |
---|
722 | // Insert new partons and update history. |
---|
723 | int i1New, i2New; |
---|
724 | if (event[systemPtr->iParton[0]].acol() == 0) { |
---|
725 | i1New = event.append( -idNewQ, 101, i1Old, 0, 0, 0, |
---|
726 | 0, event[i1Old].acol(), p1New, mEff, 0.); |
---|
727 | i2New = event.append( idNewQ, 101, i2Old, 0, 0, 0, |
---|
728 | event[i2Old].col(), 0, p2New, mEff, 0.); |
---|
729 | } else { |
---|
730 | i1New = event.append( idNewQ, 101, i1Old, 0, 0, 0, |
---|
731 | event[i1Old].col(), 0, p1New, mEff, 0.); |
---|
732 | i2New = event.append( -idNewQ, 101, i2Old, 0, 0, 0, |
---|
733 | 0, event[i2Old].acol(), p2New, mEff, 0.); |
---|
734 | } |
---|
735 | event[i1Old].statusNeg(); |
---|
736 | event[i2Old].statusNeg(); |
---|
737 | event[i1Old].daughters( i1New, 0); |
---|
738 | event[i2Old].daughters( i2New, 0); |
---|
739 | |
---|
740 | // Partons in the two new systems. |
---|
741 | for (int i = 0; i < iMin; ++i) |
---|
742 | iNewSys1.push_back( systemPtr->iParton[i]); |
---|
743 | iNewSys1.push_back( i1New); |
---|
744 | iNewSys2.push_back( i2New); |
---|
745 | for (int i = iMin + 2; i < int(systemPtr->size()); ++i) |
---|
746 | iNewSys2.push_back( systemPtr->iParton[i]); |
---|
747 | } |
---|
748 | |
---|
749 | // Erase the old system and insert the two new ones instead. |
---|
750 | colConfig.erase(iSys); |
---|
751 | colConfig.insert( iNewSys1, event); |
---|
752 | colConfig.insert( iNewSys2, event); |
---|
753 | |
---|
754 | // Done. |
---|
755 | return true; |
---|
756 | |
---|
757 | } |
---|
758 | |
---|
759 | //-------------------------------------------------------------------------- |
---|
760 | |
---|
761 | // Produce a R-hadron from a squark and another string end. |
---|
762 | |
---|
763 | bool RHadrons::produceSquark( ColConfig& colConfig, Event& event) { |
---|
764 | |
---|
765 | // Initial values. |
---|
766 | int nBody = 0; |
---|
767 | int iRNow = 0; |
---|
768 | int iNewQ = 0; |
---|
769 | int iNewL = 0; |
---|
770 | |
---|
771 | // Check at which end of the string the squark is located. |
---|
772 | int idAbsTop = event[ systemPtr->iParton[0] ].idAbs(); |
---|
773 | bool sqAtTop = (allowRSb && idAbsTop == idRSb) |
---|
774 | || (allowRSt && idAbsTop == idRSt); |
---|
775 | |
---|
776 | // Copy down system consecutively from squark end. |
---|
777 | int iBeg = event.size(); |
---|
778 | iCreRHad[iRHad] = iBeg; |
---|
779 | if (sqAtTop) for (int i = 0; i < systemPtr->size(); ++i) |
---|
780 | event.copy( systemPtr->iParton[i], 102); |
---|
781 | else for (int i = systemPtr->size() - 1; i >= 0; --i) |
---|
782 | event.copy( systemPtr->iParton[i], 102); |
---|
783 | int iEnd = event.size() - 1; |
---|
784 | |
---|
785 | // Input flavours. From now on H = heavy and L = light string ends. |
---|
786 | int idOldH = event[iBeg].id(); |
---|
787 | int idOldL = event[iEnd].id(); |
---|
788 | |
---|
789 | // Pick new flavour to form R-hadron. |
---|
790 | FlavContainer flavOld( idOldH%10); |
---|
791 | int idNewQ = flavSelPtr->pick(flavOld).id; |
---|
792 | int idRHad = toIdWithSquark( idOldH, idNewQ); |
---|
793 | if (idRHad == 0) { |
---|
794 | infoPtr->errorMsg("Error in RHadrons::produceSquark: " |
---|
795 | "cannot form R-hadron code"); |
---|
796 | return false; |
---|
797 | } |
---|
798 | |
---|
799 | // Target mass of R-hadron and z value of fragmentation function. |
---|
800 | double mRHad = particleDataPtr->m0(idRHad) + event[iBeg].m() |
---|
801 | - ( (abs(idOldH) == idRSb) ? m0Sb : m0St ); |
---|
802 | double z = zSelPtr->zFrag( idOldH, idNewQ, mRHad*mRHad); |
---|
803 | |
---|
804 | // Basic kinematics of string piece where break is to occur. |
---|
805 | Vec4 pOldH = event[iBeg].p(); |
---|
806 | int iOldL = iBeg + 1; |
---|
807 | Vec4 pOldL = event[iOldL].p(); |
---|
808 | double mOldL = event[iOldL].m(); |
---|
809 | double mNewH = mRHad / z; |
---|
810 | double sSys = (pOldH + pOldL).m2Calc(); |
---|
811 | double sRem = (1. - z) * (sSys - mNewH*mNewH); |
---|
812 | double sMin = pow2(mOldL + mCollapseRH); |
---|
813 | |
---|
814 | // If too low remaining mass in system then add one more parton to it. |
---|
815 | while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) ) |
---|
816 | && iOldL < iEnd ) { |
---|
817 | ++iOldL; |
---|
818 | pOldL += event[iOldL].p(); |
---|
819 | mOldL = event[iOldL].m(); |
---|
820 | sSys = (pOldH + pOldL).m2Calc(); |
---|
821 | sRem = (1. - z) * (sSys - mNewH*mNewH); |
---|
822 | sMin = pow2(mOldL + mCollapseRH); |
---|
823 | } |
---|
824 | |
---|
825 | // If enough mass then split off R-hadron and reduced system. |
---|
826 | if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) { |
---|
827 | Vec4 pNewH, pNewL; |
---|
828 | if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) { |
---|
829 | infoPtr->errorMsg("Error in RHadrons::produceSquark: " |
---|
830 | "failed to construct kinematics with reduced system"); |
---|
831 | return false; |
---|
832 | } |
---|
833 | |
---|
834 | // Insert R-hadron with new momentum. |
---|
835 | iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0, |
---|
836 | z * pNewH, mRHad, 0.); |
---|
837 | |
---|
838 | // Reduced system with new string endpoint and modified recoiler. |
---|
839 | idNewQ = -idNewQ; |
---|
840 | bool hasCol = (idNewQ > 0 && idNewQ < 10) || idNewQ < -10; |
---|
841 | int col = (hasCol) ? event[iOldL].acol() : 0; |
---|
842 | int acol = (hasCol) ? 0 : event[iOldL].col(); |
---|
843 | iNewQ = event.append( idNewQ, 105, iBeg, iOldL, 0, 0, col, acol, |
---|
844 | (1. - z) * pNewH, (1. - z) * mNewH, 0.); |
---|
845 | iNewL = event.copy( iOldL, 105); |
---|
846 | event[iNewL].mothers( iBeg, iOldL); |
---|
847 | event[iNewL].p( pNewL); |
---|
848 | |
---|
849 | // Done with processing of split to R-hadron and reduced system. |
---|
850 | nBody = 3; |
---|
851 | } |
---|
852 | |
---|
853 | // Two-body final state: form light hadron from endpoint and new flavour. |
---|
854 | if (nBody == 0) { |
---|
855 | FlavContainer flav1( idOldL); |
---|
856 | FlavContainer flav2( -idNewQ); |
---|
857 | int iTry = 0; |
---|
858 | int idNewL = flavSelPtr->combine( flav1, flav2); |
---|
859 | while (++iTry < NTRYMAX && idNewL == 0) |
---|
860 | idNewL = flavSelPtr->combine( flav1, flav2); |
---|
861 | if (idNewL == 0) { |
---|
862 | infoPtr->errorMsg("Error in RHadrons::produceSquark: " |
---|
863 | "cannot form light hadron code"); |
---|
864 | return false; |
---|
865 | } |
---|
866 | double mNewL = particleDataPtr->mass( idNewL); |
---|
867 | |
---|
868 | // Check that energy enough for light hadron and R-hadron. |
---|
869 | if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) { |
---|
870 | Vec4 pRHad, pNewL; |
---|
871 | if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) { |
---|
872 | infoPtr->errorMsg("Error in RHadrons::produceSquark: " |
---|
873 | "failed to construct kinematics for two-hadron decay"); |
---|
874 | return false; |
---|
875 | } |
---|
876 | |
---|
877 | // Insert R-hadron and light hadron. |
---|
878 | iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0, |
---|
879 | pRHad, mRHad, 0.); |
---|
880 | event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0, |
---|
881 | pNewL, mNewL, 0.); |
---|
882 | |
---|
883 | // Done for two-body case. |
---|
884 | nBody = 2; |
---|
885 | } |
---|
886 | } |
---|
887 | |
---|
888 | // Final case: for very low mass collapse to one single R-hadron. |
---|
889 | if (nBody == 0) { |
---|
890 | idRHad = toIdWithSquark( idOldH, idOldL); |
---|
891 | if (idRHad == 0) { |
---|
892 | infoPtr->errorMsg("Error in RHadrons::produceSquark: " |
---|
893 | "cannot form R-hadron code"); |
---|
894 | return false; |
---|
895 | } |
---|
896 | |
---|
897 | // Insert R-hadron with new momentum. |
---|
898 | iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0, |
---|
899 | systemPtr->pSum, systemPtr->mass, 0.); |
---|
900 | |
---|
901 | // Done with one-body case. |
---|
902 | nBody = 1; |
---|
903 | } |
---|
904 | |
---|
905 | // History bookkeeping: the R-hadron and processed partons. |
---|
906 | iRHadron[iRHad] = iRNow; |
---|
907 | int iLast = event.size() - 1; |
---|
908 | for (int i = iBeg; i <= iOldL; ++i) { |
---|
909 | event[i].statusNeg(); |
---|
910 | event[i].daughters( iRNow, iLast); |
---|
911 | } |
---|
912 | |
---|
913 | // Remove old string system and, if needed, insert a new one. |
---|
914 | colConfig.erase(iSys); |
---|
915 | if (nBody == 3) { |
---|
916 | vector<int> iNewSys; |
---|
917 | iNewSys.push_back( iNewQ); |
---|
918 | iNewSys.push_back( iNewL); |
---|
919 | for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys.push_back( i); |
---|
920 | colConfig.insert( iNewSys, event); |
---|
921 | } |
---|
922 | |
---|
923 | // Copy lifetime and vertex from sparticle to R-hadron. |
---|
924 | event[iRNow].tau( event[iBef].tau() ); |
---|
925 | if (event[iBef].hasVertex()) event[iRNow].vProd( event[iBef].vProd() ); |
---|
926 | |
---|
927 | // Done with production of a R-hadron from a squark. |
---|
928 | return true; |
---|
929 | |
---|
930 | } |
---|
931 | |
---|
932 | //-------------------------------------------------------------------------- |
---|
933 | |
---|
934 | // Produce a R-hadron from a gluino and two string ends (or a gluon). |
---|
935 | |
---|
936 | bool RHadrons::produceGluino( ColConfig& colConfig, Event& event) { |
---|
937 | |
---|
938 | // Initial values. |
---|
939 | int iGlui = 0; |
---|
940 | int idSave = 0; |
---|
941 | int idQLeap = 0; |
---|
942 | bool isDiq1 = false; |
---|
943 | vector<int> iSide1, iSide2, iSideTmp, iNewSys1, iNewSys2; |
---|
944 | Vec4 pGlui, pSide1, pSide2; |
---|
945 | |
---|
946 | // Decide whether to produce a gluinoball or not. |
---|
947 | bool isGBall = (rndmPtr->flat() < probGluinoballRH); |
---|
948 | |
---|
949 | // Extract one string system on either side of the gluino. |
---|
950 | for (int i = 0; i < systemPtr->size(); ++i) { |
---|
951 | int iTmp = systemPtr->iParton[i]; |
---|
952 | if (iGlui == 0 && event[ iTmp ].id() == idRGo) { |
---|
953 | iGlui = iTmp; |
---|
954 | pGlui = event[ iTmp ].p(); |
---|
955 | } else if (iGlui == 0) { |
---|
956 | iSideTmp.push_back( iTmp); |
---|
957 | pSide1 += event[ iTmp ].p(); |
---|
958 | } else { |
---|
959 | iSide2.push_back( iTmp); |
---|
960 | pSide2 += event[ iTmp ].p(); |
---|
961 | } |
---|
962 | } |
---|
963 | |
---|
964 | // Order sides from gluino outwards and with lowest relative mass first. |
---|
965 | for (int i = int(iSideTmp.size()) - 1; i >= 0; --i) |
---|
966 | iSide1.push_back( iSideTmp[i]); |
---|
967 | double m1H = (pGlui + pSide1).mCalc() - event[ iSide1.back() ].m(); |
---|
968 | double m2H = (pGlui + pSide2).mCalc() - event[ iSide2.back() ].m(); |
---|
969 | if (m2H < m1H) { |
---|
970 | swap( iSide1, iSide2); |
---|
971 | swap( pSide1, pSide2); |
---|
972 | } |
---|
973 | |
---|
974 | // Begin loop over two sides of gluino, with lowest relative mass first. |
---|
975 | for (int iSide = 1; iSide < 3; ++iSide) { |
---|
976 | |
---|
977 | // Begin processing of lower-mass string side. |
---|
978 | int idRHad = 0; |
---|
979 | double mRHad = 0.; |
---|
980 | int nBody = 0; |
---|
981 | int iRNow = 0; |
---|
982 | int iNewQ = 0; |
---|
983 | int iNewL = 0; |
---|
984 | int statusRHad = 0; |
---|
985 | |
---|
986 | // Copy down system consecutively from gluino end. |
---|
987 | int iBeg = event.size(); |
---|
988 | if (iSide == 1) { |
---|
989 | iCreRHad[iRHad] = iBeg; |
---|
990 | event.copy( iGlui, 102); |
---|
991 | for (int i = 0; i < int(iSide1.size()); ++i) |
---|
992 | event.copy( iSide1[i], 102); |
---|
993 | } else { |
---|
994 | event.copy( iGlui, 102); |
---|
995 | for (int i = 0; i < int(iSide2.size()); ++i) |
---|
996 | event.copy( iSide2[i], 102); |
---|
997 | } |
---|
998 | int iEnd = event.size() - 1; |
---|
999 | |
---|
1000 | // Pick new flavour to help form R-hadron. Initial colour values. |
---|
1001 | int idOldL = event[iEnd].id(); |
---|
1002 | int idNewQ = 21; |
---|
1003 | if (!isGBall) { |
---|
1004 | do { |
---|
1005 | FlavContainer flavOld( idOldL); |
---|
1006 | idNewQ = -flavSelPtr->pick(flavOld).id; |
---|
1007 | } while (iSide == 2 && isDiq1 && abs(idNewQ) > 10); |
---|
1008 | if (iSide == 1) isDiq1 = (abs(idNewQ) > 10); |
---|
1009 | } |
---|
1010 | bool hasCol = (event[iEnd].col() > 0); |
---|
1011 | int colR = 0; |
---|
1012 | int acolR = 0; |
---|
1013 | |
---|
1014 | // Target intermediary mass or R-hadron mass. |
---|
1015 | if (iSide == 1) { |
---|
1016 | idSave = idNewQ; |
---|
1017 | idRHad = (hasCol) ? 1009002 : -1009002; |
---|
1018 | if (hasCol) colR = event[iBeg].col(); |
---|
1019 | else acolR = event[iBeg].acol(); |
---|
1020 | statusRHad = 103; |
---|
1021 | double mNewQ = particleDataPtr->constituentMass( idNewQ); |
---|
1022 | if (isGBall) mNewQ *= 0.5; |
---|
1023 | mRHad = event[iBeg].m() + mOffsetCloudRH + mNewQ; |
---|
1024 | } else { |
---|
1025 | idRHad = toIdWithGluino( idSave, idNewQ); |
---|
1026 | if (idRHad == 0) { |
---|
1027 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1028 | "cannot form R-hadron code"); |
---|
1029 | return false; |
---|
1030 | } |
---|
1031 | statusRHad = 104; |
---|
1032 | mRHad = particleDataPtr->m0(idRHad) + event[iBeg].m() - m0Go; |
---|
1033 | } |
---|
1034 | |
---|
1035 | // z value of fragmentation function. |
---|
1036 | double z = zSelPtr->zFrag( idRGo, idNewQ, mRHad*mRHad); |
---|
1037 | |
---|
1038 | // Basic kinematics of string piece where break is to occur. |
---|
1039 | Vec4 pOldH = event[iBeg].p(); |
---|
1040 | int iOldL = iBeg + 1; |
---|
1041 | Vec4 pOldL = event[iOldL].p(); |
---|
1042 | double mOldL = event[iOldL].m(); |
---|
1043 | double mNewH = mRHad / z; |
---|
1044 | double sSys = (pOldH + pOldL).m2Calc(); |
---|
1045 | double sRem = (1. - z) * (sSys - mNewH*mNewH); |
---|
1046 | double sMin = pow2(mOldL + mCollapseRH); |
---|
1047 | |
---|
1048 | // If too low remaining mass in system then add one more parton to it. |
---|
1049 | while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) ) |
---|
1050 | && iOldL < iEnd ) { |
---|
1051 | ++iOldL; |
---|
1052 | pOldL += event[iOldL].p(); |
---|
1053 | mOldL = event[iOldL].m(); |
---|
1054 | sSys = (pOldH + pOldL).m2Calc(); |
---|
1055 | sRem = (1. - z) * (sSys - mNewH*mNewH); |
---|
1056 | sMin = pow2(mOldL + mCollapseRH); |
---|
1057 | } |
---|
1058 | |
---|
1059 | // If enough mass then split off R-hadron and reduced system. |
---|
1060 | if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) { |
---|
1061 | Vec4 pNewH, pNewL; |
---|
1062 | if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) { |
---|
1063 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1064 | "failed to construct kinematics with reduced system"); |
---|
1065 | return false; |
---|
1066 | } |
---|
1067 | |
---|
1068 | // Insert intermediary or R-hadron with new momentum, less colour. |
---|
1069 | iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, |
---|
1070 | 0, 0, colR, acolR, z * pNewH, mRHad, 0.); |
---|
1071 | |
---|
1072 | // Reduced system with new string endpoint and modified recoiler. |
---|
1073 | if (!isGBall) idNewQ = -idNewQ; |
---|
1074 | int colN = (hasCol) ? 0 : event[iOldL].acol(); |
---|
1075 | int acolN = (hasCol) ? event[iOldL].col() : 0; |
---|
1076 | iNewQ = event.append( idNewQ, 105, iBeg, iOldL, 0, 0, |
---|
1077 | colN, acolN, (1. - z) * pNewH, (1. - z) * mNewH, 0.); |
---|
1078 | iNewL = event.copy( iOldL, 105); |
---|
1079 | event[iNewL].mothers( iBeg, iOldL); |
---|
1080 | event[iNewL].p( pNewL); |
---|
1081 | |
---|
1082 | // Collect partons of new string system. |
---|
1083 | if (iSide == 1) { |
---|
1084 | iNewSys1.push_back( iNewQ); |
---|
1085 | iNewSys1.push_back( iNewL); |
---|
1086 | for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys1.push_back( i); |
---|
1087 | } else { |
---|
1088 | iNewSys2.push_back( iNewQ); |
---|
1089 | iNewSys2.push_back( iNewL); |
---|
1090 | for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys2.push_back( i); |
---|
1091 | } |
---|
1092 | |
---|
1093 | // Done with processing of split to R-hadron and reduced system. |
---|
1094 | nBody = 3; |
---|
1095 | } |
---|
1096 | |
---|
1097 | // For side-1 low-mass glueball system reabsorb full momentum. |
---|
1098 | if (nBody == 0 && isGBall && iSide == 1) { |
---|
1099 | idQLeap = event[iOldL].id(); |
---|
1100 | Vec4 pNewH = event[iBeg].p() + pOldL; |
---|
1101 | |
---|
1102 | // Insert intermediary R-hadron with new momentum, less colour. |
---|
1103 | iRNow = event.append( idRHad, statusRHad, iBeg, iEnd, |
---|
1104 | 0, 0, colR, acolR, pNewH, pNewH.mCalc(), 0.); |
---|
1105 | nBody = 1; |
---|
1106 | } |
---|
1107 | |
---|
1108 | // Two-body final state: form light hadron from endpoint and new flavour. |
---|
1109 | // Also for side 2 if gluinoball formation gives quark from side 1. |
---|
1110 | if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) { |
---|
1111 | if (isGBall) idNewQ = -idQLeap; |
---|
1112 | FlavContainer flav1( idOldL); |
---|
1113 | FlavContainer flav2( -idNewQ); |
---|
1114 | int iTry = 0; |
---|
1115 | int idNewL = flavSelPtr->combine( flav1, flav2); |
---|
1116 | while (++iTry < NTRYMAX && idNewL == 0) |
---|
1117 | idNewL = flavSelPtr->combine( flav1, flav2); |
---|
1118 | if (idNewL == 0) { |
---|
1119 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1120 | "cannot form light hadron code"); |
---|
1121 | return false; |
---|
1122 | } |
---|
1123 | double mNewL = particleDataPtr->mass( idNewL); |
---|
1124 | |
---|
1125 | // Check that energy enough for light hadron and R-hadron. |
---|
1126 | if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) { |
---|
1127 | Vec4 pRHad, pNewL; |
---|
1128 | if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) { |
---|
1129 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1130 | "failed to construct kinematics for two-hadron decay"); |
---|
1131 | return false; |
---|
1132 | } |
---|
1133 | |
---|
1134 | // Insert intermediary or R-hadron and light hadron. |
---|
1135 | iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0, |
---|
1136 | colR, acolR, pRHad, mRHad, 0.); |
---|
1137 | event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0, |
---|
1138 | pNewL, mNewL, 0.); |
---|
1139 | |
---|
1140 | // Done for two-body case. |
---|
1141 | nBody = 2; |
---|
1142 | // For this case gluinoball should be handled as normal flavour. |
---|
1143 | isGBall = false; |
---|
1144 | } |
---|
1145 | } |
---|
1146 | |
---|
1147 | // Final case: for very low mass collapse to one single R-hadron. |
---|
1148 | if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) { |
---|
1149 | if (isGBall) idSave = idQLeap; |
---|
1150 | if (iSide == 1) idSave = idOldL; |
---|
1151 | else idRHad = toIdWithGluino( idSave, idOldL); |
---|
1152 | if (idRHad == 0) { |
---|
1153 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1154 | "cannot form R-hadron code"); |
---|
1155 | return false; |
---|
1156 | } |
---|
1157 | |
---|
1158 | // Insert R-hadron with new momentum. |
---|
1159 | iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0, |
---|
1160 | colR, acolR, pOldH + pOldL, (pOldH + pOldL).mCalc(), 0.); |
---|
1161 | |
---|
1162 | // Done with one-body case. |
---|
1163 | // Even if hoped-for, it was not possible to create a gluinoball. |
---|
1164 | isGBall = false; |
---|
1165 | } |
---|
1166 | |
---|
1167 | // History bookkeeping: the processed partons. |
---|
1168 | int iLast = event.size() - 1; |
---|
1169 | for (int i = iBeg; i <= iOldL; ++i) { |
---|
1170 | event[i].statusNeg(); |
---|
1171 | event[i].daughters( iRNow, iLast); |
---|
1172 | } |
---|
1173 | |
---|
1174 | // End of loop over two sides of the gluino. |
---|
1175 | iGlui = iRNow; |
---|
1176 | } |
---|
1177 | |
---|
1178 | // History bookkeeping: insert R-hadron; delete old string system. |
---|
1179 | if (iGlui == 0) { |
---|
1180 | infoPtr->errorMsg("Error in RHadrons::produceGluino: " |
---|
1181 | "could not handle gluinoball kinematics"); |
---|
1182 | return false; |
---|
1183 | } |
---|
1184 | iRHadron[iRHad] = iGlui; |
---|
1185 | colConfig.erase(iSys); |
---|
1186 | |
---|
1187 | // Non-gluinoball: insert (up to) two new string systems. |
---|
1188 | if (!isGBall) { |
---|
1189 | if (iNewSys1.size() > 0) colConfig.insert( iNewSys1, event); |
---|
1190 | if (iNewSys2.size() > 0) colConfig.insert( iNewSys2, event); |
---|
1191 | |
---|
1192 | // Gluinoball with enough energy in first string: |
---|
1193 | // join two temporary gluons into one. |
---|
1194 | } else if (idQLeap == 0) { |
---|
1195 | int iG1 = iNewSys1[0]; |
---|
1196 | int iG2 = iNewSys2[0]; |
---|
1197 | int colG = event[iG1].col() + event[iG2].col(); |
---|
1198 | int acolG = event[iG1].acol() + event[iG2].acol(); |
---|
1199 | Vec4 pG = event[iG1].p() + event[iG2].p(); |
---|
1200 | int iG12 = event.append( 21, 107, iG1, iG2, 0, 0, colG, acolG, |
---|
1201 | pG, pG.mCalc(), 0.); |
---|
1202 | |
---|
1203 | // Temporary gluons no longer needed, but new colour to avoid warnings. |
---|
1204 | event[iG1].id( 21); |
---|
1205 | event[iG2].id( 21); |
---|
1206 | event[iG1].statusNeg(); |
---|
1207 | event[iG2].statusNeg(); |
---|
1208 | event[iG1].daughter1( iG12); |
---|
1209 | event[iG2].daughter1( iG12); |
---|
1210 | int colBridge = event.nextColTag(); |
---|
1211 | if (event[iG1].col() == 0) { |
---|
1212 | event[iG1].col( colBridge); |
---|
1213 | event[iG2].acol( colBridge); |
---|
1214 | } else { |
---|
1215 | event[iG1].acol( colBridge); |
---|
1216 | event[iG2].col( colBridge); |
---|
1217 | } |
---|
1218 | |
---|
1219 | // Form the remnant system from which the R-hadron has been split off. |
---|
1220 | vector<int> iNewSys12; |
---|
1221 | for (int i = int(iNewSys1.size()) - 1; i > 0; --i) |
---|
1222 | iNewSys12.push_back( iNewSys1[i]); |
---|
1223 | iNewSys12.push_back( iG12); |
---|
1224 | for (int i = 1; i < int(iNewSys2.size()); ++i) |
---|
1225 | iNewSys12.push_back( iNewSys2[i]); |
---|
1226 | colConfig.insert( iNewSys12, event); |
---|
1227 | |
---|
1228 | // Gluinoball where side 1 was fully eaten, and its flavour content |
---|
1229 | // should leap over to the other side, to a gluon there. |
---|
1230 | } else { |
---|
1231 | int iG2 = iNewSys2[0]; |
---|
1232 | event[iG2].id( idQLeap); |
---|
1233 | colConfig.insert( iNewSys2, event); |
---|
1234 | } |
---|
1235 | |
---|
1236 | // Copy lifetime and vertex from sparticle to R-hadron. |
---|
1237 | event[iGlui].tau( event[iBef].tau() ); |
---|
1238 | if (event[iBef].hasVertex()) event[iGlui].vProd( event[iBef].vProd() ); |
---|
1239 | |
---|
1240 | // Done with production of a R-hadron from a gluino. |
---|
1241 | return true; |
---|
1242 | |
---|
1243 | } |
---|
1244 | |
---|
1245 | //-------------------------------------------------------------------------- |
---|
1246 | |
---|
1247 | // Form a R-hadron code from a squark and a (di)quark code. |
---|
1248 | // First argument is the (anti)squark, second the (anti)(di)quark. |
---|
1249 | |
---|
1250 | int RHadrons::toIdWithSquark( int id1, int id2) { |
---|
1251 | |
---|
1252 | // Check that physical combination; return 0 if not. |
---|
1253 | int id1Abs = abs(id1); |
---|
1254 | int id2Abs = abs(id2); |
---|
1255 | if (id2Abs < 10 && id1 > 0 && id2 > 0) return 0; |
---|
1256 | if (id2Abs < 10 && id1 < 0 && id2 < 0) return 0; |
---|
1257 | if (id2Abs > 10 && id1 > 0 && id2 < 0) return 0; |
---|
1258 | if (id2Abs > 10 && id1 < 0 && id2 > 0) return 0; |
---|
1259 | |
---|
1260 | // Form R-hadron code. Flip sign for antisquark. |
---|
1261 | bool isSt = (id1Abs == idRSt); |
---|
1262 | int idRHad = 1000000; |
---|
1263 | if (id2Abs < 10) idRHad += ((isSt) ? 600 : 500) + 10 * id2Abs + 2; |
---|
1264 | else idRHad += ((isSt) ? 6000 : 5000) + 10 * (id2Abs/100) + id2Abs%10; |
---|
1265 | if (id1 < 0) idRHad = -idRHad; |
---|
1266 | |
---|
1267 | // Done. |
---|
1268 | return idRHad; |
---|
1269 | |
---|
1270 | } |
---|
1271 | |
---|
1272 | //-------------------------------------------------------------------------- |
---|
1273 | |
---|
1274 | // Resolve a R-hadron code into a squark and a (di)quark. |
---|
1275 | // Return a pair, where first is the squark and the second the (di)quark. |
---|
1276 | |
---|
1277 | pair<int,int> RHadrons::fromIdWithSquark( int idRHad) { |
---|
1278 | |
---|
1279 | // Find squark flavour content. |
---|
1280 | int idLight = (abs(idRHad) - 1000000) / 10; |
---|
1281 | int idSq = (idLight < 100) ? idLight/10 : idLight/100; |
---|
1282 | int id1 = (idSq == 6) ? idRSt : idRSb; |
---|
1283 | if (idRHad < 0) id1 = -id1; |
---|
1284 | |
---|
1285 | // Find light (di)quark flavour content. |
---|
1286 | int id2 = (idLight < 100) ? idLight%10 : idLight%100; |
---|
1287 | if (id2 > 10) id2 = 100 * id2 + abs(idRHad)%10; |
---|
1288 | if ((id2 < 10 && idRHad > 0) || (id2 > 10 && idRHad < 0)) id2 = -id2; |
---|
1289 | |
---|
1290 | // Done. |
---|
1291 | return make_pair( id1, id2); |
---|
1292 | |
---|
1293 | } |
---|
1294 | |
---|
1295 | //-------------------------------------------------------------------------- |
---|
1296 | |
---|
1297 | // Form a R-hadron code from two quark/diquark endpoints and a gluino. |
---|
1298 | |
---|
1299 | int RHadrons::toIdWithGluino( int id1, int id2) { |
---|
1300 | |
---|
1301 | // Check that physical combination; return 0 if not. Handle gluinoball. |
---|
1302 | int id1Abs = abs(id1); |
---|
1303 | int id2Abs = abs(id2); |
---|
1304 | if (id1Abs == 21 && id2Abs == 21) return 1000993; |
---|
1305 | int idMax = max( id1Abs, id2Abs); |
---|
1306 | int idMin = min( id1Abs, id2Abs); |
---|
1307 | if (idMin > 10) return 0; |
---|
1308 | if (idMax > 10 && id1 > 0 && id2 < 0) return 0; |
---|
1309 | if (idMax > 10 && id1 < 0 && id2 > 0) return 0; |
---|
1310 | if (idMax < 10 && id1 > 0 && id2 > 0) return 0; |
---|
1311 | if (idMax < 10 && id1 < 0 && id2 < 0) return 0; |
---|
1312 | |
---|
1313 | // Form R-meson code. Flip sign for antiparticle. |
---|
1314 | int idRHad = 0; |
---|
1315 | if (idMax < 10) { |
---|
1316 | idRHad = 1009003 + 100 * idMax + 10 * idMin; |
---|
1317 | if (idMin != idMax && idMax%2 == 1) { |
---|
1318 | if (id1Abs == idMax && id1 > 0) idRHad = -idRHad; |
---|
1319 | if (id2Abs == idMax && id2 > 0) idRHad = -idRHad; |
---|
1320 | } |
---|
1321 | if (idMin != idMax && idMax%2 == 0) { |
---|
1322 | if (id1Abs == idMax && id1 < 0) idRHad = -idRHad; |
---|
1323 | if (id2Abs == idMax && id2 < 0) idRHad = -idRHad; |
---|
1324 | } |
---|
1325 | |
---|
1326 | // Form R-baryon code. Flip sign for antiparticle. |
---|
1327 | } else { |
---|
1328 | int idA = idMax/1000; |
---|
1329 | int idB = (idMax/100)%10; |
---|
1330 | int idC = idMin; |
---|
1331 | if (idC > idB) swap( idB, idC); |
---|
1332 | if (idB > idA) swap( idA, idB); |
---|
1333 | if (idC > idB) swap( idB, idC); |
---|
1334 | idRHad = 1090004 + 1000 * idA + 100 * idB + 10 * idC; |
---|
1335 | if (id1 < 0) idRHad = -idRHad; |
---|
1336 | } |
---|
1337 | |
---|
1338 | // Done. |
---|
1339 | return idRHad; |
---|
1340 | |
---|
1341 | } |
---|
1342 | |
---|
1343 | //-------------------------------------------------------------------------- |
---|
1344 | |
---|
1345 | // Resolve a R-hadron code into two quark/diquark endpoints (and a gluino). |
---|
1346 | // Return a pair, where first carries colour and second anticolour. |
---|
1347 | |
---|
1348 | pair<int,int> RHadrons::fromIdWithGluino( int idRHad) { |
---|
1349 | |
---|
1350 | // Find light flavour content of R-hadron. |
---|
1351 | int idLight = (abs(idRHad) - 1000000) / 10; |
---|
1352 | int id1, id2, idTmp, idA, idB, idC; |
---|
1353 | |
---|
1354 | // Gluinoballs: split g into d dbar or u ubar. |
---|
1355 | if (idLight < 100) { |
---|
1356 | id1 = (rndmPtr->flat() < 0.5) ? 1 : 2; |
---|
1357 | id2 = -id1; |
---|
1358 | |
---|
1359 | // Gluino-meson: split into q + qbar. |
---|
1360 | } else if (idLight < 1000) { |
---|
1361 | id1 = (idLight / 10) % 10; |
---|
1362 | id2 = -(idLight % 10); |
---|
1363 | // Flip signs when first quark of down-type. |
---|
1364 | if (id1%2 == 1) { |
---|
1365 | idTmp = id1; |
---|
1366 | id1 = -id2; |
---|
1367 | id2 = -idTmp; |
---|
1368 | } |
---|
1369 | |
---|
1370 | // Gluino-baryon: split to q + qq (diquark). |
---|
1371 | // Pick diquark at random, except if c or b involved. |
---|
1372 | } else { |
---|
1373 | idA = (idLight / 100) % 10; |
---|
1374 | idB = (idLight / 10) % 10; |
---|
1375 | idC = idLight % 10; |
---|
1376 | double rndmQ = 3. * rndmPtr->flat(); |
---|
1377 | if (idA > 3) rndmQ = 0.5; |
---|
1378 | if (rndmQ < 1.) { |
---|
1379 | id1 = idA; |
---|
1380 | id2 = 1000 * idB + 100 * idC + 3; |
---|
1381 | if (idB != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2; |
---|
1382 | } else if (rndmQ < 2.) { |
---|
1383 | id1 = idB; |
---|
1384 | id2 = 1000 * idA + 100 * idC + 3; |
---|
1385 | if (idA != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2; |
---|
1386 | } else { |
---|
1387 | id1 = idC; |
---|
1388 | id2 = 1000 * idA + 100 * idB +3; |
---|
1389 | if (idA != idB && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2; |
---|
1390 | } |
---|
1391 | } |
---|
1392 | |
---|
1393 | // Flip signs for anti-R-hadron. |
---|
1394 | if (idRHad < 0) { |
---|
1395 | idTmp = id1; |
---|
1396 | id1 = -id2; |
---|
1397 | id2 = -idTmp; |
---|
1398 | } |
---|
1399 | |
---|
1400 | // Done. |
---|
1401 | return make_pair( id1, id2); |
---|
1402 | |
---|
1403 | } |
---|
1404 | |
---|
1405 | //-------------------------------------------------------------------------- |
---|
1406 | |
---|
1407 | // Construct modified four-vectors to match modified masses: |
---|
1408 | // minimal reshuffling of momentum along common axis. |
---|
1409 | // Note that last two arguments are used to provide output! |
---|
1410 | |
---|
1411 | bool RHadrons::newKin( Vec4 pOld1, Vec4 pOld2, double mNew1, double mNew2, |
---|
1412 | Vec4& pNew1, Vec4& pNew2, bool checkMargin) { |
---|
1413 | |
---|
1414 | // Squared masses in initial and final kinematics. |
---|
1415 | double sSum = (pOld1 + pOld2).m2Calc(); |
---|
1416 | double sOld1 = pOld1.m2Calc(); |
---|
1417 | double sOld2 = pOld2.m2Calc(); |
---|
1418 | double sNew1 = mNew1 * mNew1; |
---|
1419 | double sNew2 = mNew2 * mNew2; |
---|
1420 | |
---|
1421 | // Check that kinematically possible. |
---|
1422 | if (checkMargin && pow2(mNew1 + mNew2 + MSAFETY) > sSum) return false; |
---|
1423 | |
---|
1424 | // Transfer coefficients to give four-vectors with the new masses. |
---|
1425 | double lamOld = sqrt( pow2(sSum - sOld1 - sOld2) - 4. * sOld1 * sOld2 ); |
---|
1426 | double lamNew = sqrt( pow2(sSum - sNew1 - sNew2) - 4. * sNew1 * sNew2 ); |
---|
1427 | double move1 = (lamNew * (sSum - sOld1 + sOld2) |
---|
1428 | - lamOld * (sSum - sNew1 + sNew2)) / (2. * sSum * lamOld); |
---|
1429 | double move2 = (lamNew * (sSum + sOld1 - sOld2) |
---|
1430 | - lamOld * (sSum + sNew1 - sNew2)) / (2. * sSum * lamOld); |
---|
1431 | |
---|
1432 | // Construct final vectors. Done. |
---|
1433 | pNew1 = (1. + move1) * pOld1 - move2 * pOld2; |
---|
1434 | pNew2 = (1. + move2) * pOld2 - move1 * pOld1; |
---|
1435 | return true; |
---|
1436 | |
---|
1437 | } |
---|
1438 | |
---|
1439 | //========================================================================== |
---|
1440 | |
---|
1441 | } // end namespace Pythia8 |
---|