Line data Source code
1 : // SusyCouplings.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // Main authors of this file: N. Desai, P. Skands
4 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 :
7 : // Function definitions (not found in the header) for the
8 : // supersymmetric couplings class.
9 :
10 : #include "Pythia8/ParticleData.h"
11 : #include "Pythia8/SusyCouplings.h"
12 :
13 : namespace Pythia8 {
14 :
15 : //==========================================================================
16 :
17 : // The CoupSUSY class.
18 :
19 : //--------------------------------------------------------------------------
20 :
21 : // Constants: could be changed here if desired, but normally should not.
22 : // These are of technical nature, as described for each.
23 :
24 : // Allow verbose printout for debug purposes.
25 : const bool CoupSUSY::DBSUSY = false;
26 :
27 : //--------------------------------------------------------------------------
28 :
29 : // Initialize SM+SUSY couplings (only performed once).
30 :
31 : void CoupSUSY::initSUSY (SusyLesHouches* slhaPtrIn, Info* infoPtrIn,
32 : ParticleData* particleDataPtrIn, Settings* settingsPtrIn) {
33 :
34 : // Save pointers.
35 0 : infoPtr = infoPtrIn;
36 0 : slhaPtr = slhaPtrIn;
37 0 : settingsPtr = settingsPtrIn;
38 0 : particleDataPtr = particleDataPtrIn;
39 :
40 : // Only initialize SUSY parts if SUSY is actually switched on
41 0 : if (!slhaPtr->modsel.exists()) return;
42 :
43 : // Is NMSSM switched on?
44 0 : isNMSSM = (slhaPtr->modsel(3) != 1 ? false : true);
45 0 : settingsPtr->flag("SLHA:NMSSM",isNMSSM);
46 0 : int nNeut = (isNMSSM ? 5 : 4);
47 : int nChar = 2;
48 :
49 : // Initialize pole masses
50 0 : mZpole = particleDataPtr->m0(23);
51 0 : wZpole = particleDataPtr->mWidth(23);
52 0 : mWpole = particleDataPtr->m0(24);
53 0 : wWpole = particleDataPtr->mWidth(24);
54 :
55 : // Running masses and weak mixing angle
56 : // (default to pole values if no running available)
57 0 : mW = mWpole;
58 0 : mZ = mZpole;
59 0 : sin2W = sin2thetaW();
60 :
61 0 : if (settingsPtr->mode("SUSY:sin2thetaWMode") == 3) {
62 : // Possibility to force on-shell sin2W definition, mostly intended for
63 : // cross-checks
64 0 : sin2W = 1.0 - pow(mW/mZ,2);
65 0 : } else if (settingsPtr->mode("SUSY:sin2thetaWMode") == 2){
66 : // Possibility to use running sin2W definition, derived from gauge
67 : // couplings in running SLHA blocks (normally defined at SUSY scale)
68 0 : if(slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
69 0 : && slhaPtr->hmix.exists(3)) {
70 0 : double gp=slhaPtr->gauge(1);
71 0 : double g =slhaPtr->gauge(2);
72 0 : double v =slhaPtr->hmix(3);
73 0 : mW = g * v / 2.0;
74 0 : mZ = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
75 : //double tan2W = pow2(gp)/pow2(g);
76 : //if (DBSUSY) cout << " tan2W = " << tan2W << endl;
77 0 : sin2W = pow2(gp)/(pow2(g)+pow2(gp));
78 0 : } else {
79 0 : infoPtr->errorMsg("Warning in CoupSUSY::initSUSY: Block GAUGE"
80 : " not found or incomplete; using sin(thetaW) at mZ");
81 : }
82 : }
83 :
84 : // Overload the SM value of sin2thetaW
85 0 : s2tW = sin2W;
86 :
87 0 : sinW = sqrt(sin2W);
88 0 : cosW = sqrt(1.0-sin2W);
89 :
90 : // Tan(beta)
91 : // By default, use the running one in HMIX (if not found, use MINPAR)
92 :
93 0 : if(slhaPtr->hmix.exists(2))
94 0 : tanb = slhaPtr->hmix(2);
95 : else{
96 0 : tanb = slhaPtr->minpar(3);
97 0 : infoPtr->errorMsg("Warning in CoupSUSY::initSUSY: Block HMIX"
98 : " not found or incomplete; using MINPAR tan(beta)");
99 : }
100 0 : cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
101 0 : sinb = sqrt( max(0.0, 1.0 - cosb*cosb));
102 :
103 0 : double beta = acos(cosb);
104 :
105 : // Verbose output
106 : if (DBSUSY) {
107 : cout << " sin2W(Q) = " << sin2W << " mW(Q) = " << mW
108 : << " mZ(Q) = " << mZ << endl;
109 : cout << " vev(Q) = " << slhaPtr->hmix(3) << " tanb(Q) = " << tanb
110 : << endl;
111 : for (int i=1;i<=3;i++) {
112 : for (int j=1;j<=3;j++) {
113 : cout << " VCKM [" << i << "][" << j << "] = "
114 : << scientific << setw(10) << VCKMgen(i,j) << endl;
115 : }
116 : }
117 : }
118 :
119 : // Higgs sector
120 0 : if(slhaPtr->alpha.exists()) {
121 0 : alphaHiggs = slhaPtr->alpha();
122 0 : }
123 : // If RPV, assume alpha = asin(RVHMIX(1,2)) (ignore Higgs-Sneutrino mixing)
124 0 : else if (slhaPtr->modsel(4) == 1) {
125 0 : alphaHiggs = asin(slhaPtr->rvhmix(1,2));
126 0 : infoPtr->errorMsg("Info from CoupSUSY::initSUSY:","Extracting angle"
127 : " alpha from RVHMIX", true);
128 0 : } else {
129 0 : infoPtr->errorMsg("Info from CoupSUSY::initSUSY:","Block ALPHA"
130 : " not found; using alpha = beta.", true);
131 : // Define approximate alpha by simple SM limit
132 0 : alphaHiggs = atan(tanb);
133 : }
134 :
135 0 : if(slhaPtr->hmix.exists(1) && slhaPtr->hmix.exists(4)){
136 0 : muHiggs = slhaPtr->hmix(1);
137 0 : mAHiggs = sqrt(slhaPtr->hmix(4));
138 0 : } else if (slhaPtr->rvamix.exists()){
139 0 : mAHiggs = particleDataPtr->m0(36);
140 0 : muHiggs = 0.0;
141 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: Block HMIX"
142 0 : " not found or incomplete","; setting mu = 0.");
143 0 : } else{
144 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: Block HMIX"
145 0 : " not found or incomplete","; setting mu = mA = 0.");
146 0 : muHiggs = 0.0;
147 0 : mAHiggs = 0.0;
148 : }
149 :
150 : // Pass SLHA input to 2HDM sector
151 :
152 0 : double sba = sin(beta-alphaHiggs);
153 0 : double cba = cos(beta-alphaHiggs);
154 0 : double cosa = cos(alphaHiggs);
155 0 : double sina = sin(alphaHiggs);
156 :
157 : // h0
158 0 : settingsPtr->parm("HiggsH1:coup2d", -sina/cosb);
159 0 : settingsPtr->parm("HiggsH1:coup2u", cosa/sinb);
160 0 : settingsPtr->parm("HiggsH1:coup2l", cosa/sinb);
161 0 : settingsPtr->parm("HiggsH1:coup2Z", sba);
162 0 : settingsPtr->parm("HiggsH1:coup2W", sba);
163 : // H0
164 0 : settingsPtr->parm("HiggsH2:coup2d", cosa/cosb);
165 0 : settingsPtr->parm("HiggsH2:coup2u", sina/sinb);
166 0 : settingsPtr->parm("HiggsH2:coup2l", sina/sinb);
167 0 : settingsPtr->parm("HiggsH2:coup2Z", cba);
168 0 : settingsPtr->parm("HiggsH2:coup2W", cba);
169 0 : settingsPtr->parm("HiggsH2:coup2H1Z", 0.0);
170 0 : settingsPtr->parm("HiggsH2:coup2HchgW", sba);
171 : // A0
172 0 : settingsPtr->parm("HiggsA3:coup2d", tanb);
173 0 : settingsPtr->parm("HiggsA3:coup2u", cosb/sinb);
174 0 : settingsPtr->parm("HiggsA3:coup2l", cosb/sinb);
175 0 : settingsPtr->parm("HiggsA3:coup2Z", 0.0);
176 0 : settingsPtr->parm("HiggsA3:coup2W", 0.0);
177 0 : settingsPtr->parm("HiggsA3:coup2H1Z", cba);
178 0 : settingsPtr->parm("HiggsA3:coup2H2Z", sba);
179 0 : settingsPtr->parm("HiggsA3:coup2HchgW", 1.0);
180 :
181 : // H^+
182 0 : settingsPtr->parm("HiggsHchg:tanBeta", tanb);
183 0 : settingsPtr->parm("HiggsHchg:coup2H1W", cba);
184 0 : settingsPtr->parm("HiggsHchg:coup2H2W", sba);
185 :
186 : // Triple higgs couplings
187 :
188 0 : double cbpa = cos(beta+alphaHiggs);
189 0 : double sbpa = sin(beta+alphaHiggs);
190 :
191 0 : settingsPtr->parm("HiggsH1:coup2Hchg", cos(2*beta)*sbpa + 2*pow2(cosW)*sba);
192 0 : settingsPtr->parm("HiggsH2:coup2Hchg", -cos(2*beta)*cbpa + 2*pow2(cosW)*cba);
193 0 : settingsPtr->parm("HiggsH2:coup2H1H1", 2*sin(2*alphaHiggs)*sbpa
194 0 : - cos(2*alphaHiggs)*cbpa);
195 0 : settingsPtr->parm("HiggsH2:coup2A3A3", -cos(2*beta)*cbpa);
196 0 : settingsPtr->parm("HiggsH2:coup2A3H1", 0.0);
197 0 : settingsPtr->parm("HiggsA3:coup2H1H1", 0.0);
198 0 : settingsPtr->parm("HiggsA3:coup2Hchg", 0.0);
199 :
200 : // Shorthand for squark mixing matrices
201 0 : LHmatrixBlock<6> Ru(slhaPtr->usqmix);
202 0 : LHmatrixBlock<6> Rd(slhaPtr->dsqmix);
203 0 : LHmatrixBlock<6> imRu(slhaPtr->imusqmix);
204 0 : LHmatrixBlock<6> imRd(slhaPtr->imdsqmix);
205 :
206 : // Construct ~g couplings
207 0 : for (int i=1; i<=6; i++) {
208 0 : for (int j=1; j<=3; j++) {
209 0 : LsddG[i][j] = complex( Rd(i,j) , imRd(i,j));
210 0 : RsddG[i][j] = complex(-Rd(i,j+3), -imRd(i,j+3));
211 0 : LsuuG[i][j] = complex( Ru(i,j) , imRu(i,j));
212 0 : RsuuG[i][j] = complex(-Ru(i,j+3), -imRu(i,j+3));
213 :
214 : if (DBSUSY) {
215 : cout << " Lsddg [" << i << "][" << j << "] = "
216 : << scientific << setw(10) << LsddG[i][j]
217 : << " RsddG [" << i << "][" << j << "] = "
218 : << scientific << setw(10) << RsddG[i][j] << endl;
219 : cout << " Lsuug [" << i << "][" << j << "] = "
220 : << scientific << setw(10) << LsuuG[i][j]
221 : << " RsuuG [" << i << "][" << j << "] = "
222 : << scientific << setw(10) << RsuuG[i][j] << endl;
223 : }
224 : }
225 : }
226 :
227 : // Construct qqZ couplings
228 0 : for (int i=1 ; i<=6 ; i++) {
229 :
230 : // q[i] q[i] Z (def with extra factor 2 compared to [Okun])
231 0 : LqqZ[i] = af(i) - 2.0*ef(i)*sin2W ;
232 0 : RqqZ[i] = - 2.0*ef(i)*sin2W ;
233 :
234 : // tmp: verbose output
235 : if (DBSUSY) {
236 : cout << " LqqZ [" << i << "][" << i << "] = "
237 : << scientific << setw(10) << LqqZ[i]
238 : << " RqqZ [" << i << "][" << i << "] = "
239 : << scientific << setw(10) << RqqZ[i] << endl;
240 : }
241 : }
242 :
243 : // Construct ~q~qZ couplings
244 0 : for (int i=1 ; i<=6 ; i++) {
245 :
246 : // Squarks can have off-diagonal couplings as well
247 0 : for (int j=1 ; j<=6 ; j++) {
248 :
249 : // ~d[i] ~d[j] Z
250 0 : LsdsdZ[i][j] = 0.0;
251 0 : RsdsdZ[i][j] = 0.0;
252 0 : for (int k=1;k<=3;k++) {
253 0 : complex Rdik = complex(Rd(i,k), imRd(i,k) );
254 0 : complex Rdjk = complex(Rd(j,k), imRd(j,k) );
255 0 : complex Rdik3 = complex(Rd(i,k+3),imRd(i,k+3));
256 0 : complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
257 0 : LsdsdZ[i][j] += LqqZ[1] * (Rdik*conj(Rdjk));
258 0 : RsdsdZ[i][j] += RqqZ[1] * (Rdik3*conj(Rdjk3));
259 0 : }
260 :
261 : // ~u[i] ~u[j] Z
262 0 : LsusuZ[i][j] = 0.0;
263 0 : RsusuZ[i][j] = 0.0;
264 0 : for (int k=1;k<=3;k++) {
265 0 : complex Ruik = complex(Ru(i,k) ,imRu(i,k) );
266 0 : complex Rujk = complex(Ru(j,k) ,imRu(j,k) );
267 0 : complex Ruik3 = complex(Ru(i,k+3),imRu(i,k+3));
268 0 : complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
269 0 : LsusuZ[i][j] += LqqZ[2] * (Ruik*conj(Rujk));
270 0 : RsusuZ[i][j] += RqqZ[2] * (Ruik3*conj(Rujk3));
271 0 : }
272 :
273 : // tmp: verbose output
274 : if (DBSUSY) {
275 : if (max(abs(LsdsdZ[i][j]),abs(RsdsdZ[i][j])) > 1e-6) {
276 : cout << " LsdsdZ[" << i << "][" << j << "] = "
277 : << scientific << setw(10) << LsdsdZ[i][j]
278 : << " RsdsdZ[" << i << "][" << j << "] = "
279 : << scientific << setw(10) << RsdsdZ[i][j] << endl;
280 : }
281 : if (max(abs(LsusuZ[i][j]),abs(RsusuZ[i][j]))> 1e-6) {
282 : cout << " LsusuZ[" << i << "][" << j << "] = "
283 : << scientific << setw(10) << LsusuZ[i][j]
284 : << " RsusuZ[" << i << "][" << j << "] = "
285 : << scientific << setw(10) << RsusuZ[i][j] << endl;
286 : }
287 : }
288 : }
289 : }
290 :
291 0 : for(int i = 1; i < 7; i++)
292 0 : for(int j = 1; j < 7; j++){
293 0 : Rsl[i][j] = slhaPtr->selmix(i,j);
294 : }
295 :
296 :
297 0 : for(int i = 1; i < 7; i++)
298 0 : for(int j = 1; j < 7; j++){
299 0 : if (i < 4 && j < 4) Rsv[i][j] = slhaPtr->snumix(i,j);
300 0 : else Rsv[i][j] = 0.0;
301 : }
302 :
303 : // In RPV, the slepton mixing matrices include Higgs bosons
304 : // Here we just extract the entries corresponding to the slepton PDG codes
305 : // I.e., slepton-Higgs mixing is not supported.
306 :
307 : // Charged sleptons
308 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvlmix.exists()) {
309 0 : for (int i=1; i<=6; ++i)
310 0 : for (int j=1; j<=6; ++j)
311 0 : Rsl[i][j] = slhaPtr->rvlmix(i+1,j+2);
312 : // Check for Higgs-slepton mixing in RVLMIX
313 : bool hasCrossTerms = false;
314 0 : for (int i=2; i<=7; ++i)
315 0 : for (int j=1; j<=2; ++j)
316 0 : if (abs(slhaPtr->rvlmix(i,j)) > 1e-5) {
317 : hasCrossTerms = true;
318 0 : break;
319 : }
320 0 : if (hasCrossTerms)
321 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
322 : "slepton-Higgs mixing not supported internally in PYTHIA");
323 0 : }
324 :
325 : // Neutral sleptons
326 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvhmix.exists()) {
327 0 : for (int i=1; i<=3; ++i)
328 0 : for (int j=1; j<=3; ++j)
329 0 : Rsv[i][j] = slhaPtr->rvhmix(i+2,j+2);
330 : // Check for Higgs-sneutrino mixing in RVHMIX
331 : bool hasCrossTerms = false;
332 0 : for (int i=3; i<=7; ++i)
333 0 : for (int j=1; j<=2; ++j)
334 0 : if (abs(slhaPtr->rvhmix(i,j)) > 1e-5) {
335 : hasCrossTerms = true;
336 0 : break;
337 : }
338 0 : if (hasCrossTerms)
339 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
340 : "sneutrino-Higgs mixing not supported internally in PYTHIA");
341 0 : }
342 :
343 : if(DBSUSY){
344 : cout<<"Rsl"<<endl;
345 : for(int i=1;i<=6;i++){
346 : for(int j=1;j<=6;j++){
347 : cout << scientific << setw(10) << Rsl[i][j]<<" ";
348 : }
349 : cout<<endl;
350 : }
351 : cout<<"Rsv"<<endl;
352 : for(int i=1;i<=6;i++){
353 : for(int j=1;j<=6;j++){
354 : cout << scientific << setw(10) << Rsv[i][j]<<" ";
355 : }
356 : cout<<endl;
357 : }
358 : }
359 :
360 :
361 : // Construct llZ couplings;
362 0 : for (int i=11 ; i<=16 ; i++) {
363 :
364 0 : LllZ[i-10] = af(i) - 2.0*ef(i)*sin2W ;
365 0 : RllZ[i-10] = - 2.0*ef(i)*sin2W ;
366 :
367 : // tmp: verbose output
368 : if (DBSUSY) {
369 : cout << " LllZ [" << i-10 << "][" << i-10 << "] = "
370 : << scientific << setw(10) << LllZ[i-10]
371 : << " RllZ [" << i-10 << "][" << i-10 << "] = "
372 : << scientific << setw(10) << RllZ[i-10] << endl;
373 : }
374 : }
375 :
376 : // Construct ~l~lZ couplings
377 : // Initialize
378 0 : for(int i=1;i<=6;i++){
379 0 : for(int j=1;j<=6;j++){
380 0 : LslslZ[i][j] = 0.0;
381 0 : RslslZ[i][j] = 0.0;
382 :
383 0 : for (int k=1;k<=3;k++) {
384 0 : LslslZ[i][j] += LllZ[1] * Rsl[i][k] * Rsl[j][k];
385 0 : RslslZ[i][j] += RllZ[1] * Rsl[i][k+3] * Rsl[j][k+3];
386 : }
387 :
388 : // ~v[i] ~v[j] Z
389 0 : LsvsvZ[i][j] = 0.0;
390 0 : RsvsvZ[i][j] = 0.0;
391 0 : for (int k=1;k<=3;k++)
392 0 : LsvsvZ[i][j] += LllZ[2] * Rsv[i][k]*Rsv[j][k];
393 : }
394 : }
395 :
396 0 : for(int i=1;i<=6;i++){
397 0 : for(int j=1;j<=6;j++){
398 : if (DBSUSY) {
399 : if (max(abs(LsvsvZ[i][j]),abs(RsvsvZ[i][j])) > 1e-6) {
400 : cout << " LsvsvZ[" << i << "][" << j << "] = "
401 : << scientific << setw(10) << LsvsvZ[i][j]
402 : << " RsvsvZ[" << i << "][" << j << "] = "
403 : << scientific << setw(10) << RsvsvZ[i][j] << endl;
404 : }
405 : if (max(abs(LslslZ[i][j]),abs(RslslZ[i][j]))> 1e-6) {
406 : cout << " LslslZ[" << i << "][" << j << "] = "
407 : << scientific << setw(10) << LslslZ[i][j]
408 : << " RslslZ[" << i << "][" << j << "] = "
409 : << scientific << setw(10) << RslslZ[i][j] << endl;
410 : }
411 : }
412 : }
413 : }
414 :
415 : // Construct udW couplings
416 : // Loop over up [i] and down [j] quark generation
417 0 : for (int i=1;i<=3;i++) {
418 0 : for (int j=1;j<=3;j++) {
419 :
420 : // CKM matrix (use Pythia one if no SLHA)
421 : // (NB: could also try input one if no running one found, but
422 : // would then need to compute from Wolfenstein)
423 0 : complex Vij=VCKMgen(i,j);
424 0 : if (slhaPtr->vckm.exists()) {
425 0 : Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
426 0 : }
427 :
428 : // u[i] d[j] W
429 0 : LudW[i][j] = sqrt(2.0) * cosW * Vij;
430 0 : RudW[i][j] = 0.0;
431 :
432 : // tmp: verbose output
433 : if (DBSUSY) {
434 : cout << " LudW [" << i << "][" << j << "] = "
435 : << scientific << setw(10) << LudW[i][j]
436 : << " RudW [" << i << "][" << j << "] = "
437 : << scientific << setw(10) << RudW[i][j] << endl;
438 : }
439 0 : }
440 : }
441 :
442 :
443 : // Construct ~u~dW couplings
444 : // Loop over ~u[k] and ~d[l] flavours
445 0 : for (int k=1;k<=6;k++) {
446 0 : for (int l=1;l<=6;l++) {
447 :
448 0 : LsusdW[k][l]=0.0;
449 0 : RsusdW[k][l]=0.0;
450 :
451 : // Loop over u[i] and d[j] flavours
452 0 : for (int i=1;i<=3;i++) {
453 0 : for (int j=1;j<=3;j++) {
454 :
455 : // CKM matrix (use Pythia one if no SLHA)
456 : // (NB: could also try input one if no running one found, but
457 : // would then need to compute from Wolfenstein)
458 0 : complex Vij=VCKMgen(i,j);
459 0 : if (slhaPtr->vckm.exists()) {
460 0 : Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
461 0 : }
462 :
463 : // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
464 0 : complex Ruki = complex(Ru(k,i),imRu(k,i));
465 0 : complex Rdlj = complex(Rd(l,j),imRd(l,j));
466 0 : LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
467 0 : RsusdW[k][l] += 0.0;
468 :
469 0 : }
470 : }
471 :
472 : // tmp: verbose output
473 : if (DBSUSY) {
474 : if (max(abs(LsusdW[k][l]),abs(RsusdW[k][l]))> 1e-6) {
475 : cout << " LsusdW[" << k << "][" << l << "] = "
476 : << scientific << setw(10) << LsusdW[k][l]
477 : << " RsusdW[" << k << "][" << l << "] = "
478 : << scientific << setw(10) << RsusdW[k][l] << endl;
479 : }
480 : }
481 :
482 : }
483 : }
484 :
485 :
486 :
487 : // Construct lvW couplings
488 0 : for (int i=1;i<=3;i++){
489 0 : for (int j=1;j<=3;++j){
490 0 : LlvW[i][j] = (i==j) ? sqrt(2.0) * cosW : 0.0 ;
491 0 : RlvW[i][j] = 0.0;
492 :
493 : // tmp: verbose output
494 : if (DBSUSY) {
495 : cout << " LlvW [" << i << "][" << j << "] = "
496 : << scientific << setw(10) << LlvW[i][j]
497 : << " RlvW [" << i << "][" << j << "] = "
498 : << scientific << setw(10) << RlvW[i][j] << endl;
499 : }
500 : }
501 : }
502 :
503 : // Construct ~l~vW couplings
504 0 : for (int k=1;k<=6;k++) {
505 0 : for (int l=1;l<=6;l++) {
506 0 : LslsvW[k][l]=0.0;
507 0 : RslsvW[k][l]=0.0;
508 :
509 0 : if(l<=3) // Only left-handed sneutrinos
510 0 : for(int i=1;i<=3;i++)
511 0 : LslsvW[k][l] += sqrt(2.0) * cosW * Rsl[k][i] * Rsv[l][i];
512 :
513 :
514 : // tmp: verbose output
515 : if (DBSUSY) {
516 : cout << " LslsvW [" << k << "][" << l << "] = "
517 : << scientific << setw(10) << LslsvW[k][l]
518 : << " RslsvW [" << k << "][" << l << "] = "
519 : << scientific << setw(10) << RslsvW[k][l] << endl;
520 : }
521 : }
522 : }
523 :
524 : // Now we come to the ones with really many indices
525 :
526 : // RPV: check for lepton-neutralino mixing (not supported)
527 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
528 : bool hasCrossTerms = false;
529 0 : for (int i=1; i<=3; ++i)
530 0 : for (int j=4; j<=7; ++j)
531 0 : if (abs(slhaPtr->rvnmix(i,j)) > 1e-5) {
532 : hasCrossTerms = true;
533 0 : break;
534 : }
535 0 : if (hasCrossTerms)
536 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
537 : "Neutrino-Neutralino mixing not supported internally in PYTHIA");
538 0 : }
539 :
540 : // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
541 0 : for (int i=1;i<=nNeut;i++) {
542 :
543 : // Ni1, Ni2, Ni3, Ni4, Ni5
544 0 : complex ni1,ni2,ni3,ni4,ni5;
545 :
546 : // In RPV, ignore neutrino-neutralino mixing
547 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
548 0 : ni1=complex( slhaPtr->rvnmix(i+3,4), 0.0 );
549 0 : ni2=complex( slhaPtr->rvnmix(i+3,5), 0.0 );
550 0 : ni3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
551 0 : ni4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
552 0 : ni5=complex( 0.0, 0.0);
553 0 : }
554 0 : else if (!isNMSSM) {
555 0 : ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
556 0 : ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
557 0 : ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
558 0 : ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
559 0 : ni5=complex( 0.0, 0.0);
560 0 : } else {
561 0 : ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
562 0 : ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
563 0 : ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
564 0 : ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
565 0 : ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
566 : }
567 :
568 : // Change to positive mass convention
569 0 : complex iRot( 0., 1.);
570 0 : if (slhaPtr->mass(idNeut(i)) < 0.) {
571 0 : ni1 *= iRot;
572 0 : ni2 *= iRot;
573 0 : ni3 *= iRot;
574 0 : ni4 *= iRot;
575 0 : ni5 *= iRot;
576 0 : }
577 :
578 : // ~chi0 [i] ~chi0 [j] Z : loop over [j]
579 0 : for (int j=1; j<=nNeut; j++) {
580 :
581 : // neutralino [j] higgsino components
582 0 : complex nj3, nj4;
583 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
584 0 : nj3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
585 0 : nj4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
586 0 : }
587 0 : else if (!isNMSSM) {
588 0 : nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
589 0 : nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
590 0 : } else {
591 0 : nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
592 0 : nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
593 : }
594 : // Change to positive mass convention
595 0 : if (slhaPtr->mass(idNeut(j)) < 0.) {
596 0 : nj3 *= iRot;
597 0 : nj4 *= iRot;
598 0 : }
599 :
600 : // ~chi0 [i] ~chi0 [j] Z : couplings
601 0 : OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
602 0 : ORpp[i][j] = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
603 :
604 : // tmp: verbose output
605 : if (DBSUSY) {
606 : cout << " OL'' [" << i << "][" << j << "] = "
607 : << scientific << setw(10) << OLpp[i][j]
608 : << " OR'' [" << i << "][" << j << "] = "
609 : << scientific << setw(10) << ORpp[i][j] << endl;
610 : }
611 :
612 0 : }
613 :
614 : // ~chi0 [i] ~chi+ [j] W : loop over [j]
615 0 : for (int j=1; j<=nChar; j++) {
616 :
617 : // Chargino mixing
618 0 : complex uj1, uj2, vj1, vj2;
619 0 : if (slhaPtr->modsel(4)<1 || !slhaPtr->rvumix.exists()) {
620 0 : uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
621 0 : uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
622 0 : vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
623 0 : vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
624 0 : }
625 : // RPV: ignore lepton-chargino mixing
626 : else {
627 0 : uj1=complex( slhaPtr->rvumix(j+3,4), 0.0 );
628 0 : uj2=complex( slhaPtr->rvumix(j+3,5), 0.0 );
629 0 : vj1=complex( slhaPtr->rvvmix(j+3,4), 0.0 );
630 0 : vj2=complex( slhaPtr->rvvmix(j+3,5), 0.0 );
631 : }
632 :
633 : // ~chi0 [i] ~chi+ [j] W : couplings
634 0 : OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
635 0 : OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
636 :
637 : // tmp: verbose output
638 : if (DBSUSY) {
639 : cout << " OL [" << i << "][" << j << "] = "
640 : << scientific << setw(10) << OL[i][j]
641 : << " OR [" << i << "][" << j << "] = "
642 : << scientific << setw(10) << OR[i][j] << endl;
643 : }
644 0 : }
645 :
646 :
647 : // ~qqX couplings
648 : // Quark Charges
649 : double ed = -1.0/3.0;
650 : double T3d = -0.5;
651 : double eu = 2.0/3.0;
652 : double T3u = 0.5;
653 :
654 : // Loop over quark [k] generation
655 0 : for (int k=1;k<=3;k++) {
656 :
657 : // Set quark masses
658 : // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
659 0 : double mu = particleDataPtr->m0(2*k);
660 0 : double md = particleDataPtr->m0(2*k-1);
661 0 : if (k == 1) { mu=0.0 ; md=0.0; }
662 0 : if (k == 2) { md=0.0 ; mu=0.0; }
663 :
664 : // Compute running mass from Yukawas and vevs if possible.
665 0 : if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
666 0 : double ykk=slhaPtr->yd(k,k);
667 0 : double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
668 0 : if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
669 0 : }
670 0 : if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
671 0 : double ykk=slhaPtr->yu(k,k);
672 0 : double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
673 0 : if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
674 0 : }
675 :
676 : // tmp: verbose output
677 : if (DBSUSY) {
678 : cout << " Gen = " << k << " mu = " << mu << " md = " << md
679 : << " yUU,DD = " << slhaPtr->yu(k,k) << ","
680 : << slhaPtr->yd(k,k) << endl;
681 : }
682 :
683 : // Loop over squark [j] flavour
684 0 : for (int j=1;j<=6;j++) {
685 :
686 : // Squark mixing
687 0 : complex Rdjk = complex(Rd(j,k), imRd(j,k) );
688 0 : complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
689 0 : complex Rujk = complex(Ru(j,k), imRu(j,k) );
690 0 : complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
691 :
692 : // ~d[j] d[k] ~chi0[i]
693 : // Changed according to new notation
694 0 : LsddX[j][k][i] = ((ed-T3d)*sinW*ni1 + T3d*cosW*ni2)*conj(Rdjk)
695 0 : + md*cosW*ni3*conj(Rdjk3)/2.0/mW/cosb;
696 0 : RsddX[j][k][i] = -ed*sinW*conj(ni1)*conj(Rdjk3)
697 0 : + md*cosW*conj(ni3)*conj(Rdjk)/2.0/mW/cosb;
698 :
699 : // ~u[j] u[k] ~chi0[i]
700 0 : LsuuX[j][k][i] = ((eu-T3u)*sinW*ni1 + T3u*cosW*ni2)*conj(Rujk)
701 0 : + mu*cosW*ni4*conj(Rujk3)/2.0/mW/sinb;
702 0 : RsuuX[j][k][i] = -eu*sinW*conj(ni1)*conj(Rujk3)
703 0 : + mu*cosW*conj(ni4)*conj(Rujk)/2.0/mW/sinb;
704 :
705 : if (DBSUSY) {
706 : if (abs(LsddX[j][k][i]) > 1e-6) {
707 : // tmp: verbose output
708 : cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
709 : << scientific << setw(10) << LsddX[j][k][i] << endl;
710 : }
711 : if (abs(RsddX[j][k][i]) > 1e-6) {
712 : // tmp: verbose output
713 : cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
714 : << scientific << setw(10) << RsddX[j][k][i] << endl;
715 : }
716 : if (abs(LsuuX[j][k][i]) > 1e-6) {
717 : // tmp: verbose output
718 : cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
719 : << scientific << setw(10) << LsuuX[j][k][i] << endl;
720 : }
721 : if (abs(RsuuX[j][k][i]) > 1e-6) {
722 : // tmp: verbose output
723 : cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
724 : << scientific << setw(10) << RsuuX[j][k][i] << endl;
725 : }
726 : }
727 0 : }
728 : }
729 :
730 : // Start slepton couplings
731 : // Lepton Charges
732 : double el = -1.0;
733 : double T3l = -0.5;
734 : double ev = 0.0;
735 : double T3v = 0.5;
736 :
737 : // Need to define lepton mass
738 0 : for (int k=1;k<=3;k++) {
739 : // Set lepton masses
740 : double ml(0.0);
741 0 : if(k==3) ml = particleDataPtr->m0(15);
742 :
743 0 : for(int j=1;j<=6;j++){
744 0 : LsllX[j][k][i] = 0;
745 0 : RsllX[j][k][i] = 0;
746 0 : LsvvX[j][k][i] = 0;
747 0 : RsvvX[j][k][i] = 0;
748 :
749 : // ~l[j] l[k] ~chi0[i]
750 : // Changed according to new notation
751 0 : LsllX[j][k][i] = ((el-T3l)*sinW*ni1 + T3l*cosW*ni2)*Rsl[j][k]
752 0 : + ml*cosW*ni3/2.0/mW/cosb*Rsl[j][k+3];
753 0 : RsllX[j][k][i] = -el*sinW*conj(ni1)*Rsl[j][k+3]
754 0 : + ml*cosW*conj(ni3)/2.0/mW/cosb*Rsl[j][k];
755 :
756 0 : if(j<3 && j==k){
757 : // No sneutrino mixing; only left handed
758 : // ~v[j] v[k] ~chi0[i]
759 0 : LsvvX[j][k][i] = ((ev-T3v)*sinW*ni1 + T3v*cosW*ni2);
760 0 : }
761 :
762 : if (DBSUSY) {
763 : if (abs(LsllX[j][k][i]) > 1e-6) {
764 : // tmp: verbose output
765 : cout << " LsllX[" << j << "][" << k << "][" << i << "] = "
766 : << scientific << setw(10) << LsllX[j][k][i] << endl;
767 : }
768 : if (abs(RsllX[j][k][i]) > 1e-6) {
769 : // tmp: verbose output
770 : cout << " RsllX[" << j << "][" << k << "][" << i << "] = "
771 : << scientific << setw(10) << RsllX[j][k][i] << endl;
772 : }
773 : if (abs(LsvvX[j][k][i]) > 1e-6) {
774 : // tmp: verbose output
775 : cout << " LsvvX[" << j << "][" << k << "][" << i << "] = "
776 : << scientific << setw(10) << LsvvX[j][k][i] << endl;
777 : }
778 : }
779 : }
780 : }
781 0 : }
782 :
783 : // RPV: check for lepton-chargino mixing (not supported)
784 0 : if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvumix.exists()) {
785 : bool hasCrossTerms = false;
786 0 : for (int i=1; i<=3; ++i)
787 0 : for (int j=4; j<=5; ++j)
788 0 : if (abs(slhaPtr->rvumix(i,j)) > 1e-5
789 0 : || abs(slhaPtr->rvvmix(i,j)) > 1e-5) {
790 : hasCrossTerms = true;
791 0 : break;
792 : }
793 0 : if (hasCrossTerms)
794 0 : infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
795 : "Lepton-Chargino mixing not supported internally in PYTHIA");
796 0 : }
797 :
798 : // Construct ~chi+ couplings
799 : // sqrt(2)
800 0 : double rt2 = sqrt(2.0);
801 :
802 0 : for (int i=1;i<=nChar;i++) {
803 :
804 : // Ui1, Ui2, Vi1, Vi2
805 0 : complex ui1,ui2,vi1,vi2;
806 0 : ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
807 0 : ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
808 0 : vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
809 0 : vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
810 :
811 : // ~chi+ [i] ~chi- [j] Z : loop over [j]
812 0 : for (int j=1; j<=nChar; j++) {
813 :
814 : // Chargino mixing
815 0 : complex uj1, uj2, vj1, vj2;
816 0 : uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
817 0 : uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
818 0 : vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
819 0 : vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
820 :
821 : // ~chi+ [i] ~chi- [j] Z : couplings
822 0 : OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
823 0 : + ( (i == j) ? sin2W : 0.0);
824 0 : ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
825 0 : + ( (i == j) ? sin2W : 0.0);
826 :
827 : if (DBSUSY) {
828 : // tmp: verbose output
829 : cout << " OL' [" << i << "][" << j << "] = "
830 : << scientific << setw(10) << OLp[i][j]
831 : << " OR' [" << i << "][" << j << "] = "
832 : << scientific << setw(10) << ORp[i][j] << endl;
833 : }
834 0 : }
835 :
836 : // Loop over quark [l] flavour
837 0 : for (int l=1;l<=3;l++) {
838 :
839 : // Set quark [l] masses
840 : // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
841 0 : double mul = particleDataPtr->m0(2*l);
842 0 : double mdl = particleDataPtr->m0(2*l-1);
843 0 : if (l == 1 || l == 2) { mul=0.0 ; mdl=0.0; }
844 :
845 : // Compute running mass from Yukawas and vevs if possible.
846 0 : if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
847 0 : double yll=slhaPtr->yd(l,l);
848 0 : double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
849 0 : if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
850 0 : }
851 0 : if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
852 0 : double yll=slhaPtr->yu(l,l);
853 0 : double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
854 0 : if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
855 0 : }
856 :
857 : // Loop over squark [j] flavour
858 0 : for (int j=1;j<=6;j++) {
859 :
860 : //Initialise to zero
861 0 : LsduX[j][l][i] = 0.0;
862 0 : RsduX[j][l][i] = 0.0;
863 0 : LsudX[j][l][i] = 0.0;
864 0 : RsudX[j][l][i] = 0.0;
865 :
866 : // Loop over off-diagonal quark [k] generation
867 0 : for (int k=1;k<=3;k++) {
868 :
869 : // Set quark [k] masses
870 : // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
871 0 : double muk = particleDataPtr->m0(2*k);
872 0 : double mdk = particleDataPtr->m0(2*k-1);
873 0 : if (k == 1) { muk=0.0 ; mdk=0.0; }
874 0 : if (k == 2) { mdk=0.0 ; muk=0.0; }
875 :
876 : // Compute running mass from Yukawas and vevs if possible.
877 0 : if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
878 0 : double ykk=slhaPtr->yd(k,k);
879 0 : double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
880 0 : if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
881 0 : }
882 0 : if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
883 0 : double ykk=slhaPtr->yu(k,k);
884 0 : double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
885 0 : if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
886 0 : }
887 :
888 : // CKM matrix (use Pythia one if no SLHA)
889 : // (NB: could also try input one if no running one found, but
890 : // would then need to compute from Wolfenstein)
891 0 : complex Vlk=VCKMgen(l,k);
892 0 : complex Vkl=VCKMgen(k,l);
893 0 : if (slhaPtr->vckm.exists()) {
894 0 : Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
895 0 : Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
896 0 : }
897 :
898 : // Squark mixing
899 0 : complex Rdjk = complex(Rd(j,k), imRd(j,k) );
900 0 : complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
901 0 : complex Rujk = complex(Ru(j,k), imRu(j,k) );
902 0 : complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
903 :
904 :
905 : // ~d[j] u[l] ~chi+[i]
906 0 : LsduX[j][l][i] += (ui1*conj(Rdjk)
907 0 : - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
908 0 : RsduX[j][l][i] -= mul*conj(vi2)*Vlk*conj(Rdjk)/rt2/mW/sinb;
909 :
910 : // ~u[j] d[l] ~chi+[i]
911 0 : LsudX[j][l][i] += (conj(vi1)*Rujk
912 0 : - muk*conj(vi2)*Rujk3/rt2/mW/sinb)*Vkl;
913 0 : RsudX[j][l][i] -= mdl*ui2*Vkl*Rujk/rt2/mW/cosb;
914 :
915 0 : }
916 :
917 : if (DBSUSY) {
918 : if (max(abs(LsduX[j][l][i]),abs(RsduX[j][l][i])) > 1e-6) {
919 : // tmp: verbose output
920 : cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
921 : << scientific << setw(10) << LsduX[j][l][i];
922 : cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
923 : << scientific << setw(10) << RsduX[j][l][i] << endl;
924 : }
925 : if (max(abs(LsudX[j][l][i]),abs(RsudX[j][l][i])) > 1e-6) {
926 : // tmp: verbose output
927 : cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
928 : << scientific << setw(10) << LsudX[j][l][i];
929 : cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
930 : << scientific << setw(10) << RsudX[j][l][i] << endl;
931 : }
932 : }
933 : }
934 0 : }
935 :
936 : // Loop over slepton [j] flavour
937 0 : for (int j=1;j<=6;j++) {
938 0 : for (int k=1;k<=3;k++) {
939 :
940 0 : LslvX[j][k][i] = 0.0;
941 0 : RslvX[j][k][i] = 0.0;
942 0 : LsvlX[j][k][i] = 0.0;
943 0 : RsvlX[j][k][i] = 0.0;
944 :
945 : // Set lepton [k] masses
946 0 : double ml = 0.0;
947 0 : if (k == 3) ml = particleDataPtr->m0(15);
948 :
949 0 : if(j==k || j==k+3){ // No lepton mixing
950 :
951 : // ~l[j] v[l] ~chi+[i]
952 0 : LslvX[j][k][i] += ui1- ml*ui2/rt2/mW/cosb;
953 0 : RslvX[j][k][i] -= ml*conj(vi2)/rt2/mW/sinb;
954 :
955 : // ~v[j] l[l] ~chi+[i]
956 0 : if(j<=3){ // No right handed sneutrinos
957 0 : LsvlX[j][k][i] += conj(vi1) - ml*conj(vi2)/rt2/mW/sinb;
958 0 : }
959 : }
960 :
961 : if (DBSUSY) {
962 : if (max(abs(LslvX[j][k][i]),abs(RslvX[j][k][i])) > 1e-6) {
963 : // tmp: verbose output
964 : cout << " LslvX[" << j << "][" << k << "][" << i << "] = "
965 : << scientific << setw(10) << LslvX[j][k][i];
966 : cout << " RslvX[" << j << "][" << k << "][" << i << "] = "
967 : << scientific << setw(10) << RslvX[j][k][i] << endl;
968 : }
969 : if (max(abs(LsvlX[j][k][i]),abs(RsvlX[j][k][i])) > 1e-6) {
970 : // tmp: verbose output
971 : cout << " LsvlX[" << j << "][" << k << "][" << i << "] = "
972 : << scientific << setw(10) << LsvlX[j][k][i];
973 : cout << " RsvlX[" << j << "][" << k << "][" << i << "] = "
974 : << scientific << setw(10) << RsvlX[j][k][i] << endl;
975 : }
976 : }
977 0 : }
978 : }
979 0 : }
980 :
981 : // Shorthand for RPV couplings
982 : // The input LNV lambda couplings
983 0 : LHtensor3Block<3> rvlle(slhaPtr->rvlamlle);
984 : // The input LNV lambda' couplings
985 0 : LHtensor3Block<3> rvlqd(slhaPtr->rvlamlqd);
986 : // The input BNV lambda'' couplings
987 0 : LHtensor3Block<3> rvudd(slhaPtr->rvlamudd);
988 :
989 0 : isLLE = false;
990 0 : isLQD = false;
991 0 : isUDD = false;
992 :
993 : // Symmetry properties
994 0 : if(rvlle.exists()){
995 0 : isLLE = true;
996 0 : for(int i=1;i<=3;i++){
997 0 : for(int j=i;j<=3;j++){
998 : //lambda(i,j,k)=-lambda(j,i,k)
999 0 : for(int k=1;k<=3;k++){
1000 0 : if(i==j){
1001 0 : rvLLE[i][j][k] = 0.0;
1002 0 : }else{
1003 0 : rvLLE[i][j][k] = rvlle(i,j,k);
1004 0 : rvLLE[j][i][k] = -rvlle(i,j,k);
1005 : }
1006 : }
1007 : }
1008 : }
1009 0 : }
1010 :
1011 0 : if(rvlqd.exists()){
1012 0 : isLQD=true;
1013 0 : for(int i=1;i<=3;i++){
1014 0 : for(int j=1;j<=3;j++){
1015 0 : for(int k=1;k<=3;k++){
1016 0 : rvLQD[i][j][k] = rvlqd(i,j,k);
1017 : }
1018 : }
1019 : }
1020 0 : }
1021 :
1022 : //lambda''(k,j,i)=-lambda''(k,i,j)
1023 :
1024 0 : if(rvudd.exists()){
1025 0 : isUDD = true;
1026 0 : for(int i=1;i<=3;i++){
1027 0 : for(int j=i;j<=3;j++){
1028 0 : for(int k=1;k<=3;k++){
1029 0 : if(i==j){
1030 0 : rvUDD[k][i][j] = 0.0;
1031 0 : }else{
1032 0 : rvUDD[k][i][j] = rvudd(k,i,j);
1033 0 : rvUDD[k][j][i] = -rvudd(k,i,j);
1034 : }
1035 : }
1036 : }
1037 : }
1038 0 : }
1039 :
1040 : if(DBSUSY){
1041 : for(int i=1;i<=3;i++){
1042 : for(int j=1;j<=3;j++){
1043 : for(int k=1;k<=3;k++){
1044 : if(rvlle.exists())
1045 : cout<<"LambdaLLE["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLLE[i][j][k]<<" ";
1046 : if(rvlqd.exists())
1047 : cout<<"LambdaLQD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLQD[i][j][k]<<" ";
1048 : if(rvudd.exists())
1049 : cout<<"LambdaUDD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvUDD[i][j][k]
1050 : <<"\n";
1051 : }
1052 : }
1053 : }
1054 : }
1055 :
1056 : // Store the squark mixing matrix
1057 0 : for(int i=1;i<=6;i++){
1058 0 : for(int j=1;j<=3;j++){
1059 0 : Rusq[i][j] = complex(Ru(i,j), imRu(i,j) );
1060 0 : Rusq[i][j+3] = complex(Ru(i,j+3),imRu(i,j+3));
1061 0 : Rdsq[i][j] = complex(Rd(i,j), imRd(i,j) );
1062 0 : Rdsq[i][j+3] = complex(Rd(i,j+3),imRd(i,j+3));
1063 : }
1064 : }
1065 :
1066 : if(DBSUSY){
1067 : cout<<"Ru"<<endl;
1068 : for(int i=1;i<=6;i++){
1069 : for(int j=1;j<=6;j++){
1070 : cout << real(Rusq[i][j])<<" ";
1071 : }
1072 : cout<<endl;
1073 : }
1074 : cout<<"Rd"<<endl;
1075 : for(int i=1;i<=6;i++){
1076 : for(int j=1;j<=6;j++){
1077 : cout << real(Rdsq[i][j])<<" ";
1078 : }
1079 : cout<<endl;
1080 : }
1081 : }
1082 :
1083 : // Let everyone know we are ready
1084 0 : isInit = true;
1085 0 : }
1086 :
1087 : //--------------------------------------------------------------------------
1088 :
1089 : // Return neutralino flavour codes.
1090 :
1091 : int CoupSUSY::idNeut(int idChi) {
1092 :
1093 : int id = 0;
1094 0 : if (idChi == 1) id = 1000022;
1095 0 : else if (idChi == 2) id = 1000023;
1096 0 : else if (idChi == 3) id = 1000025;
1097 0 : else if (idChi == 4) id = 1000035;
1098 0 : else if (idChi == 5) id = 1000045;
1099 0 : return id;
1100 : }
1101 :
1102 : //--------------------------------------------------------------------------
1103 :
1104 : // Return chargino flavour codes.
1105 :
1106 : int CoupSUSY::idChar(int idChi) {
1107 :
1108 : int id = 0;
1109 0 : if (idChi == 1) id = 1000024;
1110 0 : else if (idChi == -1) id = -1000024;
1111 0 : else if (idChi == 2) id = 1000037;
1112 0 : else if (idChi == -2) id = -1000037;
1113 0 : return id;
1114 : }
1115 :
1116 : //--------------------------------------------------------------------------
1117 :
1118 : // Return sup flavour codes.
1119 :
1120 : int CoupSUSY::idSup(int iSup) {
1121 :
1122 : int id = 0;
1123 0 : int sgn = ( iSup > 0 ) ? 1 : -1;
1124 0 : iSup = abs(iSup);
1125 0 : if (iSup == 1) id = 1000002;
1126 0 : else if (iSup == 2) id = 1000004;
1127 0 : else if (iSup == 3) id = 1000006;
1128 0 : else if (iSup == 4) id = 2000002;
1129 0 : else if (iSup == 5) id = 2000004;
1130 0 : else if (iSup == 6) id = 2000006;
1131 0 : return sgn*id;
1132 : }
1133 :
1134 : //--------------------------------------------------------------------------
1135 :
1136 : // Return sdown flavour codes
1137 :
1138 : int CoupSUSY::idSdown(int iSdown) {
1139 :
1140 : int id = 0;
1141 0 : int sgn = ( iSdown > 0 ) ? 1 : -1;
1142 0 : iSdown = abs(iSdown);
1143 0 : if (iSdown == 1) id = 1000001;
1144 0 : else if (iSdown == 2) id = 1000003;
1145 0 : else if (iSdown == 3) id = 1000005;
1146 0 : else if (iSdown == 4) id = 2000001;
1147 0 : else if (iSdown == 5) id = 2000003;
1148 0 : else if (iSdown == 6) id = 2000005;
1149 0 : return sgn*id;
1150 : }
1151 :
1152 : //--------------------------------------------------------------------------
1153 :
1154 : // Function to return slepton flavour codes
1155 :
1156 : int CoupSUSY::idSlep(int iSlep) {
1157 :
1158 : int id = 0;
1159 0 : int sgn = ( iSlep > 0 ) ? 1 : -1;
1160 0 : iSlep = abs(iSlep);
1161 0 : if (iSlep == 1) id = 1000011;
1162 0 : else if (iSlep == 2) id = 1000013;
1163 0 : else if (iSlep == 3) id = 1000015;
1164 0 : else if (iSlep == 4) id = 2000011;
1165 0 : else if (iSlep == 5) id = 2000013;
1166 0 : else if (iSlep == 6) id = 2000015;
1167 0 : return sgn*id;
1168 : }
1169 :
1170 : //--------------------------------------------------------------------------
1171 :
1172 : // Return neutralino code; zero if not a (recognized) neutralino
1173 :
1174 : int CoupSUSY::typeNeut(int idPDG) {
1175 : int type = 0;
1176 0 : int idAbs = abs(idPDG);
1177 0 : if(idAbs == 1000022) type = 1;
1178 0 : else if(idAbs == 1000023) type = 2;
1179 0 : else if(idAbs == 1000025) type = 3;
1180 0 : else if(idAbs == 1000035) type = 4;
1181 0 : else if(isNMSSM && idAbs == 1000045) type = 5;
1182 0 : return type;
1183 :
1184 : }
1185 :
1186 : //--------------------------------------------------------------------------
1187 :
1188 : // Check whether particle is a Chargino
1189 :
1190 : int CoupSUSY::typeChar(int idPDG) {
1191 : int type = 0;
1192 0 : if(abs(idPDG) == 1000024) type = 1;
1193 0 : else if (abs(idPDG) == 1000037)type = 2;
1194 0 : return type;
1195 : }
1196 :
1197 : //==========================================================================
1198 :
1199 : } // end namespace Pythia8
|
