Line data Source code
1 : // RHadrons.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 RHadrons class.
7 :
8 : #include "Pythia8/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 0 : infoPtr = infoPtrIn;
55 0 : particleDataPtr = particleDataPtrIn;
56 0 : rndmPtr = rndmPtrIn;
57 :
58 : // Flags and parameters related to R-hadron formation and decay.
59 0 : allowRH = settings.flag("RHadrons:allow");
60 0 : maxWidthRH = settings.parm("RHadrons:maxWidth");
61 0 : idRSb = settings.mode("RHadrons:idSbottom");
62 0 : idRSt = settings.mode("RHadrons:idStop");
63 0 : idRGo = settings.mode("RHadrons:idGluino");
64 0 : setMassesRH = settings.flag("RHadrons:setMasses");
65 0 : probGluinoballRH = settings.parm("RHadrons:probGluinoball");
66 0 : mOffsetCloudRH = settings.parm("RHadrons:mOffsetCloud");
67 0 : mCollapseRH = settings.parm("RHadrons:mCollapse");
68 0 : diquarkSpin1RH = settings.parm("RHadrons:diquarkSpin1");
69 :
70 : // Check whether sbottom, stop or gluino should form R-hadrons.
71 0 : allowRSb = allowRH && idRSb > 0
72 0 : && (particleDataPtr->mWidth(idRSb) < maxWidthRH);
73 0 : allowRSt = allowRH && idRSt > 0
74 0 : && (particleDataPtr->mWidth(idRSt) < maxWidthRH);
75 0 : allowRGo = allowRH && idRGo > 0
76 0 : && (particleDataPtr->mWidth(idRGo) < maxWidthRH);
77 0 : allowSomeR = allowRSb || allowRSt || allowRGo;
78 :
79 : // Set nominal masses of sbottom R-mesons and R-baryons.
80 0 : if (allowRSb) {
81 0 : m0Sb = particleDataPtr->m0(idRSb);
82 0 : if (setMassesRH) {
83 0 : for (int i = 0; i < 14; ++i) {
84 0 : int idR = IDRHADSB[i];
85 0 : double m0RHad = m0Sb + mOffsetCloudRH;
86 0 : m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
87 0 : if (i > 4)
88 0 : m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
89 0 : particleDataPtr->m0( idR, m0RHad);
90 : }
91 0 : }
92 :
93 : // Set widths and lifetimes of sbottom R-hadrons.
94 0 : double mWidthRHad = particleDataPtr->mWidth(idRSb);
95 0 : double tau0RHad = particleDataPtr->tau0( idRSb);
96 0 : for (int i = 0; i < 14; ++i) {
97 0 : particleDataPtr->mWidth( IDRHADSB[i], mWidthRHad);
98 0 : particleDataPtr->tau0( IDRHADSB[i], tau0RHad);
99 : }
100 0 : }
101 :
102 : // Set nominal masses of stop R-mesons and R-baryons.
103 0 : if (allowRSt) {
104 0 : m0St = particleDataPtr->m0(idRSt);
105 0 : if (setMassesRH) {
106 0 : for (int i = 0; i < 14; ++i) {
107 0 : int idR = IDRHADST[i];
108 0 : double m0RHad = m0St + mOffsetCloudRH;
109 0 : m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
110 0 : if (i > 4)
111 0 : m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
112 0 : particleDataPtr->m0( idR, m0RHad);
113 : }
114 0 : }
115 :
116 : // Set widths and lifetimes of stop R-hadrons.
117 0 : double mWidthRHad = particleDataPtr->mWidth(idRSt);
118 0 : double tau0RHad = particleDataPtr->tau0( idRSt);
119 0 : for (int i = 0; i < 14; ++i) {
120 0 : particleDataPtr->mWidth( IDRHADST[i], mWidthRHad);
121 0 : particleDataPtr->tau0( IDRHADST[i], tau0RHad);
122 : }
123 0 : }
124 :
125 : // Set nominal masses of gluino R-glueballs, R-mesons and R-baryons.
126 0 : if (allowRGo) {
127 0 : m0Go = particleDataPtr->m0(idRGo);
128 0 : if (setMassesRH) {
129 0 : particleDataPtr->m0( IDRHADGO[0], m0Go + 2. * mOffsetCloudRH
130 0 : + particleDataPtr->constituentMass(21) );
131 0 : for (int i = 1; i < 38; ++i) {
132 0 : int idR = IDRHADGO[i];
133 0 : double m0RHad = m0Go + 2. * mOffsetCloudRH;
134 0 : m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
135 0 : m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
136 0 : if (i > 15)
137 0 : m0RHad += particleDataPtr->constituentMass( (idR%10000)/1000 );
138 0 : particleDataPtr->m0( idR, m0RHad);
139 : }
140 0 : }
141 :
142 : // Set widths and lifetimes of gluino R-hadrons.
143 0 : double mWidthRHad = particleDataPtr->mWidth(idRGo);
144 0 : double tau0RHad = particleDataPtr->tau0( idRGo);
145 0 : for (int i = 0; i < 38; ++i) {
146 0 : particleDataPtr->mWidth( IDRHADGO[i], mWidthRHad);
147 0 : particleDataPtr->tau0( IDRHADGO[i], tau0RHad);
148 : }
149 0 : }
150 :
151 : // Done.
152 0 : return true;
153 :
154 0 : }
155 : //--------------------------------------------------------------------------
156 :
157 : // Check if a given particle can produce R-hadrons.
158 :
159 : bool RHadrons::givesRHadron( int id) {
160 0 : if (allowRSb && abs(id) == idRSb) return true;
161 0 : if (allowRSt && abs(id) == idRSt) return true;
162 0 : if (allowRGo && id == idRGo) return true;
163 0 : return false;
164 :
165 0 : }
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 0 : if (!allowSomeR) return true;
175 :
176 : // Reset arrays and values.
177 0 : iBefRHad.resize(0);
178 0 : iCreRHad.resize(0);
179 0 : iRHadron.resize(0);
180 0 : iAftRHad.resize(0);
181 0 : isTriplet.resize(0);
182 0 : nRHad = 0;
183 :
184 : // Loop over event and identify hadronizing sparticles.
185 0 : for (int i = 0; i < event.size(); ++i)
186 0 : if (event[i].isFinal() && givesRHadron(event[i].id())) {
187 0 : iBefRHad.push_back(i);
188 0 : iCreRHad.push_back(i);
189 0 : iRHadron.push_back(0);
190 0 : iAftRHad.push_back(0);
191 0 : isTriplet.push_back(true);
192 0 : }
193 0 : nRHad = iRHadron.size();
194 :
195 : // Done if no hadronizing sparticles.
196 0 : if (nRHad == 0) return true;
197 :
198 : // Max two R-hadrons. Randomize order of processing.
199 0 : if (nRHad > 2) {
200 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
201 : "cannot handle more than two R-hadrons");
202 0 : return false;
203 : }
204 0 : 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 0 : iBef = iBefRHad[0];
208 0 : iSys = colConfig.findSinglet( iBef);
209 0 : systemPtr = &colConfig[iSys];
210 0 : if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) {
211 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
212 : "cannot handle system with junction");
213 0 : return false;
214 : }
215 0 : if (nRHad == 2) {
216 0 : iBef = iBefRHad[1];
217 0 : iSys = colConfig.findSinglet( iBefRHad[1]);
218 0 : systemPtr = &colConfig[iSys];
219 0 : if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) {
220 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
221 : "cannot handle system with junction");
222 0 : return false;
223 : }
224 : }
225 :
226 : // Open up a closed gluon/gluino loop.
227 0 : iBef = iBefRHad[0];
228 0 : iSys = colConfig.findSinglet( iBef);
229 0 : systemPtr = &colConfig[iSys];
230 0 : if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) {
231 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
232 : "cannot open up closed gluon/gluino loop");
233 0 : return false;
234 : }
235 0 : if (nRHad == 2) {
236 0 : iBef = iBefRHad[1];
237 0 : iSys = colConfig.findSinglet( iBefRHad[1]);
238 0 : systemPtr = &colConfig[iSys];
239 0 : if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) {
240 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
241 : "cannot open up closed gluon/gluino loop");
242 0 : return false;
243 : }
244 : }
245 :
246 : // Split up a colour singlet system that should give two R-hadrons.
247 0 : if (nRHad == 2) {
248 0 : int iSys1 = colConfig.findSinglet( iBefRHad[0]);
249 0 : int iSys2 = colConfig.findSinglet( iBefRHad[1]);
250 0 : if (iSys2 == iSys1) {
251 0 : iSys = iSys1;
252 0 : systemPtr = &colConfig[iSys];
253 0 : if ( !splitSystem( colConfig, event) ) {
254 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
255 : "failed to handle two sparticles in same system");
256 0 : return false;
257 : }
258 : }
259 0 : }
260 :
261 : // Loop over R-hadrons to be formed. Find its colour singlet system.
262 0 : for (iRHad = 0; iRHad < nRHad; ++iRHad) {
263 0 : iBef = iBefRHad[iRHad];
264 0 : iSys = colConfig.findSinglet( iBef);
265 0 : if (iSys < 0) {
266 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
267 : "sparticle not in any colour singlet");
268 0 : return false;
269 : }
270 0 : systemPtr = &colConfig[iSys];
271 :
272 : // For now don't handle systems involving junctions or loops.
273 0 : if (systemPtr->hasJunction) {
274 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
275 : "cannot handle system with junction");
276 0 : return false;
277 : }
278 0 : if (systemPtr->isClosed) {
279 0 : infoPtr->errorMsg("Error in RHadrons::produce: "
280 : "cannot handle closed colour loop");
281 0 : return false;
282 : }
283 :
284 : // Handle formation of R-hadron separately for gluino and squark.
285 0 : if (event[iBef].id() == idRGo) isTriplet[iRHad] = false;
286 0 : bool formed = (isTriplet[iRHad]) ? produceSquark( colConfig, event)
287 0 : : produceGluino( colConfig, event);
288 0 : if (!formed) return false;
289 :
290 : // End of loop over R-hadrons. Done.
291 0 : }
292 0 : return true;
293 :
294 0 : }
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 0 : for (iRHad = 0; iRHad < nRHad; ++iRHad) {
305 0 : int iRNow = iRHadron[iRHad];
306 0 : int iRBef = iBefRHad[iRHad];
307 0 : int idRHad = event[iRNow].id();
308 0 : double mRHad = event[iRNow].m();
309 0 : 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 0 : pair<int,int> idPair = (isTriplet[iRHad])
315 0 : ? 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 0 : double fracR = mRBef / mRHad;
322 0 : if (fracR >= 1.) {
323 0 : infoPtr->errorMsg("Error in RHadrons::decay: "
324 : "too low R-hadron mass for decay");
325 0 : return false;
326 : }
327 :
328 : // Squark: new colour needed in the breakup.
329 0 : if (isTriplet[iRHad]) {
330 0 : int colNew = event.nextColTag();
331 0 : int col = (event[iRBef].col() != 0) ? colNew : 0;
332 0 : int acol = (col == 0) ? colNew : 0;
333 :
334 : // Store the constituents of a squark R-hadron.
335 0 : iR0 = event.append( id1, 106, iRNow, 0, 0, 0, col, acol,
336 0 : fracR * event[iRNow].p(), fracR * mRHad, 0.);
337 0 : iR2 = event.append( id2, 106, iRNow, 0, 0, 0, acol, col,
338 0 : (1. - fracR) * event[iRNow].p(), (1. - fracR) * mRHad, 0.);
339 :
340 : // Gluino: set mass sharing between two spectators.
341 0 : } else {
342 0 : double m1Eff = particleDataPtr->constituentMass(id1) + mOffsetCloudRH;
343 0 : double m2Eff = particleDataPtr->constituentMass(id2) + mOffsetCloudRH;
344 0 : double frac1 = (1. - fracR) * m1Eff / ( m1Eff + m2Eff);
345 0 : double frac2 = (1. - fracR) * m2Eff / ( m1Eff + m2Eff);
346 :
347 : // Two new colours needed in the breakups.
348 0 : int col1 = event.nextColTag();
349 0 : int col2 = event.nextColTag();
350 :
351 : // Store the constituents of a gluino R-hadron.
352 0 : iR0 = event.append( idRGo, 106, iRNow, 0, 0, 0, col2, col1,
353 0 : fracR * event[iRNow].p(), fracR * mRHad, 0.);
354 0 : event.append( id1, 106, iRNow, 0, 0, 0, col1, 0,
355 0 : frac1 * event[iRNow].p(), frac1 * mRHad, 0.);
356 0 : iR2 = event.append( id2, 106, iRNow, 0, 0, 0, 0, col2,
357 0 : frac2 * event[iRNow].p(), frac2 * mRHad, 0.);
358 : }
359 :
360 : // Mark R-hadron as decayed and update history.
361 0 : event[iRNow].statusNeg();
362 0 : event[iRNow].daughters( iR0, iR2);
363 0 : iAftRHad[iRHad] = iR0;
364 :
365 : // Set secondary vertex for decay products, but no lifetime.
366 0 : Vec4 vDec = event[iRNow].vProd() + event[iRNow].tau()
367 0 : * event[iR0].p() / event[iR0].m();
368 0 : 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 0 : }
373 :
374 : // Done.
375 0 : return true;
376 :
377 0 : }
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 0 : vector<int> leg1, leg2, leg3;
387 : int iLegSP = 0;
388 : int iPosSP = 0;
389 : int iLeg = 0;
390 : int iPos = 0;
391 0 : for (int i = 0; i < systemPtr->size(); ++i) {
392 0 : ++iPos;
393 0 : int iP = systemPtr->iParton[i];
394 0 : if (iP < 0) {
395 0 : ++iLeg;
396 : iPos = -1;
397 0 : } else if (iLeg == 1) leg1.push_back( iP);
398 0 : else if (iLeg == 2) leg2.push_back( iP);
399 0 : else if (iLeg == 3) leg3.push_back( iP);
400 0 : if (iP == iBef) {
401 : iLegSP = iLeg;
402 : iPosSP = iPos;
403 0 : }
404 0 : }
405 0 : 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 0 : if (iLegSP == 2) swap(leg2, leg1);
410 0 : else if (iLegSP == 3) swap(leg3, leg1);
411 0 : for (int i = 0; i < iPosSP; ++i)
412 0 : if (event[leg1[i]].id() != 21) return true;
413 :
414 : // No gluon between sparticle and junction: find kinetic energy of system.
415 0 : if (iPosSP == 0) {
416 0 : Vec4 pSP = event[iBef].p();
417 0 : Vec4 pRec = 0.;
418 0 : for (int i = 0; i < int(leg2.size()); ++i) pRec += event[leg2[i]].p();
419 0 : for (int i = 0; i < int(leg3.size()); ++i) pRec += event[leg3[i]].p();
420 0 : double mSys = (pSP + pRec).mCalc();
421 0 : double mSP = pSP.mCalc();
422 0 : double mRec = pRec.mCalc();
423 0 : double eKin = mSys - mSP - mRec;
424 :
425 : // Insert new gluon that borrows part of kinetic energy.
426 0 : double mNewG = EGBORROWMAX * eKin / (EGBORROWMAX + eKin) ;
427 0 : Vec4 pNewG = (mNewG / mSys) * (pSP + pRec);
428 0 : int newCol = event.nextColTag();
429 0 : bool isCol = (event[leg1.back()].col() > 0);
430 0 : int colNG = (isCol)? newCol : event[iBef].acol();
431 0 : int acolNG = (isCol) ? event[iBef].col() : newCol;
432 0 : int iNewG = event.append( 21, 101, iBef, 0, 0, 0, colNG, acolNG,
433 0 : pNewG, mNewG, 0.);
434 0 : leg1.insert( leg1.begin(), iNewG);
435 0 : ++iPosSP;
436 :
437 : // Boosts for remainder systems that gave up energy.
438 0 : double mNewSys = mSys - mNewG;
439 0 : double pAbsOld = 0.5 * sqrtpos( pow2(mSys*mSys - mSP*mSP - mRec*mRec)
440 0 : - pow2(2. * mSP * mRec) ) / mSys;
441 0 : double pAbsNew = 0.5 * sqrtpos( pow2(mNewSys*mNewSys - mSP*mSP - mRec*mRec)
442 0 : - pow2(2. * mSP * mRec) ) / mNewSys;
443 0 : RotBstMatrix MtoCM, MfromCM, MSP, MRec;
444 0 : MtoCM.toCMframe( pSP, pRec);
445 0 : MfromCM = MtoCM;
446 0 : MfromCM.invert();
447 0 : MSP = MtoCM;
448 0 : MSP.bst( 0., 0., -pAbsOld / sqrt(mSP * mSP + pAbsOld * pAbsOld));
449 0 : MSP.bst( 0., 0., pAbsNew / sqrt(mSP * mSP + pAbsNew * pAbsNew));
450 0 : MSP.rotbst( MfromCM);
451 0 : MRec = MtoCM;
452 0 : MRec.bst( 0., 0., pAbsOld / sqrt(mRec * mRec + pAbsOld * pAbsOld));
453 0 : MRec.bst( 0., 0., -pAbsNew / sqrt(mRec * mRec + pAbsNew * pAbsNew));
454 0 : MRec.rotbst( MfromCM);
455 :
456 : // Copy down recoling partons and boost their momenta.
457 0 : int iNewSP = event.copy( iBef, 101);
458 0 : event[iNewSP].mother2(0);
459 0 : event[iBef].daughter1(iNewG);
460 0 : event[iNewSP].rotbst( MSP);
461 0 : leg1[iPosSP] = iNewSP;
462 0 : if (iBefRHad[0] == iBef) iBefRHad[0] = iNewSP;
463 0 : else if (nRHad > 1 && iBefRHad[1] == iBef) iBefRHad[1] = iNewSP;
464 0 : iBef = iNewSP;
465 0 : for (int i = 0; i < int(leg2.size()); ++i) {
466 0 : int iCopy = event.copy( leg2[i], 101);
467 0 : event[iCopy].rotbst( MRec);
468 0 : if (iBefRHad[0] == leg2[i]) iBefRHad[0] = iCopy;
469 0 : else if (nRHad > 1 && iBefRHad[1] == leg2[i]) iBefRHad[1] = iCopy;
470 0 : leg2[i] = iCopy;
471 : }
472 0 : for (int i = 0; i < int(leg3.size()); ++i) {
473 0 : int iCopy = event.copy( leg3[i], 101);
474 0 : event[iCopy].rotbst( MRec);
475 0 : if (iBefRHad[0] == leg3[i]) iBefRHad[0] = iCopy;
476 0 : else if (nRHad > 1 && iBefRHad[1] == leg3[i]) iBefRHad[1] = iCopy;
477 0 : leg3[i] = iCopy;
478 : }
479 :
480 : // Now always at least one gluon between sparticle and junction.
481 0 : }
482 :
483 : // Find gluon with largest energy in sparticle rest frame.
484 : int iGspl = 0;
485 0 : double eGspl = event[leg1[0]].p() * event[iBef].p();
486 0 : for (int i = 1; i < iPosSP; ++i) {
487 0 : double eGnow = event[leg1[i]].p() * event[iBef].p();
488 0 : if (eGnow > eGspl) {
489 : iGspl = i;
490 : eGspl = eGnow;
491 0 : }
492 : }
493 0 : int iG = leg1[iGspl];
494 :
495 : // Split this gluon into a collinear quark.antiquark pair.
496 0 : int idNewQ = flavSelPtr->pickLightQ();
497 0 : int iNewQ = event.append( idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0,
498 0 : 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
499 0 : int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(),
500 0 : 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
501 0 : event[iG].statusNeg();
502 0 : event[iG].daughters( iNewQ, iNewQb);
503 0 : if (event[leg1.back()].col() == 0) swap( iNewQ, iNewQb);
504 :
505 : // Collect two new systems after split.
506 0 : vector<int> iNewSys1, iNewSys2;
507 0 : iNewSys1.push_back( iNewQb);
508 0 : for (int i = iGspl + 1; i < int(leg1.size()); ++i)
509 0 : iNewSys1.push_back( leg1[i]);
510 0 : iNewSys2.push_back( -10);
511 0 : for (int i = 0; i < iGspl; ++i) iNewSys2.push_back( leg1[i]);
512 0 : iNewSys2.push_back( iNewQ);
513 0 : iNewSys2.push_back( -11);
514 0 : for (int i = 0; i < int(leg2.size()); ++i) iNewSys2.push_back( leg2[i]);
515 0 : iNewSys2.push_back( -12);
516 0 : 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 0 : colConfig.erase(iSys);
520 0 : colConfig.insert( iNewSys1, event);
521 0 : colConfig.insert( iNewSys2, event);
522 :
523 : // Done.
524 : return true;
525 :
526 0 : }
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 0 : for (int i = 0; i < systemPtr->size(); ++i) {
538 0 : int iGNow = systemPtr->iParton[i];
539 0 : if (event[iGNow].id() == 21) {
540 0 : double eGnow = event[iGNow].p() * event[iBef].p();
541 0 : if (eGnow > eGspl) {
542 : iGspl = i;
543 : eGspl = eGnow;
544 0 : }
545 0 : }
546 : }
547 0 : if (iGspl == -1) return false;
548 :
549 : // Split this gluon into a collinear quark.antiquark pair.
550 0 : int iG = systemPtr->iParton[iGspl];
551 0 : int idNewQ = flavSelPtr->pickLightQ();
552 0 : int iNewQ = event.append( idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0,
553 0 : 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
554 0 : int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(),
555 0 : 0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
556 0 : event[iG].statusNeg();
557 0 : event[iG].daughters( iNewQ, iNewQb);
558 :
559 : // Order partons in new system. Note order of colour flow.
560 0 : int iNext = iGspl + 1;
561 0 : if (iNext == systemPtr->size()) iNext = 0;
562 0 : if (event[ systemPtr->iParton[iNext]].acol() != event[iNewQ].col())
563 0 : swap( iNewQ, iNewQb);
564 0 : vector<int> iNewSys;
565 0 : iNewSys.push_back( iNewQ);
566 0 : for (int i = iGspl + 1; i < systemPtr->size(); ++i)
567 0 : iNewSys.push_back( systemPtr->iParton[i]);
568 0 : for (int i = 0; i < iGspl; ++i)
569 0 : iNewSys.push_back( systemPtr->iParton[i]);
570 0 : iNewSys.push_back( iNewQb);
571 :
572 : // Erase the old system and insert the new one instead.
573 0 : colConfig.erase(iSys);
574 0 : colConfig.insert( iNewSys, event);
575 :
576 : // Done.
577 : return true;
578 :
579 0 : }
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 0 : for (int i = 0; i < int(systemPtr->size()); ++i) {
592 0 : int iTmp = systemPtr->iParton[i];
593 0 : if ( givesRHadron( event[iTmp].id()) ) {
594 0 : if (iFirst == -1) iFirst = i;
595 : else iSecond = i;
596 : }
597 : }
598 0 : int nLeg = iSecond - iFirst;
599 :
600 : // New flavour pair for breaking the string, and its mass.
601 0 : int idNewQ = flavSelPtr->pickLightQ();
602 0 : double mNewQ = particleDataPtr->constituentMass( idNewQ);
603 0 : vector<int> iNewSys1, iNewSys2;
604 :
605 : // If sparticles are neighbours then need new q-qbar pair in between.
606 0 : if (nLeg == 1) {
607 :
608 : // The location of the two sparticles.
609 0 : int i1Old = systemPtr->iParton[iFirst];
610 0 : int i2Old = systemPtr->iParton[iSecond];
611 :
612 : // Take a fraction of momentum to give breakup quark pair.
613 0 : double mHat = (event[i1Old].p() + event[i2Old].p()).mCalc();
614 0 : double mMax = mHat - event[i1Old].m() - event[i2Old].m();
615 0 : if (mMax < 2. * (mNewQ + MSAFETY)) return false;
616 0 : double mEff = min( 2. * (mNewQ + mOffsetCloudRH), mMax - 2. * MSAFETY);
617 0 : double frac = mEff / mHat;
618 0 : 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 0 : Vec4 p1New, p2New;
623 0 : if ( !newKin( event[i1Old].p(), event[i2Old].p(),
624 0 : event[i1Old].m() / (1. - frac), event[i2Old].m() / (1. - frac),
625 0 : p1New, p2New) ) return false;
626 0 : p1New *= 1. - frac;
627 0 : p2New *= 1. - frac;
628 :
629 : // Fill in new partons after branching, with correct colour/flavour sign.
630 0 : int i1New, i2New, i3New, i4New;
631 0 : int newCol = event.nextColTag();
632 0 : i1New = event.copy( i1Old, 101);
633 0 : if (event[i2Old].acol() == event[i1Old].col()) {
634 0 : i3New = event.append( -idNewQ, 101, i1Old, 0, 0, 0,
635 0 : 0, event[i2Old].acol(), 0.5 * pEff, 0.5 * mEff, 0.);
636 0 : i2New = event.copy( i2Old, 101);
637 0 : event[i2New].acol( newCol);
638 0 : i4New = event.append( idNewQ, 101, i2Old, 0, 0, 0,
639 0 : newCol, 0, 0.5 * pEff, 0.5 * mEff, 0.);
640 0 : } else {
641 0 : i3New = event.append( idNewQ, 101, i1Old, 0, 0, 0,
642 0 : event[i2Old].col(), 0, 0.5 * pEff, 0.5 * mEff, 0.);
643 0 : i2New = event.copy( i2Old, 101);
644 0 : event[i2New].col( newCol);
645 0 : i4New = event.append( -idNewQ, 101, i2Old, 0, 0, 0,
646 0 : 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 0 : event[i1New].p( p1New);
652 0 : event[i2New].p( p2New);
653 0 : event[i1Old].daughters( i1New, i3New);
654 0 : event[i1New].mother2( 0);
655 0 : event[i2Old].daughters( i2New, i4New);
656 0 : event[i2New].mother2( 0);
657 0 : iBefRHad[0] = i1New;
658 0 : iBefRHad[1] = i2New;
659 :
660 : // Partons in the two new systems.
661 0 : for (int i = 0; i < iFirst; ++i)
662 0 : iNewSys1.push_back( systemPtr->iParton[i]);
663 0 : iNewSys1.push_back( i1New);
664 0 : iNewSys1.push_back( i3New);
665 0 : iNewSys2.push_back( i4New);
666 0 : iNewSys2.push_back( i2New);
667 0 : for (int i = iSecond + 1; i < int(systemPtr->size()); ++i)
668 0 : iNewSys2.push_back( systemPtr->iParton[i]);
669 :
670 : // If one gluon between sparticles then split it into a
671 : // collinear q-qbar pair.
672 0 : } else if (nLeg == 2) {
673 :
674 : // Gluon to be split and its two daughters.
675 0 : int iOld = systemPtr->iParton[iFirst + 1];
676 0 : int i1New = event.append( idNewQ, 101, iOld, 0, 0, 0,
677 0 : event[iOld].col(), 0, 0.5 * event[iOld].p(),
678 0 : 0.5 * event[iOld].m(), 0.);
679 0 : int i2New = event.append( -idNewQ, 101, iOld, 0, 0, 0,
680 0 : 0, event[iOld].acol(), 0.5 * event[iOld].p(),
681 0 : 0.5 * event[iOld].m(), 0.);
682 0 : event[iOld].statusNeg();
683 0 : event[iOld].daughters( i1New, i2New);
684 :
685 : // Partons in the two new systems.
686 0 : if (event[systemPtr->iParton[iFirst]].col() == event[i2New].acol())
687 0 : swap( i1New, i2New);
688 0 : for (int i = 0; i <= iFirst; ++i)
689 0 : iNewSys1.push_back( systemPtr->iParton[i]);
690 0 : iNewSys1.push_back( i1New);
691 0 : iNewSys2.push_back( i2New);
692 0 : for (int i = iSecond; i < int(systemPtr->size()); ++i)
693 0 : 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 0 : } 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 0 : for (int i = iFirst + 1; i < iSecond - 1; ++i) {
705 0 : int i1Tmp = systemPtr->iParton[i];
706 0 : int i2Tmp = systemPtr->iParton[i + 1];
707 0 : double mTmp = (event[i1Tmp].p() + event[i2Tmp].p()).mCalc();
708 0 : if (mTmp < mMin) {
709 : iMin = i;
710 : i1Old = i1Tmp;
711 : i2Old = i2Tmp;
712 : mMin = mTmp;
713 0 : }
714 : }
715 0 : double mEff = min( mNewQ + mOffsetCloudRH, 0.4 * mMin);
716 :
717 : // New kinematics by sharing mHat appropriately.
718 0 : Vec4 p1New, p2New;
719 0 : if ( !newKin( event[i1Old].p(), event[i2Old].p(),
720 0 : mEff, mEff, p1New, p2New, false) ) return false;
721 :
722 : // Insert new partons and update history.
723 0 : int i1New, i2New;
724 0 : if (event[systemPtr->iParton[0]].acol() == 0) {
725 0 : i1New = event.append( -idNewQ, 101, i1Old, 0, 0, 0,
726 0 : 0, event[i1Old].acol(), p1New, mEff, 0.);
727 0 : i2New = event.append( idNewQ, 101, i2Old, 0, 0, 0,
728 0 : event[i2Old].col(), 0, p2New, mEff, 0.);
729 0 : } else {
730 0 : i1New = event.append( idNewQ, 101, i1Old, 0, 0, 0,
731 0 : event[i1Old].col(), 0, p1New, mEff, 0.);
732 0 : i2New = event.append( -idNewQ, 101, i2Old, 0, 0, 0,
733 0 : 0, event[i2Old].acol(), p2New, mEff, 0.);
734 : }
735 0 : event[i1Old].statusNeg();
736 0 : event[i2Old].statusNeg();
737 0 : event[i1Old].daughters( i1New, 0);
738 0 : event[i2Old].daughters( i2New, 0);
739 :
740 : // Partons in the two new systems.
741 0 : for (int i = 0; i < iMin; ++i)
742 0 : iNewSys1.push_back( systemPtr->iParton[i]);
743 0 : iNewSys1.push_back( i1New);
744 0 : iNewSys2.push_back( i2New);
745 0 : for (int i = iMin + 2; i < int(systemPtr->size()); ++i)
746 0 : iNewSys2.push_back( systemPtr->iParton[i]);
747 0 : }
748 :
749 : // Erase the old system and insert the two new ones instead.
750 0 : colConfig.erase(iSys);
751 0 : colConfig.insert( iNewSys1, event);
752 0 : colConfig.insert( iNewSys2, event);
753 :
754 : // Done.
755 0 : return true;
756 :
757 0 : }
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 0 : int iNewQ = 0;
769 0 : int iNewL = 0;
770 :
771 : // Check at which end of the string the squark is located.
772 0 : int idAbsTop = event[ systemPtr->iParton[0] ].idAbs();
773 0 : bool sqAtTop = (allowRSb && idAbsTop == idRSb)
774 0 : || (allowRSt && idAbsTop == idRSt);
775 :
776 : // Copy down system consecutively from squark end.
777 0 : int iBeg = event.size();
778 0 : iCreRHad[iRHad] = iBeg;
779 0 : if (sqAtTop) for (int i = 0; i < systemPtr->size(); ++i)
780 0 : event.copy( systemPtr->iParton[i], 102);
781 0 : else for (int i = systemPtr->size() - 1; i >= 0; --i)
782 0 : event.copy( systemPtr->iParton[i], 102);
783 0 : int iEnd = event.size() - 1;
784 :
785 : // Input flavours. From now on H = heavy and L = light string ends.
786 0 : int idOldH = event[iBeg].id();
787 0 : int idOldL = event[iEnd].id();
788 :
789 : // Pick new flavour to form R-hadron.
790 0 : FlavContainer flavOld( idOldH%10);
791 0 : int idNewQ = flavSelPtr->pick(flavOld).id;
792 0 : int idRHad = toIdWithSquark( idOldH, idNewQ);
793 0 : if (idRHad == 0) {
794 0 : infoPtr->errorMsg("Error in RHadrons::produceSquark: "
795 : "cannot form R-hadron code");
796 0 : return false;
797 : }
798 :
799 : // Target mass of R-hadron and z value of fragmentation function.
800 0 : double mRHad = particleDataPtr->m0(idRHad) + event[iBeg].m()
801 0 : - ( (abs(idOldH) == idRSb) ? m0Sb : m0St );
802 0 : double z = zSelPtr->zFrag( idOldH, idNewQ, mRHad*mRHad);
803 :
804 : // Basic kinematics of string piece where break is to occur.
805 0 : Vec4 pOldH = event[iBeg].p();
806 0 : int iOldL = iBeg + 1;
807 0 : Vec4 pOldL = event[iOldL].p();
808 0 : double mOldL = event[iOldL].m();
809 0 : double mNewH = mRHad / z;
810 0 : double sSys = (pOldH + pOldL).m2Calc();
811 0 : double sRem = (1. - z) * (sSys - mNewH*mNewH);
812 0 : double sMin = pow2(mOldL + mCollapseRH);
813 :
814 : // If too low remaining mass in system then add one more parton to it.
815 0 : while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) )
816 0 : && iOldL < iEnd ) {
817 0 : ++iOldL;
818 0 : pOldL += event[iOldL].p();
819 0 : mOldL = event[iOldL].m();
820 0 : sSys = (pOldH + pOldL).m2Calc();
821 0 : sRem = (1. - z) * (sSys - mNewH*mNewH);
822 0 : sMin = pow2(mOldL + mCollapseRH);
823 : }
824 :
825 : // If enough mass then split off R-hadron and reduced system.
826 0 : if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) {
827 0 : Vec4 pNewH, pNewL;
828 0 : if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) {
829 0 : infoPtr->errorMsg("Error in RHadrons::produceSquark: "
830 : "failed to construct kinematics with reduced system");
831 0 : return false;
832 : }
833 :
834 : // Insert R-hadron with new momentum.
835 0 : iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
836 0 : z * pNewH, mRHad, 0.);
837 :
838 : // Reduced system with new string endpoint and modified recoiler.
839 0 : idNewQ = -idNewQ;
840 0 : bool hasCol = (idNewQ > 0 && idNewQ < 10) || idNewQ < -10;
841 0 : int col = (hasCol) ? event[iOldL].acol() : 0;
842 0 : int acol = (hasCol) ? 0 : event[iOldL].col();
843 0 : iNewQ = event.append( idNewQ, 105, iBeg, iOldL, 0, 0, col, acol,
844 0 : (1. - z) * pNewH, (1. - z) * mNewH, 0.);
845 0 : iNewL = event.copy( iOldL, 105);
846 0 : event[iNewL].mothers( iBeg, iOldL);
847 0 : event[iNewL].p( pNewL);
848 :
849 : // Done with processing of split to R-hadron and reduced system.
850 : nBody = 3;
851 0 : }
852 :
853 : // Two-body final state: form light hadron from endpoint and new flavour.
854 0 : if (nBody == 0) {
855 0 : FlavContainer flav1( idOldL);
856 0 : FlavContainer flav2( -idNewQ);
857 : int iTry = 0;
858 0 : int idNewL = flavSelPtr->combine( flav1, flav2);
859 0 : while (++iTry < NTRYMAX && idNewL == 0)
860 0 : idNewL = flavSelPtr->combine( flav1, flav2);
861 0 : if (idNewL == 0) {
862 0 : infoPtr->errorMsg("Error in RHadrons::produceSquark: "
863 : "cannot form light hadron code");
864 0 : return false;
865 : }
866 0 : double mNewL = particleDataPtr->mSel( idNewL);
867 :
868 : // Check that energy enough for light hadron and R-hadron.
869 0 : if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) {
870 0 : Vec4 pRHad, pNewL;
871 0 : if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) {
872 0 : infoPtr->errorMsg("Error in RHadrons::produceSquark: "
873 : "failed to construct kinematics for two-hadron decay");
874 0 : return false;
875 : }
876 :
877 : // Insert R-hadron and light hadron.
878 0 : iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
879 0 : pRHad, mRHad, 0.);
880 0 : event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0,
881 0 : pNewL, mNewL, 0.);
882 :
883 : // Done for two-body case.
884 : nBody = 2;
885 0 : }
886 0 : }
887 :
888 : // Final case: for very low mass collapse to one single R-hadron.
889 0 : if (nBody == 0) {
890 0 : idRHad = toIdWithSquark( idOldH, idOldL);
891 0 : if (idRHad == 0) {
892 0 : infoPtr->errorMsg("Error in RHadrons::produceSquark: "
893 : "cannot form R-hadron code");
894 0 : return false;
895 : }
896 :
897 : // Insert R-hadron with new momentum.
898 0 : iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
899 0 : systemPtr->pSum, systemPtr->mass, 0.);
900 :
901 : // Done with one-body case.
902 : nBody = 1;
903 0 : }
904 :
905 : // History bookkeeping: the R-hadron and processed partons.
906 0 : iRHadron[iRHad] = iRNow;
907 0 : int iLast = event.size() - 1;
908 0 : for (int i = iBeg; i <= iOldL; ++i) {
909 0 : event[i].statusNeg();
910 0 : event[i].daughters( iRNow, iLast);
911 : }
912 :
913 : // Remove old string system and, if needed, insert a new one.
914 0 : colConfig.erase(iSys);
915 0 : if (nBody == 3) {
916 0 : vector<int> iNewSys;
917 0 : iNewSys.push_back( iNewQ);
918 0 : iNewSys.push_back( iNewL);
919 0 : for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys.push_back( i);
920 0 : colConfig.insert( iNewSys, event);
921 0 : }
922 :
923 : // Copy lifetime and vertex from sparticle to R-hadron.
924 0 : event[iRNow].tau( event[iBef].tau() );
925 0 : 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 0 : }
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 0 : vector<int> iSide1, iSide2, iSideTmp, iNewSys1, iNewSys2;
944 0 : Vec4 pGlui, pSide1, pSide2;
945 :
946 : // Decide whether to produce a gluinoball or not.
947 0 : bool isGBall = (rndmPtr->flat() < probGluinoballRH);
948 :
949 : // Extract one string system on either side of the gluino.
950 0 : for (int i = 0; i < systemPtr->size(); ++i) {
951 0 : int iTmp = systemPtr->iParton[i];
952 0 : if (iGlui == 0 && event[ iTmp ].id() == idRGo) {
953 0 : iGlui = iTmp;
954 0 : pGlui = event[ iTmp ].p();
955 0 : } else if (iGlui == 0) {
956 0 : iSideTmp.push_back( iTmp);
957 0 : pSide1 += event[ iTmp ].p();
958 0 : } else {
959 0 : iSide2.push_back( iTmp);
960 0 : pSide2 += event[ iTmp ].p();
961 : }
962 0 : }
963 :
964 : // Order sides from gluino outwards and with lowest relative mass first.
965 0 : for (int i = int(iSideTmp.size()) - 1; i >= 0; --i)
966 0 : iSide1.push_back( iSideTmp[i]);
967 0 : double m1H = (pGlui + pSide1).mCalc() - event[ iSide1.back() ].m();
968 0 : double m2H = (pGlui + pSide2).mCalc() - event[ iSide2.back() ].m();
969 0 : if (m2H < m1H) {
970 0 : swap( iSide1, iSide2);
971 0 : swap( pSide1, pSide2);
972 : }
973 :
974 : // Begin loop over two sides of gluino, with lowest relative mass first.
975 0 : 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 0 : int iNewQ = 0;
983 0 : int iNewL = 0;
984 : int statusRHad = 0;
985 :
986 : // Copy down system consecutively from gluino end.
987 0 : int iBeg = event.size();
988 0 : if (iSide == 1) {
989 0 : iCreRHad[iRHad] = iBeg;
990 0 : event.copy( iGlui, 102);
991 0 : for (int i = 0; i < int(iSide1.size()); ++i)
992 0 : event.copy( iSide1[i], 102);
993 0 : } else {
994 0 : event.copy( iGlui, 102);
995 0 : for (int i = 0; i < int(iSide2.size()); ++i)
996 0 : event.copy( iSide2[i], 102);
997 : }
998 0 : int iEnd = event.size() - 1;
999 :
1000 : // Pick new flavour to help form R-hadron. Initial colour values.
1001 0 : int idOldL = event[iEnd].id();
1002 : int idNewQ = 21;
1003 0 : if (!isGBall) {
1004 : do {
1005 0 : FlavContainer flavOld( idOldL);
1006 0 : idNewQ = -flavSelPtr->pick(flavOld).id;
1007 0 : } while (iSide == 2 && isDiq1 && abs(idNewQ) > 10);
1008 0 : if (iSide == 1) isDiq1 = (abs(idNewQ) > 10);
1009 : }
1010 0 : bool hasCol = (event[iEnd].col() > 0);
1011 : int colR = 0;
1012 : int acolR = 0;
1013 :
1014 : // Target intermediary mass or R-hadron mass.
1015 0 : if (iSide == 1) {
1016 : idSave = idNewQ;
1017 0 : idRHad = (hasCol) ? 1009002 : -1009002;
1018 0 : if (hasCol) colR = event[iBeg].col();
1019 0 : else acolR = event[iBeg].acol();
1020 : statusRHad = 103;
1021 0 : double mNewQ = particleDataPtr->constituentMass( idNewQ);
1022 0 : if (isGBall) mNewQ *= 0.5;
1023 0 : mRHad = event[iBeg].m() + mOffsetCloudRH + mNewQ;
1024 0 : } else {
1025 0 : idRHad = toIdWithGluino( idSave, idNewQ);
1026 0 : if (idRHad == 0) {
1027 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1028 : "cannot form R-hadron code");
1029 0 : return false;
1030 : }
1031 : statusRHad = 104;
1032 0 : mRHad = particleDataPtr->m0(idRHad) + event[iBeg].m() - m0Go;
1033 : }
1034 :
1035 : // z value of fragmentation function.
1036 0 : double z = zSelPtr->zFrag( idRGo, idNewQ, mRHad*mRHad);
1037 :
1038 : // Basic kinematics of string piece where break is to occur.
1039 0 : Vec4 pOldH = event[iBeg].p();
1040 0 : int iOldL = iBeg + 1;
1041 0 : Vec4 pOldL = event[iOldL].p();
1042 0 : double mOldL = event[iOldL].m();
1043 0 : double mNewH = mRHad / z;
1044 0 : double sSys = (pOldH + pOldL).m2Calc();
1045 0 : double sRem = (1. - z) * (sSys - mNewH*mNewH);
1046 0 : double sMin = pow2(mOldL + mCollapseRH);
1047 :
1048 : // If too low remaining mass in system then add one more parton to it.
1049 0 : while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) )
1050 0 : && iOldL < iEnd ) {
1051 0 : ++iOldL;
1052 0 : pOldL += event[iOldL].p();
1053 0 : mOldL = event[iOldL].m();
1054 0 : sSys = (pOldH + pOldL).m2Calc();
1055 0 : sRem = (1. - z) * (sSys - mNewH*mNewH);
1056 0 : sMin = pow2(mOldL + mCollapseRH);
1057 : }
1058 :
1059 : // If enough mass then split off R-hadron and reduced system.
1060 0 : if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) {
1061 0 : Vec4 pNewH, pNewL;
1062 0 : if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) {
1063 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1064 : "failed to construct kinematics with reduced system");
1065 0 : return false;
1066 : }
1067 :
1068 : // Insert intermediary or R-hadron with new momentum, less colour.
1069 0 : iRNow = event.append( idRHad, statusRHad, iBeg, iOldL,
1070 0 : 0, 0, colR, acolR, z * pNewH, mRHad, 0.);
1071 :
1072 : // Reduced system with new string endpoint and modified recoiler.
1073 0 : if (!isGBall) idNewQ = -idNewQ;
1074 0 : int colN = (hasCol) ? 0 : event[iOldL].acol();
1075 0 : int acolN = (hasCol) ? event[iOldL].col() : 0;
1076 0 : iNewQ = event.append( idNewQ, 105, iBeg, iOldL, 0, 0,
1077 0 : colN, acolN, (1. - z) * pNewH, (1. - z) * mNewH, 0.);
1078 0 : iNewL = event.copy( iOldL, 105);
1079 0 : event[iNewL].mothers( iBeg, iOldL);
1080 0 : event[iNewL].p( pNewL);
1081 :
1082 : // Collect partons of new string system.
1083 0 : if (iSide == 1) {
1084 0 : iNewSys1.push_back( iNewQ);
1085 0 : iNewSys1.push_back( iNewL);
1086 0 : for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys1.push_back( i);
1087 0 : } else {
1088 0 : iNewSys2.push_back( iNewQ);
1089 0 : iNewSys2.push_back( iNewL);
1090 0 : 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 0 : }
1096 :
1097 : // For side-1 low-mass glueball system reabsorb full momentum.
1098 0 : if (nBody == 0 && isGBall && iSide == 1) {
1099 0 : idQLeap = event[iOldL].id();
1100 0 : Vec4 pNewH = event[iBeg].p() + pOldL;
1101 :
1102 : // Insert intermediary R-hadron with new momentum, less colour.
1103 0 : iRNow = event.append( idRHad, statusRHad, iBeg, iEnd,
1104 0 : 0, 0, colR, acolR, pNewH, pNewH.mCalc(), 0.);
1105 : nBody = 1;
1106 0 : }
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 0 : if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) {
1111 0 : if (isGBall) idNewQ = -idQLeap;
1112 0 : FlavContainer flav1( idOldL);
1113 0 : FlavContainer flav2( -idNewQ);
1114 : int iTry = 0;
1115 0 : int idNewL = flavSelPtr->combine( flav1, flav2);
1116 0 : while (++iTry < NTRYMAX && idNewL == 0)
1117 0 : idNewL = flavSelPtr->combine( flav1, flav2);
1118 0 : if (idNewL == 0) {
1119 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1120 : "cannot form light hadron code");
1121 0 : return false;
1122 : }
1123 0 : double mNewL = particleDataPtr->mSel( idNewL);
1124 :
1125 : // Check that energy enough for light hadron and R-hadron.
1126 0 : if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) {
1127 0 : Vec4 pRHad, pNewL;
1128 0 : if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) {
1129 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1130 : "failed to construct kinematics for two-hadron decay");
1131 0 : return false;
1132 : }
1133 :
1134 : // Insert intermediary or R-hadron and light hadron.
1135 0 : iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0,
1136 0 : colR, acolR, pRHad, mRHad, 0.);
1137 0 : event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0,
1138 0 : 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 0 : }
1145 0 : }
1146 :
1147 : // Final case: for very low mass collapse to one single R-hadron.
1148 0 : if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) {
1149 0 : if (isGBall) idSave = idQLeap;
1150 0 : if (iSide == 1) idSave = idOldL;
1151 0 : else idRHad = toIdWithGluino( idSave, idOldL);
1152 0 : if (idRHad == 0) {
1153 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1154 : "cannot form R-hadron code");
1155 0 : return false;
1156 : }
1157 :
1158 : // Insert R-hadron with new momentum.
1159 0 : iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0,
1160 0 : 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 0 : }
1166 :
1167 : // History bookkeeping: the processed partons.
1168 0 : int iLast = event.size() - 1;
1169 0 : for (int i = iBeg; i <= iOldL; ++i) {
1170 0 : event[i].statusNeg();
1171 0 : event[i].daughters( iRNow, iLast);
1172 : }
1173 :
1174 : // End of loop over two sides of the gluino.
1175 : iGlui = iRNow;
1176 0 : }
1177 :
1178 : // History bookkeeping: insert R-hadron; delete old string system.
1179 0 : if (iGlui == 0) {
1180 0 : infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1181 : "could not handle gluinoball kinematics");
1182 0 : return false;
1183 : }
1184 0 : iRHadron[iRHad] = iGlui;
1185 0 : colConfig.erase(iSys);
1186 :
1187 : // Non-gluinoball: insert (up to) two new string systems.
1188 0 : if (!isGBall) {
1189 0 : if (iNewSys1.size() > 0) colConfig.insert( iNewSys1, event);
1190 0 : 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 0 : } else if (idQLeap == 0) {
1195 0 : int iG1 = iNewSys1[0];
1196 0 : int iG2 = iNewSys2[0];
1197 0 : int colG = event[iG1].col() + event[iG2].col();
1198 0 : int acolG = event[iG1].acol() + event[iG2].acol();
1199 0 : Vec4 pG = event[iG1].p() + event[iG2].p();
1200 0 : int iG12 = event.append( 21, 107, iG1, iG2, 0, 0, colG, acolG,
1201 0 : pG, pG.mCalc(), 0.);
1202 :
1203 : // Temporary gluons no longer needed, but new colour to avoid warnings.
1204 0 : event[iG1].id( 21);
1205 0 : event[iG2].id( 21);
1206 0 : event[iG1].statusNeg();
1207 0 : event[iG2].statusNeg();
1208 0 : event[iG1].daughter1( iG12);
1209 0 : event[iG2].daughter1( iG12);
1210 0 : int colBridge = event.nextColTag();
1211 0 : if (event[iG1].col() == 0) {
1212 0 : event[iG1].col( colBridge);
1213 0 : event[iG2].acol( colBridge);
1214 0 : } else {
1215 0 : event[iG1].acol( colBridge);
1216 0 : event[iG2].col( colBridge);
1217 : }
1218 :
1219 : // Form the remnant system from which the R-hadron has been split off.
1220 0 : vector<int> iNewSys12;
1221 0 : for (int i = int(iNewSys1.size()) - 1; i > 0; --i)
1222 0 : iNewSys12.push_back( iNewSys1[i]);
1223 0 : iNewSys12.push_back( iG12);
1224 0 : for (int i = 1; i < int(iNewSys2.size()); ++i)
1225 0 : iNewSys12.push_back( iNewSys2[i]);
1226 0 : 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 0 : } else {
1231 0 : int iG2 = iNewSys2[0];
1232 0 : event[iG2].id( idQLeap);
1233 0 : colConfig.insert( iNewSys2, event);
1234 : }
1235 :
1236 : // Copy lifetime and vertex from sparticle to R-hadron.
1237 0 : event[iGlui].tau( event[iBef].tau() );
1238 0 : if (event[iBef].hasVertex()) event[iGlui].vProd( event[iBef].vProd() );
1239 :
1240 : // Done with production of a R-hadron from a gluino.
1241 0 : return true;
1242 :
1243 0 : }
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 0 : int id1Abs = abs(id1);
1254 0 : int id2Abs = abs(id2);
1255 0 : if (id2Abs < 10 && id1 > 0 && id2 > 0) return 0;
1256 0 : if (id2Abs < 10 && id1 < 0 && id2 < 0) return 0;
1257 0 : if (id2Abs > 10 && id1 > 0 && id2 < 0) return 0;
1258 0 : if (id2Abs > 10 && id1 < 0 && id2 > 0) return 0;
1259 :
1260 : // Form R-hadron code. Flip sign for antisquark.
1261 0 : bool isSt = (id1Abs == idRSt);
1262 : int idRHad = 1000000;
1263 0 : if (id2Abs < 10) idRHad += ((isSt) ? 600 : 500) + 10 * id2Abs + 2;
1264 0 : else idRHad += ((isSt) ? 6000 : 5000) + 10 * (id2Abs/100) + id2Abs%10;
1265 0 : if (id1 < 0) idRHad = -idRHad;
1266 :
1267 : // Done.
1268 : return idRHad;
1269 :
1270 0 : }
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 0 : int idLight = (abs(idRHad) - 1000000) / 10;
1281 0 : int idSq = (idLight < 100) ? idLight/10 : idLight/100;
1282 0 : int id1 = (idSq == 6) ? idRSt : idRSb;
1283 0 : if (idRHad < 0) id1 = -id1;
1284 :
1285 : // Find light (di)quark flavour content.
1286 0 : int id2 = (idLight < 100) ? idLight%10 : idLight%100;
1287 0 : if (id2 > 10) id2 = 100 * id2 + abs(idRHad)%10;
1288 0 : if ((id2 < 10 && idRHad > 0) || (id2 > 10 && idRHad < 0)) id2 = -id2;
1289 :
1290 : // Done.
1291 0 : return make_pair( id1, id2);
1292 :
1293 0 : }
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 0 : int id1Abs = abs(id1);
1303 0 : int id2Abs = abs(id2);
1304 0 : if (id1Abs == 21 && id2Abs == 21) return 1000993;
1305 0 : int idMax = max( id1Abs, id2Abs);
1306 0 : int idMin = min( id1Abs, id2Abs);
1307 0 : if (idMin > 10) return 0;
1308 0 : if (idMax > 10 && id1 > 0 && id2 < 0) return 0;
1309 0 : if (idMax > 10 && id1 < 0 && id2 > 0) return 0;
1310 0 : if (idMax < 10 && id1 > 0 && id2 > 0) return 0;
1311 0 : if (idMax < 10 && id1 < 0 && id2 < 0) return 0;
1312 :
1313 : // Form R-meson code. Flip sign for antiparticle.
1314 : int idRHad = 0;
1315 0 : if (idMax < 10) {
1316 0 : idRHad = 1009003 + 100 * idMax + 10 * idMin;
1317 0 : if (idMin != idMax && idMax%2 == 1) {
1318 0 : if (id1Abs == idMax && id1 > 0) idRHad = -idRHad;
1319 0 : if (id2Abs == idMax && id2 > 0) idRHad = -idRHad;
1320 : }
1321 0 : if (idMin != idMax && idMax%2 == 0) {
1322 0 : if (id1Abs == idMax && id1 < 0) idRHad = -idRHad;
1323 0 : if (id2Abs == idMax && id2 < 0) idRHad = -idRHad;
1324 : }
1325 :
1326 : // Form R-baryon code. Flip sign for antiparticle.
1327 : } else {
1328 0 : int idA = idMax/1000;
1329 0 : int idB = (idMax/100)%10;
1330 0 : int idC = idMin;
1331 0 : if (idC > idB) swap( idB, idC);
1332 0 : if (idB > idA) swap( idA, idB);
1333 0 : if (idC > idB) swap( idB, idC);
1334 0 : idRHad = 1090004 + 1000 * idA + 100 * idB + 10 * idC;
1335 0 : if (id1 < 0) idRHad = -idRHad;
1336 0 : }
1337 :
1338 : // Done.
1339 : return idRHad;
1340 :
1341 0 : }
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 0 : int idLight = (abs(idRHad) - 1000000) / 10;
1352 0 : int id1, id2, idTmp, idA, idB, idC;
1353 :
1354 : // Gluinoballs: split g into d dbar or u ubar.
1355 0 : if (idLight < 100) {
1356 0 : id1 = (rndmPtr->flat() < 0.5) ? 1 : 2;
1357 0 : id2 = -id1;
1358 :
1359 : // Gluino-meson: split into q + qbar.
1360 0 : } else if (idLight < 1000) {
1361 0 : id1 = (idLight / 10) % 10;
1362 0 : id2 = -(idLight % 10);
1363 : // Flip signs when first quark of down-type.
1364 0 : if (id1%2 == 1) {
1365 : idTmp = id1;
1366 0 : id1 = -id2;
1367 0 : id2 = -idTmp;
1368 0 : }
1369 :
1370 : // Gluino-baryon: split to q + qq (diquark).
1371 : // Pick diquark at random, except if c or b involved.
1372 : } else {
1373 0 : idA = (idLight / 100) % 10;
1374 0 : idB = (idLight / 10) % 10;
1375 0 : idC = idLight % 10;
1376 0 : double rndmQ = 3. * rndmPtr->flat();
1377 0 : if (idA > 3) rndmQ = 0.5;
1378 0 : if (rndmQ < 1.) {
1379 0 : id1 = idA;
1380 0 : id2 = 1000 * idB + 100 * idC + 3;
1381 0 : if (idB != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1382 0 : } else if (rndmQ < 2.) {
1383 0 : id1 = idB;
1384 0 : id2 = 1000 * idA + 100 * idC + 3;
1385 0 : if (idA != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1386 : } else {
1387 0 : id1 = idC;
1388 0 : id2 = 1000 * idA + 100 * idB +3;
1389 0 : if (idA != idB && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1390 : }
1391 : }
1392 :
1393 : // Flip signs for anti-R-hadron.
1394 0 : if (idRHad < 0) {
1395 0 : idTmp = id1;
1396 0 : id1 = -id2;
1397 0 : id2 = -idTmp;
1398 0 : }
1399 :
1400 : // Done.
1401 0 : return make_pair( id1, id2);
1402 :
1403 0 : }
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 0 : double sSum = (pOld1 + pOld2).m2Calc();
1416 0 : double sOld1 = pOld1.m2Calc();
1417 0 : double sOld2 = pOld2.m2Calc();
1418 0 : double sNew1 = mNew1 * mNew1;
1419 0 : double sNew2 = mNew2 * mNew2;
1420 :
1421 : // Check that kinematically possible.
1422 0 : if (checkMargin && pow2(mNew1 + mNew2 + MSAFETY) > sSum) return false;
1423 :
1424 : // Transfer coefficients to give four-vectors with the new masses.
1425 0 : double lamOld = sqrt( pow2(sSum - sOld1 - sOld2) - 4. * sOld1 * sOld2 );
1426 0 : double lamNew = sqrt( pow2(sSum - sNew1 - sNew2) - 4. * sNew1 * sNew2 );
1427 0 : double move1 = (lamNew * (sSum - sOld1 + sOld2)
1428 0 : - lamOld * (sSum - sNew1 + sNew2)) / (2. * sSum * lamOld);
1429 0 : double move2 = (lamNew * (sSum + sOld1 - sOld2)
1430 0 : - lamOld * (sSum + sNew1 - sNew2)) / (2. * sSum * lamOld);
1431 :
1432 : // Construct final vectors. Done.
1433 0 : pNew1 = (1. + move1) * pOld1 - move2 * pOld2;
1434 0 : pNew2 = (1. + move2) * pOld2 - move1 * pOld1;
1435 : return true;
1436 :
1437 0 : }
1438 :
1439 : //==========================================================================
1440 :
1441 : } // end namespace Pythia8
|
