1 // SigmaCompositeness.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 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.
6 // Function definitions (not found in the header) for the
7 // compositeness simulation classes.
9 #include "SigmaCompositeness.h"
13 //==========================================================================
15 // Sigma1qg2qStar class.
16 // Cross section for q g -> q^* (excited quark state).
17 // Note: for simplicity decay is assumed isotropic.
19 //--------------------------------------------------------------------------
21 // Initialize process.
23 void Sigma1qg2qStar::initProc() {
25 // Set up process properties from the chosen quark flavour.
26 idRes = 4000000 + idq;
27 codeSave = 4000 + idq;
28 if (idq == 1) nameSave = "d g -> d^*";
29 else if (idq == 2) nameSave = "u g -> u^*";
30 else if (idq == 3) nameSave = "s g -> s^*";
31 else if (idq == 4) nameSave = "c g -> c^*";
32 else nameSave = "b g -> b^*";
34 // Store q* mass and width for propagator.
35 mRes = particleDataPtr->m0(idRes);
36 GammaRes = particleDataPtr->mWidth(idRes);
38 GamMRat = GammaRes / mRes;
40 // Locally stored properties and couplings.
41 Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
42 coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol");
44 // Set pointer to particle properties and decay table.
45 qStarPtr = particleDataPtr->particleDataEntryPtr(idRes);
49 //--------------------------------------------------------------------------
51 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
53 void Sigma1qg2qStar::sigmaKin() {
55 // Incoming width for correct quark.
56 widthIn = pow3(mH) * alpS * pow2(coupFcol) / (3. * pow2(Lambda));
58 // Set up Breit-Wigner.
59 sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
63 //--------------------------------------------------------------------------
65 // Evaluate sigmaHat(sHat) for specific incoming flavours.
67 double Sigma1qg2qStar::sigmaHat() {
69 // Identify whether correct incoming flavours.
70 int idqNow = (id2 == 21) ? id1 : id2;
71 if (abs(idqNow) != idq) return 0.;
73 // Outgoing width and total sigma. Done.
74 return widthIn * sigBW * qStarPtr->resWidthOpen(idqNow, mH);
78 //--------------------------------------------------------------------------
80 // Select identity, colour and anticolour.
82 void Sigma1qg2qStar::setIdColAcol() {
85 int idqNow = (id2 == 21) ? id1 : id2;
86 int idqStar = (idqNow > 0) ? idRes : -idRes;
87 setId( id1, id2, idqStar);
89 // Colour flow topology.
90 if (id1 == idqNow) setColAcol( 1, 0, 2, 1, 2, 0);
91 else setColAcol( 2, 1, 1, 0, 2, 0);
92 if (idqNow < 0) swapColAcol();
96 //--------------------------------------------------------------------------
98 // Evaluate weight for q* decay angle.
100 double Sigma1qg2qStar::weightDecay( Event& process, int iResBeg,
103 // q* should sit in entry 5. Sequential Z/W decay assumed isotropic.
104 if (iResBeg != 5 || iResEnd != 5) return 1.;
106 // Sign of asymmetry.
107 int sideIn = (process[3].idAbs() < 20) ? 1 : 2;
108 int sideOut = (process[6].idAbs() < 20) ? 1 : 2;
109 double eps = (sideIn == sideOut) ? 1. : -1.;
111 // Phase space factors.
112 double mr1 = pow2(process[6].m()) / sH;
113 double mr2 = pow2(process[7].m()) / sH;
114 double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
116 // Reconstruct decay angle. Default isotropic decay.
117 double cosThe = (process[3].p() - process[4].p())
118 * (process[7].p() - process[6].p()) / (sH * betaf);
122 // Decay q* -> q (g/gamma) or q (Z^0/W^+-).
123 int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
124 if (idBoson == 21 || idBoson == 22) {
125 wt = 1. + eps * cosThe;
127 } else if (idBoson == 23 || idBoson == 24) {
128 double mrB = (sideOut == 1) ? mr2 : mr1;
129 double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
130 wt = 1. + eps * cosThe * ratB;
139 //==========================================================================
141 // Sigma1lgm2lStar class.
142 // Cross section for l gamma -> l^* (excited lepton state).
143 // Note: for simplicity decay is assumed isotropic.
145 //--------------------------------------------------------------------------
147 // Initialize process.
149 void Sigma1lgm2lStar::initProc() {
151 // Set up process properties from the chosen lepton flavour.
152 idRes = 4000000 + idl;
153 codeSave = 4000 + idl;
154 if (idl == 11) nameSave = "e gamma -> e^*";
155 else if (idl == 13) nameSave = "mu gamma -> mu^*";
156 else nameSave = "tau gamma -> tau^*";
158 // Store l* mass and width for propagator.
159 mRes = particleDataPtr->m0(idRes);
160 GammaRes = particleDataPtr->mWidth(idRes);
162 GamMRat = GammaRes / mRes;
164 // Locally stored properties and couplings.
165 Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
166 double coupF = settingsPtr->parm("ExcitedFermion:coupF");
167 double coupFp = settingsPtr->parm("ExcitedFermion:coupFprime");
168 coupChg = -0.5 * coupF - 0.5 * coupFp;
170 // Set pointer to particle properties and decay table.
171 qStarPtr = particleDataPtr->particleDataEntryPtr(idRes);
175 //--------------------------------------------------------------------------
177 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
179 void Sigma1lgm2lStar::sigmaKin() {
181 // Incoming width for correct lepton.
182 widthIn = pow3(mH) * alpEM * pow2(coupChg) / pow2(Lambda);
184 // Set up Breit-Wigner.
185 sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
189 //--------------------------------------------------------------------------
191 // Evaluate sigmaHat(sHat) for specific incoming flavours.
193 double Sigma1lgm2lStar::sigmaHat() {
195 // Identify whether correct incoming flavours.
196 int idlNow = (id2 == 22) ? id1 : id2;
197 if (abs(idlNow) != idl) return 0.;
199 // Outgoing width and total sigma. Done.
200 return widthIn * sigBW * qStarPtr->resWidthOpen(idlNow, mH);
204 //--------------------------------------------------------------------------
206 // Select identity, colour and anticolour.
208 void Sigma1lgm2lStar::setIdColAcol() {
211 int idlNow = (id2 == 22) ? id1 : id2;
212 int idlStar = (idlNow > 0) ? idRes : -idRes;
213 setId( id1, id2, idlStar);
216 setColAcol( 0, 0, 0, 0, 0, 0);
220 //--------------------------------------------------------------------------
222 // Evaluate weight for l* decay angle.
224 double Sigma1lgm2lStar::weightDecay( Event& process, int iResBeg,
227 // l* should sit in entry 5. Sequential Z/W decay assumed isotropic.
228 if (iResBeg != 5 || iResEnd != 5) return 1.;
230 // Sign of asymmetry.
231 int sideIn = (process[3].idAbs() < 20) ? 1 : 2;
232 int sideOut = (process[6].idAbs() < 20) ? 1 : 2;
233 double eps = (sideIn == sideOut) ? 1. : -1.;
235 // Phase space factors.
236 double mr1 = pow2(process[6].m()) / sH;
237 double mr2 = pow2(process[7].m()) / sH;
238 double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
240 // Reconstruct decay angle. Default isotropic decay.
241 double cosThe = (process[3].p() - process[4].p())
242 * (process[7].p() - process[6].p()) / (sH * betaf);
246 // Decay l* -> l gamma or l (Z^0/W^+-).
247 int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
249 wt = 1. + eps * cosThe;
251 } else if (idBoson == 23 || idBoson == 24) {
252 double mrB = (sideOut == 1) ? mr2 : mr1;
253 double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
254 wt = 1. + eps * cosThe * ratB;
263 //==========================================================================
265 // Sigma2qq2qStarq class.
266 // Cross section for q q' -> q^* q' (excited quark state).
267 // Note: for simplicity decay is assumed isotropic.
269 //--------------------------------------------------------------------------
271 // Initialize process.
273 void Sigma2qq2qStarq::initProc() {
275 // Set up process properties from the chosen quark flavour.
276 idRes = 4000000 + idq;
277 codeSave = 4020 + idq;
278 if (idq == 1) nameSave = "q q -> d^* q";
279 else if (idq == 2) nameSave = "q q -> u^* q";
280 else if (idq == 3) nameSave = "q q -> s^* q";
281 else if (idq == 4) nameSave = "q q -> c^* q";
282 else nameSave = "q q -> b^* q";
284 // Locally stored properties and couplings.
285 Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
286 preFac = M_PI / pow4(Lambda);
288 // Secondary open width fractions.
289 openFracPos = particleDataPtr->resOpenFrac( idRes);
290 openFracNeg = particleDataPtr->resOpenFrac(-idRes);
294 //--------------------------------------------------------------------------
296 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
298 void Sigma2qq2qStarq::sigmaKin() {
300 // Two possible expressions, for like or unlike sign.
301 sigmaA = preFac * (1. - s3 / sH);
302 sigmaB = preFac * (-uH) * (sH + tH) / sH2;
306 //--------------------------------------------------------------------------
308 // Evaluate sigmaHat(sHat) for specific incoming flavours.
310 double Sigma2qq2qStarq::sigmaHat() {
312 // Identify different allowed incoming flavour combinations.
313 int id1Abs = abs(id1);
314 int id2Abs = abs(id2);
315 double open1 = (id1 > 0) ? openFracPos : openFracNeg;
316 double open2 = (id2 > 0) ? openFracPos : openFracNeg;
319 if (id1Abs == idq) sigma += (4./3.) * sigmaA * open1;
320 if (id2Abs == idq) sigma += (4./3.) * sigmaA * open2;
321 } else if (id1Abs == idq && id2 == -id1)
322 sigma = (8./3.) * sigmaB * (open1 + open2);
323 else if (id2 == -id1) sigma = sigmaB * (open1 + open2);
324 else if (id1Abs == idq) sigma = sigmaB * open1;
325 else if (id2Abs == idq) sigma = sigmaB * open2;
332 //--------------------------------------------------------------------------
334 // Select identity, colour and anticolour.
336 void Sigma2qq2qStarq::setIdColAcol() {
338 // Flavours: either side may have been excited.
341 if (abs(id1) == idq) open1 = (id1 > 0) ? openFracPos : openFracNeg;
342 if (abs(id2) == idq) open2 = (id2 > 0) ? openFracPos : openFracNeg;
343 if (open1 == 0. && open2 == 0.) {
344 open1 = (id1 > 0) ? openFracPos : openFracNeg;
345 open2 = (id2 > 0) ? openFracPos : openFracNeg;
347 bool excite1 = (open1 > 0.);
348 if (open1 > 0. && open2 > 0.) excite1
349 = (rndmPtr->flat() * (open1 + open2) < open1);
351 // Always excited quark in slot 3 so colour flow flipped or not.
353 id3 = (id1 > 0) ? idq : -idq;
355 if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
356 else setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
357 if (id1 < 0) swapColAcol();
359 id3 = (id2 > 0) ? idq : -idq;
362 if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
363 else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
364 if (id1 < 0) swapColAcol();
366 setId( id1, id2, id3, id4);
370 //==========================================================================
372 // Sigma2qqbar2lStarlbar class.
373 // Cross section for q qbar -> l^* lbar (excited lepton state).
374 // Note: for simplicity decay is assumed isotropic.
376 //--------------------------------------------------------------------------
378 // Initialize process.
380 void Sigma2qqbar2lStarlbar::initProc() {
382 // Set up process properties from the chosen lepton flavour.
383 idRes = 4000000 + idl;
384 codeSave = 4020 + idl;
385 if (idl == 11) nameSave = "q qbar -> e^*+- e^-+";
386 else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_ebar";
387 else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^-+";
388 else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mubar";
389 else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^-+";
390 else nameSave = "q qbar -> nu_tau^* nu_taubar";
392 // Secondary open width fractions.
393 openFracPos = particleDataPtr->resOpenFrac( idRes);
394 openFracNeg = particleDataPtr->resOpenFrac(-idRes);
396 // Locally stored properties and couplings.
397 Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
398 preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 3.;
402 //--------------------------------------------------------------------------
404 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
406 void Sigma2qqbar2lStarlbar::sigmaKin() {
408 // Only one possible expressions
409 sigma = preFac * (-uH) * (sH + tH) / sH2;
413 //--------------------------------------------------------------------------
415 // Select identity, colour and anticolour.
417 void Sigma2qqbar2lStarlbar::setIdColAcol() {
419 // Flavours: either lepton or antilepton may be excited.
420 if (rndmPtr->flat() * (openFracPos + openFracNeg) < openFracPos) {
421 setId( id1, id2, idRes, -idl);
422 if (id1 < 0) swapTU = true;
424 setId( id1, id2, -idRes, idl);
425 if (id1 > 0) swapTU = true;
428 // Colour flow trivial.
429 if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
430 else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
434 //==========================================================================
436 // Sigma2QCqq2qq class.
437 // Cross section for q q -> q q (quark contact interactions).
439 //--------------------------------------------------------------------------
441 // Initialize process.
443 void Sigma2QCqq2qq::initProc() {
445 m_Lambda2 = settingsPtr->parm("ContactInteractions:Lambda");
446 m_etaLL = settingsPtr->mode("ContactInteractions:etaLL");;
447 m_etaRR = settingsPtr->mode("ContactInteractions:etaRR");;
448 m_etaLR = settingsPtr->mode("ContactInteractions:etaLR");;
449 m_Lambda2 *= m_Lambda2;
453 //--------------------------------------------------------------------------
455 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
457 void Sigma2QCqq2qq::sigmaKin() {
459 // Calculate kinematics dependence for different terms.
460 sigT = (4./9.) * (sH2 + uH2) / tH2;
461 sigU = (4./9.) * (sH2 + tH2) / uH2;
462 sigTU = - (8./27.) * sH2 / (tH * uH);
463 sigST = - (8./27.) * uH2 / (sH * tH);
465 sigQCSTU = sH2 * (1 / tH + 1 / uH);
466 sigQCUTS = uH2 * (1 / tH + 1 / sH);
470 //--------------------------------------------------------------------------
472 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
474 double Sigma2QCqq2qq::sigmaHat() {
476 // Terms from QC contact interactions.
481 // Combine cross section terms; factor 1/2 when identical quarks.
485 sigSum = 0.5 * (sigT + sigU + sigTU); // SM terms.
488 sigQCLL = (8./9.) * alpS * (m_etaLL/m_Lambda2) * sigQCSTU
489 + (8./3.) * pow2(m_etaLL/m_Lambda2) * sH2;
490 sigQCRR = (8./9.) * alpS * (m_etaRR/m_Lambda2) * sigQCSTU
491 + (8./3.) * pow2(m_etaRR/m_Lambda2) * sH2;
492 sigQCLR = 2. * (uH2 + tH2) * pow2(m_etaLR/m_Lambda2);
498 // q qbar -> q qbar, without pure s-channel term.
499 } else if (id2 == -id1) {
501 sigSum = sigT + sigST; // SM terms.
503 // Contact terms, minus the terms included in qqbar2qqbar.
504 sigQCLL = (8./9.) * alpS * (m_etaLL/m_Lambda2) * sigQCUTS
505 + (5./3.) * pow2(m_etaLL/m_Lambda2) * uH2;
506 sigQCRR = (8./9.) * alpS * (m_etaRR/m_Lambda2) * sigQCUTS
507 + (5./3.) * pow2(m_etaRR/m_Lambda2) * uH2;
508 sigQCLR = 2. * sH2 * pow2(m_etaLR/m_Lambda2);
510 // q q' -> q q' or q qbar' -> q qbar'
513 sigSum = sigT; // SM terms.
517 sigQCLL = pow2(m_etaLL/m_Lambda2) * sH2;
518 sigQCRR = pow2(m_etaRR/m_Lambda2) * sH2;
519 sigQCLR = 2 * pow2(m_etaLR/m_Lambda2) * uH2;
521 sigQCLL = pow2(m_etaLL/m_Lambda2) * uH2;
522 sigQCRR = pow2(m_etaRR/m_Lambda2) * uH2;
523 sigQCLR = 2 * pow2(m_etaLR/m_Lambda2) * sH2;
528 double sigma = (M_PI/sH2) * ( pow2(alpS) * sigSum
529 + sigQCLL + sigQCRR + sigQCLR );
534 //--------------------------------------------------------------------------
536 // Select identity, colour and anticolour.
538 void Sigma2QCqq2qq::setIdColAcol() {
540 // Outgoing = incoming flavours.
541 setId( id1, id2, id1, id2);
543 // Colour flow topologies. Swap when antiquarks.
544 if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
545 else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
546 if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
547 setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
548 if (id1 < 0) swapColAcol();
552 //==========================================================================
554 // Sigma2QCqqbar2qqbar class.
555 // Cross section for q qbar -> q' qbar' (quark contact interactions).
557 //--------------------------------------------------------------------------
559 // Initialize process.
561 void Sigma2QCqqbar2qqbar::initProc() {
563 m_nQuarkNew = settingsPtr->mode("ContactInteractions:nQuarkNew");
564 m_Lambda2 = settingsPtr->parm("ContactInteractions:Lambda");
565 m_etaLL = settingsPtr->mode("ContactInteractions:etaLL");
566 m_etaRR = settingsPtr->mode("ContactInteractions:etaRR");
567 m_etaLR = settingsPtr->mode("ContactInteractions:etaLR");
568 m_Lambda2 *= m_Lambda2;
572 //--------------------------------------------------------------------------
574 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
576 void Sigma2QCqqbar2qqbar::sigmaKin() {
579 idNew = 1 + int( m_nQuarkNew * rndmPtr->flat() );
580 mNew = particleDataPtr->m0(idNew);
583 // Calculate kinematics dependence.
586 if (sH > 4. * m2New) {
587 sigS = (4./9.) * (tH2 + uH2) / sH2;
588 sigQC = pow2(m_etaLL/m_Lambda2) * uH2
589 + pow2(m_etaRR/m_Lambda2) * uH2
590 + 2 * pow2(m_etaLR/m_Lambda2) * tH2;
593 // Answer is proportional to number of outgoing flavours.
594 sigma = (M_PI / sH2) * m_nQuarkNew * ( pow2(alpS) * sigS + sigQC);
598 //--------------------------------------------------------------------------
600 // Select identity, colour and anticolour.
602 void Sigma2QCqqbar2qqbar::setIdColAcol() {
604 // Set outgoing flavours ones.
605 id3 = (id1 > 0) ? idNew : -idNew;
606 setId( id1, id2, id3, -id3);
608 // Colour flow topologies. Swap when antiquarks.
609 setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
610 if (id1 < 0) swapColAcol();
614 //==========================================================================
616 } // end namespace Pythia8