PYTHIA8/pythia8130/src/SigmaSUSY.cxx

   1 // SigmaSUSY.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the
   7 // supersymmetry simulation classes.
   8
   9 #include "SigmaSUSY.h"
  10
  11 namespace Pythia8 {
  12
  13 //**************************************************************************
  14
  15 // Sigma2qqbar2gauginogaugino "mother" class.
  16 // Cross section for gaugino pair production: neutralino pair,
  17 // neutralino-chargino, and chargino pair production all inherit from this.
  18
  19 //*********
  20
  21 // Initialize process.
  22
  23 void Sigma2qqbar2gauginogaugino::initProc() {
  24
  25   // Construct name of process.
  26   nameSave = "q qbar -> " + ParticleDataTable::name(id3) + " "
  27     + ParticleDataTable::name(id4);
  28
  29   // Count number of final-state charged particles (default = 0)
  30   // (useful since charginos inherit from this class.)
  31   nCharged=0;
  32   if (abs(id3) == 1000024 || abs(id3) == 1000037) nCharged++;
  33   if (abs(id4) == 1000024 || abs(id4) == 1000037) nCharged++;
  34
  35   // Set up couplings
  36   mZpole    = ParticleDataTable::m0(23);
  37   wZpole    = ParticleDataTable::mWidth(23);
  38   double mWpole = ParticleDataTable::m0(24);
  39
  40   // Running masses and weak mixing angle
  41   // (default to pole values if no running available)
  42   double mW = mWpole;
  43   double mZ = mZpole;
  44   sin2W     = 1.0 - pow(mW/mZ,2);
  45   if (slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
  46       && slhaPtr->hmix.exists(3)) {
  47     double gp=slhaPtr->gauge(1);
  48     double g =slhaPtr->gauge(2);
  49     double v =slhaPtr->hmix(3);
  50     mW      = g * v / 2.0;
  51     mZ      = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
  52     sin2W   = pow2(gp)/(pow2(g)+pow2(gp));
  53   }
  54
  55   // Shorthand for SUSY couplings
  56   // By default, use the running one in HMIX
  57   // If not found, use the MZ one in MINPAR
  58   double tanb = slhaPtr->hmix.exists(2) ?
  59     slhaPtr->hmix(2) : slhaPtr->minpar(3);
  60   double sinW=sqrt(sin2W);
  61   double cosW=sqrt(1.0-sin2W);
  62   double cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
  63   double sinb = sqrt(max(0.0,1.0-cosb*cosb));
  64   SusyLesHouches::matrixblock<6> Ru(slhaPtr->usqmix);
  65   SusyLesHouches::matrixblock<6> Rd(slhaPtr->dsqmix);
  66   SusyLesHouches::matrixblock<6> imRu(slhaPtr->imusqmix);
  67   SusyLesHouches::matrixblock<6> imRd(slhaPtr->imusqmix);
  68
  69   // Local complex copies of neutralino mixing matrix entries.
  70   // (incl NMSSM for future use)
  71   complex ni1,nj1,ni2,nj2,ni3,nj3,ni4,nj4,ni5,nj5;
  72   if (slhaPtr->modsel(3) != 1) {
  73     ni1=complex( slhaPtr->nmix(id3chi,1), slhaPtr->imnmix(id3chi,1) );
  74     nj1=complex( slhaPtr->nmix(id4chi,1), slhaPtr->imnmix(id4chi,1) );
  75     ni2=complex( slhaPtr->nmix(id3chi,2), slhaPtr->imnmix(id3chi,2) );
  76     nj2=complex( slhaPtr->nmix(id4chi,2), slhaPtr->imnmix(id4chi,2) );
  77     ni3=complex( slhaPtr->nmix(id3chi,3), slhaPtr->imnmix(id3chi,3) );
  78     nj3=complex( slhaPtr->nmix(id4chi,3), slhaPtr->imnmix(id4chi,3) );
  79     ni4=complex( slhaPtr->nmix(id3chi,4), slhaPtr->imnmix(id3chi,4) );
  80     nj4=complex( slhaPtr->nmix(id4chi,4), slhaPtr->imnmix(id4chi,4) );
  81     ni5=complex( 0.,0.);
  82     nj5=complex( 0.,0.);
  83   } else {
  84     ni1=complex( slhaPtr->nmnmix(id3chi,1), slhaPtr->imnmnmix(id3chi,1) );
  85     nj1=complex( slhaPtr->nmnmix(id4chi,1), slhaPtr->imnmnmix(id4chi,1) );
  86     ni2=complex( slhaPtr->nmnmix(id3chi,2), slhaPtr->imnmnmix(id3chi,2) );
  87     nj2=complex( slhaPtr->nmnmix(id4chi,2), slhaPtr->imnmnmix(id4chi,2) );
  88     ni3=complex( slhaPtr->nmnmix(id3chi,3), slhaPtr->imnmnmix(id3chi,3) );
  89     nj3=complex( slhaPtr->nmnmix(id4chi,3), slhaPtr->imnmnmix(id4chi,3) );
  90     ni4=complex( slhaPtr->nmnmix(id3chi,4), slhaPtr->imnmnmix(id3chi,4) );
  91     nj4=complex( slhaPtr->nmnmix(id4chi,4), slhaPtr->imnmnmix(id4chi,4) );
  92     ni5=complex( slhaPtr->nmnmix(id3chi,5), slhaPtr->imnmnmix(id3chi,5) );
  93     nj5=complex( slhaPtr->nmnmix(id4chi,5), slhaPtr->imnmnmix(id4chi,5) );
  94   }
  95
  96   // Change to positive mass convention.
  97   complex iRot( 0., 1.);
  98   if (slhaPtr->mass(abs(id3)) < 0.) {
  99     ni1 *= iRot;
 100     ni2 *= iRot;
 101     ni3 *= iRot;
 102     ni4 *= iRot;
 103     ni5 *= iRot;
 104   };
 105   if (slhaPtr->mass(abs(id4)) < 0.) {
 106     nj1 *= iRot;
 107     nj2 *= iRot;
 108     nj3 *= iRot;
 109     nj4 *= iRot;
 110     nj5 *= iRot;
 111   };
 112
 113   // Local copies of Chargino mixing
 114   complex ui1,ui2,vi1,vi2,uj1,uj2,vj1,vj2;
 115   if (id3chi <= 2) {
 116     ui1=complex( slhaPtr->umix(id3chi,1), slhaPtr->imumix(id3chi,1) );
 117     ui2=complex( slhaPtr->umix(id3chi,2), slhaPtr->imumix(id3chi,2) );
 118     vi1=complex( slhaPtr->vmix(id3chi,1), slhaPtr->imvmix(id3chi,1) );
 119     vi2=complex( slhaPtr->vmix(id3chi,2), slhaPtr->imvmix(id3chi,2) );
 120   }
 121   if (id4chi <= 2) {
 122     uj1=complex( slhaPtr->umix(id4chi,1), slhaPtr->imumix(id4chi,1) );
 123     uj2=complex( slhaPtr->umix(id4chi,2), slhaPtr->imumix(id4chi,2) );
 124     vj1=complex( slhaPtr->vmix(id4chi,1), slhaPtr->imvmix(id4chi,1) );
 125     vj2=complex( slhaPtr->vmix(id4chi,2), slhaPtr->imvmix(id4chi,2) );
 126   }
 127
 128   // Z chi_i chi_j
 129   OLpp = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
 130   ORpp =  0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
 131
 132   // Z cha_i cha_j
 133   OLp = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
 134     + ( (id3chi == id4chi) ? sin2W : 0.0);
 135   ORp = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
 136     + ( (id3chi == id4chi) ? sin2W : 0.0);
 137
 138   // W chi_i cha_j
 139   OL = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
 140   OR = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
 141
 142   // Z q_{idq} q_{idq} (def with extra factor 2 compared to [Okun])
 143   for (int idq = 1; idq <= 5; ++idq) {
 144     // No FCNC in Zqq, so here sum over diagonal only
 145     // No t in beam, so only sum up to 5
 146     LqqZ[idq] = CoupEW::lf(idq);
 147     RqqZ[idq] = CoupEW::rf(idq);
 148   }
 149
 150   // ~chi^0_i ~q_jsq idq
 151   for (int jsq = 1; jsq<=6; jsq++) {
 152     // Sum jsq over all squarks
 153     for (int idq = 1; idq<=5; idq++) {
 154       // No t in beam, so only sum iq over 5 lightest quarks.
 155
 156       // quark index
 157       int k=(idq+1)/2;
 158
 159       // Set quark mass for use in couplings
 160       // Initial guess 0,0,0,mc,mb, with the latter from the PDT
 161       double mq = ParticleDataTable::m0(idq);
 162       if (idq <= 3) mq=0.0;
 163
 164       // Treat u and d quarks separately
 165       if (idq % 2 == 1) {
 166
 167         // idq = d quark
 168         double eq = -1.0/3.0;
 169         double T3q = -0.5;
 170
 171         // Compute running mass from Yukawas and vevs if possible.
 172         if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
 173           double ykk=slhaPtr->yd(k,k);
 174           double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
 175           if (ykk != 0.0) mq = ykk * v1 / sqrt(2.0) ;
 176           //      cout <<scientific<<" q "<<idq << "  (k="<<k<<")  Y = "
 177           //     <<ykk<<"  v = "<<v1<<"  m = "<<mq<<endl;
 178         }
 179
 180         // Shorthand for squark mixing matrices
 181         complex RdjkL(Rd(jsq,k),imRd(jsq,k));
 182         complex RdjkR(Rd(jsq,k+3),imRd(jsq,k+3));
 183
 184         // L(~d_jsq,d_k,~chi^0_i)
 185         LsqXi[jsq][idq]=((eq-T3q)*sinW*ni1+T3q*cosW*ni2) * conj(RdjkL)
 186           + mq*cosW*ni3*conj(RdjkR)/2.0/mW/cosb ;
 187
 188         // L(~d_jsq,d_k,~chi^0_j)
 189         LsqXj[jsq][idq]=((eq-T3q)*sinW*nj1+T3q*cosW*nj2) * conj(RdjkL)
 190           + mq*cosW*nj3*conj(RdjkR)/2.0/mW/cosb ;
 191
 192         // R(~d_jsq,d_k,~chi^0_i)
 193         RsqXi[jsq][idq]=eq*sinW*ni1*RdjkR - mq*cosW*ni3*RdjkL/2.0/mW/cosb ;
 194         RsqXi[jsq][idq]=-conj(RsqXi[jsq][idq]) ;
 195
 196         // R(~d_jsq,d_k,~chi^0_j)
 197         RsqXj[jsq][idq]=eq*sinW*nj1*RdjkR - mq*cosW*nj3*RdjkL/2.0/mW/cosb ;
 198         RsqXj[jsq][idq]=-conj(RsqXj[jsq][idq]) ;
 199
 200         // Chargino couplings (initialize)
 201         LsqCi[jsq][idq]=0.0;
 202         RsqCi[jsq][idq]=0.0;
 203         LsqCj[jsq][idq]=0.0;
 204         RsqCj[jsq][idq]=0.0;
 205
 206         // Sum over flavours (off-diagonal mixings)
 207         // (note l <-> k with respect to [Klasen et al])
 208         for (int l=1; l<=3; l++) {
 209
 210           // Set quark masses for use in couplings
 211           // Initial guess 0,0,0,mc,mb, with the latter from the PDT
 212           double mu = ParticleDataTable::m0(2*l);
 213           if (2*l <= 3) mu=0.0;
 214
 215           // Compute running u mass from Yukawas and vevs if possible.
 216           if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
 217             double yll=slhaPtr->yu(l,l);
 218             double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
 219             if (yll != 0.0) mu = yll * v2 / sqrt(2.0) ;
 220           }
 221
 222           // Shorthand for u squark mixing
 223           complex RujlL(Ru(jsq,l),imRu(jsq,l));
 224           complex RujlR(Ru(jsq,l+3),imRu(jsq,l+3));
 225
 226           // CKM matrix (use Pythia-8 one if no SLHA)
 227           complex Vlk=VCKM::Vid(l,k);
 228           if (slhaPtr->vckm.exists())
 229             Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
 230
 231           // L(~u_jsq,d_k,~chi^+/-_i)
 232           LsqCi[jsq][idq] += Vlk*(conj(vi1)*RujlL
 233                                   -mu*conj(vi2)*RujlR/sqrt(2.0)/mW/sinb );
 234
 235           // L(~u_jsq,d_k,~chi^+/-_j)
 236           LsqCj[jsq][idq] += Vlk*(conj(vj1)*RujlL
 237                                   -mu*conj(vj2)*RujlR/sqrt(2.0)/mW/sinb );
 238
 239           // R(~u_jsq,d_k,~chi^+/-_i)
 240           RsqCi[jsq][idq] += - Vlk*mq*ui2*RujlL/sqrt(2.0)/mW/cosb ;
 241
 242           // R(~u_jsq,d_k,~chi^+/-_j)
 243           RsqCj[jsq][idq] += - Vlk*mq*uj2*RujlL/sqrt(2.0)/mW/cosb ;
 244
 245         }
 246
 247       } else {
 248
 249         // idq = u quark
 250         double eq = 2.0/3.0;
 251         double T3q = 0.5;
 252
 253         // Compute mass from Yukawas and vevs if possible.
 254         if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
 255           double ykk=slhaPtr->yu(k,k);
 256           double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
 257           if (ykk != 0.0) mq = ykk * v2 / sqrt(2.0) ;
 258           //      cout <<scientific<<" q "<<idq << "  (k="<<k<<")  Y = "
 259           //     <<ykk<<"  v = "<<v2<<"  m = "<<mq<<endl;
 260         }
 261
 262         // Shorthand for squark mixing matrices
 263         complex RujkL(Ru(jsq,k),imRu(jsq,k));
 264         complex RujkR(Ru(jsq,k+3),imRu(jsq,k+3));
 265
 266         // L(~u_jsq,u_k,~chi^0_i)
 267         LsqXi[jsq][idq]=((eq-T3q)*sinW*ni1+T3q*cosW*ni2) * conj(RujkL)
 268           + mq*cosW*ni4*conj(RujkR)/2.0/mW/sinb ;
 269
 270         // L(~u_jsq,u_k,~chi^0_j)
 271         LsqXj[jsq][idq]=((eq-T3q)*sinW*nj1+T3q*cosW*nj2) * conj(RujkL)
 272           + mq*cosW*nj4*conj(RujkR)/2.0/mW/sinb ;
 273
 274         // R(~u_jsq,u_k,~chi^0_i)
 275         RsqXi[jsq][idq]=eq*sinW*ni1*RujkR - mq*cosW*ni4*RujkL/2.0/mW/sinb ;
 276         RsqXi[jsq][idq]=-conj(RsqXi[jsq][idq]);
 277
 278         // R(~u_jsq,u_k,~chi^0_i)
 279         RsqXj[jsq][idq]=eq*sinW*nj1*RujkR - mq*cosW*nj4*RujkL/2.0/mW/sinb ;
 280         RsqXj[jsq][idq]=-conj(RsqXj[jsq][idq]);
 281
 282         // Chargino couplings (initialize)
 283         LsqCi[jsq][idq]=0.0;
 284         RsqCi[jsq][idq]=0.0;
 285         LsqCj[jsq][idq]=0.0;
 286         RsqCj[jsq][idq]=0.0;
 287
 288         // Sum over flavours (off-diagonal mixings)
 289         // (note l <-> k with respect to [Klasen et al])
 290         for (int l=1; l<=3; l++) {
 291
 292           // Set quark masses for use in couplings
 293           // Initial guess 0,0,0,mc,mb, with the latter from the PDT
 294           double md = ParticleDataTable::m0(2*l-1);
 295           if (2*l-1 <= 3) md=0.0;
 296
 297           // Compute running d mass from Yukawas and vevs if possible.
 298           if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
 299             double yll=slhaPtr->yd(l,l);
 300             double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
 301             if (yll != 0.0) md = yll * v1 / sqrt(2.0) ;
 302           }
 303
 304           // Shorthand for d squark mixing
 305           complex RdjlL(Rd(jsq,l),imRd(jsq,l));
 306           complex RdjlR(Rd(jsq,l+3),imRd(jsq,l+3));
 307
 308           // CKM matrix (use Pythia-8 one if no SLHA)
 309           complex Vkl=VCKM::Vid(k,l);
 310           if (slhaPtr->vckm.exists())
 311             Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
 312
 313           // L(~d_jsq,u_k,~chi^+/-_i)
 314           LsqCi[jsq][idq] += Vkl*(ui1*conj(RdjlL)
 315                                   -md*ui2*conj(RdjlR)/sqrt(2.0)/mW/cosb) ;
 316
 317           // L(~d_jsq,u_k,~chi^+/-_j)
 318           LsqCj[jsq][idq] += Vkl*(uj1*conj(RdjlL)
 319                                   -md*uj2*conj(RdjlR)/sqrt(2.0)/mW/cosb) ;
 320
 321           // R(~d_jsq,u_k,~chi^+/-_i)
 322           RsqCi[jsq][idq] += - Vkl*mq*conj(vi2*RdjlL)/sqrt(2.0)/mW/sinb ;
 323
 324           // R(~d_jsq,u_k,~chi^+/-_j)
 325           RsqCj[jsq][idq] += - Vkl*mq*conj(vi2*RdjlL)/sqrt(2.0)/mW/sinb ;
 326
 327         }
 328
 329       }
 330     }
 331   }
 332 }
 333
 334 //*********
 335
 336 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 337
 338 void Sigma2qqbar2gauginogaugino::sigmaKin() {
 339
 340   // Common flavour-independent factor.
 341   sigma0 = (M_PI / sH2) * pow2(alpEM) ;
 342
 343   // Factor 1/2 for identical final particles.
 344   if (id3 == id4 && nCharged == 0) sigma0 *= 0.5;
 345
 346   // Auxiliary factors for use below
 347   ui     = uH - s3;
 348   uj     = uH - s4;
 349   ti     = tH - s3;
 350   tj     = tH - s4;
 351   sz     = sH - pow2(mZpole);
 352   d      = pow2(sz) + pow2(mZpole * wZpole);
 353   propZ  = complex( sz / d, mZpole * wZpole / d);
 354
 355 }
 356
 357 //*********
 358
 359 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 360
 361 // neutralino pairs:
 362
 363 double Sigma2qqbar2gauginogaugino::sigmaHat() {
 364
 365   // Only allow quark-antiquark incoming states
 366   if (id1*id2 >= 0) {
 367     return 0.0;
 368   }
 369
 370   // Only allow incoming states with sum(charge) = 0
 371   if ((id1+id2) % 2 != 0) {
 372     return 0.0;
 373   }
 374
 375   // Flavour-dependent kinematics-dependent couplings.
 376   int idAbs1    = abs(id1);
 377   int idAbs2    = abs(id2);
 378   complex QuLL = (id1 == -id2) ? LqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
 379   complex QtLL = (id1 == -id2) ? LqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
 380   complex QuRR = (id1 == -id2) ? RqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
 381   complex QtRR = (id1 == -id2) ? RqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
 382   complex QuLR = 0.0;
 383   complex QtLR = 0.0;
 384   complex QuRL = 0.0;
 385   complex QtRL = 0.0;
 386
 387   // Add t-channel squark flavour sums to QmXY couplings
 388   for (int jsq=1; jsq<=6; jsq++) {
 389     int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1+1) % 2 + 1;
 390     double msq2=pow(ParticleDataTable::m0(idsq),2);
 391     double usq     = uH - msq2;
 392     double tsq     = tH - msq2;
 393     QuLL += LsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
 394     QtLL -= conj(LsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
 395     QuRR += RsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
 396     QtRR -= conj(RsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
 397     QuLR += RsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
 398     QtLR += conj(LsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
 399     QuRL += LsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
 400     QtRL += conj(RsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
 401   }
 402
 403   // Multiply by overall normalization
 404   QuLL *= 1.0/(sin2W * (1 - sin2W));
 405   QtLL *= 1.0/(sin2W * (1 - sin2W));
 406   QuRR *= 1.0/(sin2W * (1 - sin2W));
 407   QtRR *= 1.0/(sin2W * (1 - sin2W));
 408   QuLR *= 1.0/(sin2W * (1 - sin2W));
 409   QtLR *= 1.0/(sin2W * (1 - sin2W));
 410   QuRL *= 1.0/(sin2W * (1 - sin2W));
 411   QtRL *= 1.0/(sin2W * (1 - sin2W));
 412
 413   // Compute matrix element weight
 414   double weight = 0;
 415   // Average over separate helicity contributions
 416   // LL (ha = -1, hb = +1) (divided by 4 for average)
 417   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 418     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
 419   // RR (ha =  1, hb = -1) (divided by 4 for average)
 420   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 421     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
 422   // RL (ha =  1, hb =  1) (divided by 4 for average)
 423   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 424     - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
 425   // LR (ha = -1, hb = -1) (divided by 4 for average)
 426   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 427     - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
 428
 429   // Cross section, including colour factor.
 430   double sigma = sigma0 * weight;
 431   if (idAbs1 < 9) sigma /= 3.;
 432
 433   // Answer.
 434   return sigma;
 435
 436 }
 437
 438 // neutralino-chargino
 439
 440 double Sigma2qqbar2chi0char::sigmaHat() {
 441
 442   // This cross section not fully debugged yet, skip.
 443   return 0.0;
 444
 445   // Only allow quark-antiquark incoming states
 446   if (id1*id2 >= 0) {
 447     return 0.0;
 448   }
 449
 450   // Only allow incoming states with sum(charge) = +/- 1
 451   if ((id1+id2) % 2 == 0) {
 452     return 0.0;
 453   }
 454
 455   // Flavour-dependent kinematics-dependent couplings.
 456   int idAbs1    = abs(id1);
 457   int idAbs2    = abs(id2);
 458   complex QuLL = (id1 == -id2) ? LqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
 459   complex QtLL = (id1 == -id2) ? LqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
 460   complex QuRR = (id1 == -id2) ? RqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
 461   complex QtRR = (id1 == -id2) ? RqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
 462   complex QuLR = 0.0;
 463   complex QtLR = 0.0;
 464   complex QuRL = 0.0;
 465   complex QtRL = 0.0;
 466
 467   // Add t-channel squark flavour sums to QmXY couplings
 468   for (int jsq=1; jsq<=6; jsq++) {
 469     int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1+1) % 2 + 1;
 470     double msq2=pow(ParticleDataTable::m0(idsq),2);
 471     double usq     = uH - msq2;
 472     double tsq     = tH - msq2;
 473     QuLL += LsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
 474     QtLL -= conj(LsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
 475     QuRR += RsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
 476     QtRR -= conj(RsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
 477     QuLR += RsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
 478     QtLR += conj(LsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
 479     QuRL += LsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
 480     QtRL += conj(RsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
 481   }
 482
 483   // Compute matrix element weight
 484   double weight = 0;
 485   // Average over separate helicity contributions
 486   // LL (ha = -1, hb = +1) (divided by 4 for average)
 487   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 488     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
 489   // RR (ha =  1, hb = -1) (divided by 4 for average)
 490   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 491     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
 492   // RL (ha =  1, hb =  1) (divided by 4 for average)
 493   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 494     - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
 495   // LR (ha = -1, hb = -1) (divided by 4 for average)
 496   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 497     - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
 498
 499   // Cross section, including colour factor.
 500   double sigma = sigma0 * weight;
 501   if (idAbs1 < 9) sigma /= 3.;
 502
 503   // Answer.
 504   return sigma;
 505
 506 }
 507
 508 // chargino-chargino
 509
 510 double Sigma2qqbar2charchar::sigmaHat() {
 511
 512   // This cross section not fully debugged yet, skip.
 513   return 0.0;
 514
 515   // Only allow quark-antiquark incoming states
 516   if (id1*id2 >= 0) {
 517     return 0.0;
 518   }
 519
 520   // Only allow incoming states with sum(charge) = 0
 521   if ((id1+id2) % 2 != 0) {
 522     return 0.0;
 523   }
 524
 525   // Flavour-dependent kinematics-dependent couplings.
 526   int idAbs1    = abs(id1);
 527   int idAbs2    = abs(id2);
 528
 529   // Add s-channel Z
 530   complex QuLL = (id1 == -id2) ? - LqqZ[idAbs1] * conj(OLp)/2.0 * propZ
 531     / (sin2W * (1 - sin2W)) : 0.0;
 532   complex QtLL = (id1 == -id2) ? - LqqZ[idAbs1] * conj(ORp)/2.0 * propZ
 533     / (sin2W * (1 - sin2W)) : 0.0;
 534   complex QuRR = (id1 == -id2) ? - RqqZ[idAbs1] * conj(ORp)/2.0 * propZ
 535     / (sin2W * (1 - sin2W)) : 0.0;
 536   complex QtRR = (id1 == -id2) ? - RqqZ[idAbs1] * conj(OLp)/2.0 * propZ
 537     / (sin2W * (1 - sin2W)) : 0.0;
 538   complex QuLR = 0.0;
 539   complex QtLR = 0.0;
 540   complex QuRL = 0.0;
 541   complex QtRL = 0.0;
 542
 543   // Add s-channel photon
 544   if (idAbs1 == idAbs2 && abs(id3) == abs(id4)) {
 545     double e = -1.0/3.0;
 546     if (idAbs1 % 2 == 0) e = 2.0/3.0;
 547     if (idAbs1 > 10) e = -1.0;
 548     QuLL += e/sH;
 549     QtLL += e/sH;
 550     QuRR += e/sH;
 551     QtRR += e/sH;
 552   }
 553
 554   // Add t-channel squark flavour sums to QmXY couplings
 555   for (int jsq=1; jsq<=6; jsq++) {
 556     int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1 % 2) + 1;
 557     double msq2=pow(ParticleDataTable::m0(idsq),2);
 558     double usq     = uH - msq2;
 559     double tsq     = tH - msq2;
 560
 561     // up-type quarks get u-channel d-squark contributions
 562     if (idAbs1 % 2 == 0 ) {
 563
 564       // Flip 1<->2 and i<->j if id1 is antiquark
 565       if (id1 > 0) {
 566         QuLL += LsqCi[jsq][idAbs2]*conj(LsqCj[jsq][idAbs1])/usq/2.0/sin2W;
 567         QuRR += RsqCi[jsq][idAbs2]*conj(RsqCj[jsq][idAbs1])/usq/2.0/sin2W;
 568         QuLR += RsqCi[jsq][idAbs2]*conj(LsqCj[jsq][idAbs1])/usq/2.0/sin2W;
 569         QuRL += LsqCi[jsq][idAbs2]*conj(RsqCj[jsq][idAbs1])/usq/2.0/sin2W;
 570       } else {
 571         QuLL += LsqCj[jsq][idAbs1]*conj(LsqCi[jsq][idAbs2])/usq/2.0/sin2W;
 572         QuRR += RsqCj[jsq][idAbs1]*conj(RsqCi[jsq][idAbs2])/usq/2.0/sin2W;
 573         QuLR += RsqCj[jsq][idAbs1]*conj(LsqCi[jsq][idAbs2])/usq/2.0/sin2W;
 574         QuRL += LsqCj[jsq][idAbs1]*conj(RsqCi[jsq][idAbs2])/usq/2.0/sin2W;
 575       }
 576     }
 577
 578     // down-type quarks get t-channel u-squark contributions
 579     else if (idAbs1 % 2 == 1) {
 580       // Flip 1<->2 and i<->j if id1 is antiquark
 581       if (id1 > 0) {
 582         QtLL -= LsqCi[jsq][idAbs1]*conj(LsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
 583         QtRR -= RsqCi[jsq][idAbs1]*conj(RsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
 584         QtLR += LsqCi[jsq][idAbs1]*conj(RsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
 585         QtRL += RsqCi[jsq][idAbs1]*conj(LsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
 586       } else {
 587         QtLL -= LsqCj[jsq][idAbs2]*conj(LsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
 588         QtRR -= RsqCj[jsq][idAbs2]*conj(RsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
 589         QtLR += LsqCj[jsq][idAbs2]*conj(RsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
 590         QtRL += RsqCj[jsq][idAbs2]*conj(LsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
 591       }
 592
 593     }
 594   }
 595
 596   // Compute matrix element weight
 597   double weight = 0;
 598   // Average over separate helicity contributions
 599   // LL (ha = -1, hb = +1) (divided by 4 for average)
 600   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 601     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
 602   // RR (ha =  1, hb = -1) (divided by 4 for average)
 603   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 604     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
 605   // RL (ha =  1, hb =  1) (divided by 4 for average)
 606   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 607     - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
 608   // LR (ha = -1, hb = -1) (divided by 4 for average)
 609   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 610     - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
 611
 612   // Cross section, including colour factor.
 613   double sigma = sigma0 * weight;
 614   if (idAbs1 < 9) sigma /= 3.;
 615
 616   // Answer.
 617   return sigma;
 618
 619 }
 620
 621 //*********
 622
 623 // Select identity, colour and anticolour.
 624
 625 void Sigma2qqbar2gauginogaugino::setIdColAcol() {
 626
 627   // Set flavours.
 628
 629   if (nCharged == 0) {
 630     // Neutralino-Neutralino
 631     setId( id1, id2, id3, id4);
 632
 633   } else if (nCharged == 1) {
 634     // Neutralino-Chargino: set charge according to in-states
 635     if ( abs(id1%1) == 1 && id1 > 0 ) {
 636       setId( id1, id2, id3, -id4);
 637     } else {
 638       setId( id1, id2, id3, id4);
 639     }
 640
 641   } else if (nCharged == 2) {
 642     // Chargino-Chargino: charge already included in id3, id4
 643     setId( id1, id2, id3, id4);
 644   }
 645
 646   // Colour flow topologies. Swap when antiquarks.
 647   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
 648   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 649   if (id1 < 0) swapColAcol();
 650
 651 }
 652
 653 //**************************************************************************
 654
 655 } // end namespace Pythia8
 656