- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / SusyCouplings.cxx
1 // SusyCouplings.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 Peter Skands, 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 // supersymmetric couplings class. 
8
9 #include "SusyCouplings.h"
10
11 namespace Pythia8 {
12
13 //==========================================================================
14
15 // The CoupSUSY class.
16
17 //--------------------------------------------------------------------------
18
19 // Constants: could be changed here if desired, but normally should not.
20 // These are of technical nature, as described for each.
21
22 // Allow verbose printout for debug purposes.
23 const bool CoupSUSY::DEBUG = false;
24
25 //--------------------------------------------------------------------------
26
27 // Initialize SM+SUSY couplings (only performed once).
28
29 void CoupSUSY::init (SusyLesHouches* slhaPtrIn, Settings* settingsPtrIn, 
30   ParticleData* particleDataPtrIn, CoupSM* coupSMPtrIn) {
31
32   // Save pointers..
33   slhaPtr         = slhaPtrIn;
34   settingsPtr     = settingsPtrIn;
35   particleDataPtr = particleDataPtrIn;
36   coupSMPtr       = coupSMPtrIn;
37   
38   // Is NMSSM switched on?
39   isNMSSM = (slhaPtr->modsel(3) != 1 ? false : true);
40   settingsPtr->flag("SLHA:NMSSM",isNMSSM);
41   int nNeut = (isNMSSM ? 5 : 4);
42   int nChar = 2;
43
44   // Initialize pole masses 
45   mZpole    = particleDataPtr->m0(23);
46   wZpole    = particleDataPtr->mWidth(23);
47   mWpole    = particleDataPtr->m0(24);
48   wWpole    = particleDataPtr->mWidth(24);
49   
50   // Running masses and weak mixing angle 
51   // (default to pole values if no running available)
52   mW        = mWpole;
53   mZ        = mZpole;
54   sin2W     = 1.0 - pow(mW/mZ,2);  
55
56   if (slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2) 
57       && slhaPtr->hmix.exists(3)) {
58     double gp=slhaPtr->gauge(1);
59     double g =slhaPtr->gauge(2);
60     double v =slhaPtr->hmix(3);
61     mW      = g * v / 2.0;
62     mZ      = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
63     double tan2W   = pow2(gp)/pow2(g);
64     if (DEBUG) cout << " tan2W = " << tan2W << endl;
65     sin2W   = pow2(gp)/(pow2(g)+pow2(gp));  
66   }
67   sinW = sqrt(sin2W);
68   cosW = sqrt(1.0-sin2W);
69
70   // Tan(beta)
71   // By default, use the running one in HMIX (if not found, use MINPAR)
72   tanb = slhaPtr->hmix.exists(2) ? slhaPtr->hmix(2) : slhaPtr->minpar(3);
73   cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
74   sinb = sqrt(max(0.0,1.0-cosb*cosb));
75   
76   // tmp : verbose output
77   if (DEBUG) {
78     cout << " sin2W(Q) = " << sin2W << "  mW(Q) = " << mW 
79          << "  mZ(Q) = " << mZ << endl;
80     cout << " vev(Q) = " << slhaPtr->hmix(3) << " tanb(Q) = " << tanb
81          << endl;
82     for (int i=1;i<=3;i++) {
83       for (int j=1;j<=3;j++) {
84         cout << " VCKM  [" << i << "][" << j << "] = " 
85              << scientific << setw(10) << coupSMPtr->VCKMgen(i,j) << endl;
86       }
87     }
88   }  
89   
90   // Shorthand for squark mixing matrices 
91   SusyLesHouches::matrixblock<6> Ru(slhaPtr->usqmix);
92   SusyLesHouches::matrixblock<6> Rd(slhaPtr->dsqmix);
93   SusyLesHouches::matrixblock<6> imRu(slhaPtr->imusqmix);
94   SusyLesHouches::matrixblock<6> imRd(slhaPtr->imusqmix);
95   
96   // Construct ~g couplings
97   for (int i=1 ; i<=6 ; i++) {
98     for (int j=1 ; j<=3 ; j++) {
99       LsddG[i][j] = complex( Rd(i,j)  ,  imRd(i,j));
100       RsddG[i][j] = complex(-Rd(i,j+3), -imRd(i,j+3));
101       LsuuG[i][j] = complex( Ru(i,j)  ,  imRu(i,j));
102       RsuuG[i][j] = complex(-Ru(i,j+3), -imRu(i,j+3));
103     }
104   }
105   
106   // Construct qqZ couplings
107   for (int i=1 ; i<=6 ; i++) {
108     
109     // q[i] q[i] Z (def with extra factor 2 compared to [Okun])
110     LqqZ[i] = coupSMPtr->af(i) - 2.0*coupSMPtr->ef(i)*sin2W ;
111     RqqZ[i] =                  - 2.0*coupSMPtr->ef(i)*sin2W ;
112
113     // tmp: verbose output
114     if (DEBUG) {
115       cout << " LqqZ  [" << i << "][" << i << "] = " 
116            << scientific << setw(10) << LqqZ[i] 
117            << " RqqZ  [" << i << "][" << i  << "] = " 
118            << scientific << setw(10) << RqqZ[i] << endl;
119     }
120   }
121
122   // Construct ~q~qZ couplings
123   for (int i=1 ; i<=6 ; i++) {
124
125     // Squarks can have off-diagonal couplings as well
126     for (int j=1 ; j<=6 ; j++) {
127       
128       // ~d[i] ~d[j] Z
129       LsdsdZ[i][j] = 0.0;
130       RsdsdZ[i][j] = 0.0;
131       for (int k=1;k<=3;k++) {
132         complex Rdik  = complex(Rd(i,k),  imRd(i,k)  );
133         complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
134         complex Rdik3 = complex(Rd(i,k+3),imRd(i,k+3));
135         complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
136         LsdsdZ[i][j] += LqqZ[1] * (Rdik*conj(Rdjk)); 
137         RsdsdZ[i][j] += RqqZ[1] * (Rdik3*conj(Rdjk3)); 
138       }
139       
140       // ~u[i] ~u[j] Z
141       LsusuZ[i][j] = 0.0;
142       RsusuZ[i][j] = 0.0; 
143       for (int k=1;k<=3;k++) {
144         complex Ruik  = complex(Ru(i,k)  ,imRu(i,k)  );
145         complex Rujk  = complex(Ru(j,k)  ,imRu(j,k)  );
146         complex Ruik3 = complex(Ru(i,k+3),imRu(i,k+3));
147         complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
148         LsusuZ[i][j] += LqqZ[2] * (Ruik*conj(Rujk)); 
149         RsusuZ[i][j] += RqqZ[2] * (Ruik3*conj(Rujk3)); 
150       }
151       
152     // tmp: verbose output
153       if (DEBUG) {
154         if (max(abs(LsdsdZ[i][j]),abs(RsdsdZ[i][j])) > 1e-6) {
155           cout << " LsdsdZ[" << i << "][" << j << "] = " 
156                << scientific << setw(10) << LsdsdZ[i][j]
157                << " RsdsdZ[" << i << "][" << j << "] = " 
158                << scientific << setw(10) << RsdsdZ[i][j] << endl;
159         }
160         if (max(abs(LsusuZ[i][j]),abs(RsusuZ[i][j]))> 1e-6) {
161           cout << " LsusuZ[" << i << "][" << j << "] = " 
162                << scientific << setw(10) << LsusuZ[i][j]
163                << " RsusuZ[" << i << "][" << j << "] = " 
164                << scientific << setw(10) << RsusuZ[i][j] << endl;
165         }
166       }
167         
168     }
169     
170   }
171   
172   // Construct udW couplings
173   
174   // Loop over up [i] and down [j] quark generation
175   for (int i=1;i<=3;i++) {
176     for (int j=1;j<=3;j++) {
177       
178       // CKM matrix (use Pythia one if no SLHA)
179       // (NB: could also try input one if no running one found, but
180       // would then need to compute from Wolfenstein)
181       complex Vij=coupSMPtr->VCKMgen(i,j);
182       if (slhaPtr->vckm.exists()) {
183         Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
184       }
185       
186       // u[i] d[j] W  
187       LudW[i][j] = sqrt(2.0) * cosW * Vij;
188       RudW[i][j] = 0.0;
189             
190       // tmp: verbose output
191       if (DEBUG) {
192         cout << " LudW  [" << i << "][" << j << "] = " 
193              << scientific << setw(10) << LudW[i][j]
194              << " RudW  [" << i << "][" << j << "] = " 
195              << scientific << setw(10) << RudW[i][j] << endl;
196       }
197     }
198   }
199
200   // Construct ~u~dW couplings
201
202   // Loop over ~u[k] and ~d[l] flavours
203   for (int k=1;k<=6;k++) {
204     for (int l=1;l<=6;l++) {
205           
206       LsusdW[k][l]=0.0; 
207       RsusdW[k][l]=0.0;
208
209       // Loop over u[i] and d[j] flavours
210       for (int i=1;i<=3;i++) { 
211         for (int j=1;j<=3;j++) {
212           
213           // CKM matrix (use Pythia one if no SLHA)
214           // (NB: could also try input one if no running one found, but
215           // would then need to compute from Wolfenstein)
216           complex Vij=coupSMPtr->VCKMgen(i,j);
217           if (slhaPtr->vckm.exists()) {
218             Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
219           }
220       
221           // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
222           complex Ruki = complex(Ru(k,i),imRu(k,i));
223           complex Rdlj = complex(Rd(l,j),imRd(l,j));
224           LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
225           RsusdW[k][l] += 0.0;
226           
227         }
228       }
229
230       // tmp: verbose output
231       if (DEBUG) {
232         if (max(abs(LsusdW[k][l]),abs(RsusdW[k][l]))> 1e-6) {
233           cout << " LsusdW[" << k << "][" << l << "] = " 
234                << scientific << setw(10) << LsusdW[k][l]
235                << " RsusdW[" << k << "][" << l << "] = " 
236                << scientific << setw(10) << RsusdW[k][l] << endl;
237         }
238       }
239
240     }
241   }
242   
243   // Now we come to the ones with really many indices
244   
245   // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
246   for (int i=1;i<=nNeut;i++) {
247     
248     // Ni1, Ni2, Ni3, Ni4, Ni5
249     complex ni1,ni2,ni3,ni4,ni5;
250     if (not isNMSSM) {  
251       ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
252       ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
253       ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
254       ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
255       ni5=complex( 0.0, 0.0);
256     } else {
257       ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
258       ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
259       ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
260       ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
261       ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
262     }
263     
264     // Change to positive mass convention
265     complex iRot( 0., 1.);
266     if (slhaPtr->mass(idNeut(i)) < 0.) {
267       ni1 *= iRot;
268       ni2 *= iRot;
269       ni3 *= iRot;
270       ni4 *= iRot;
271       ni5 *= iRot;
272     }
273     
274     // ~chi0 [i] ~chi0 [j] Z : loop over [j]
275     for (int j=1; j<=nNeut; j++) {
276       
277       // neutralino [j] higgsino components
278       complex nj3, nj4;
279       if (not isNMSSM) {
280         nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
281         nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
282       } else {
283         nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
284         nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
285       }
286       // Change to positive mass convention
287       if (slhaPtr->mass(idNeut(j)) < 0.) {
288         nj3 *= iRot;
289         nj4 *= iRot;
290       }
291       
292       // ~chi0 [i] ~chi0 [j] Z : couplings
293       OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
294       ORpp[i][j] =  0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
295       
296     // tmp: verbose output
297       if (DEBUG) {
298         cout << " OL''  [" << i << "][" << j << "] = " 
299              << scientific << setw(10) << OLpp[i][j]
300              << " OR''  [" << i << "][" << j << "] = " 
301              << scientific << setw(10) << ORpp[i][j] << endl;
302       }
303         
304     }
305     
306     // ~chi0 [i] ~chi+ [j] W : loop over [j]
307     for (int j=1; j<=nChar; j++) {
308       
309       // Chargino mixing
310       complex uj1, uj2, vj1, vj2;      
311       uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
312       uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
313       vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
314       vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
315       
316       // ~chi0 [i] ~chi+ [j] W : couplings
317       OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
318       OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
319       
320       // tmp: verbose output
321       if (DEBUG) {
322         cout << " OL    [" << i << "][" << j << "] = " 
323              << scientific << setw(10) << OL[i][j]
324              << " OR    [" << i << "][" << j << "] = " 
325              << scientific << setw(10) << OR[i][j] << endl;
326       }
327     }
328     
329     // Charges
330     double ed  = -1.0/3.0;
331     double T3d = -0.5;
332     double eu  =  2.0/3.0;
333     double T3u =  0.5;
334     
335     // Loop over quark [k] generation
336     for (int k=1;k<=3;k++) {
337       
338       // Set quark masses
339       // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
340       double mu = particleDataPtr->m0(2*k);
341       double md = particleDataPtr->m0(2*k-1);
342       if (k == 1) { mu=0.0 ; md=0.0; }
343       if (k == 2) { md=0.0 ; mu=0.0; } 
344       
345       // Compute running mass from Yukawas and vevs if possible.
346       if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
347         double ykk=slhaPtr->yd(k,k);
348         double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
349         if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
350       }
351       if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
352         double ykk=slhaPtr->yu(k,k);
353         double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
354         if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
355       }
356       
357       // tmp: verbose output
358       if (DEBUG) {
359         cout  <<  " Gen = " << k << " mu = " << mu << " md = " << md 
360               << " yUU,DD = " << slhaPtr->yu(k,k) << "," 
361               << slhaPtr->yd(k,k) << endl;
362       }
363       
364       // Loop over squark [j] flavour
365       for (int j=1;j<=6;j++) {
366         
367         // Squark mixing
368         complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
369         complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
370         complex Rujk  = complex(Ru(j,k),  imRu(j,k)  );
371         complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
372         
373         // ~d[j] d[k] ~chi0[i]
374         double rt2 = sqrt(2.0);
375         LsddX[j][k][i] = ((ed-T3d)*sinW/cosW*ni1 + T3d*ni2)*conj(Rdjk)/rt2
376           + md*ni3*conj(Rdjk3)/2.0/rt2/mW/cosb; 
377         RsddX[j][k][i] = -ed*sinW/cosW*conj(ni1)*conj(Rdjk3)/rt2 
378           + md*conj(ni3)*conj(Rdjk)/2.0/rt2/mW/cosb;
379
380         // ~u[j] u[k] ~chi0[i]
381         LsuuX[j][k][i] = ((eu-T3u)*sinW/cosW*ni1 + T3u*ni2)*conj(Rujk)/rt2
382           + mu*ni4*conj(Rujk3)/2.0/rt2/mW/sinb;
383         RsuuX[j][k][i] = -eu*sinW/cosW*conj(ni1)*conj(Rujk3)/rt2
384           + mu*conj(ni4)*conj(Rujk)/2.0/rt2/mW/sinb;
385
386         if (DEBUG) {
387           if (abs(LsddX[j][k][i]) > 1e-6) {
388             // tmp: verbose output
389             cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
390                  << scientific << setw(10) << LsddX[j][k][i] << endl;
391           }
392           if (abs(RsddX[j][k][i]) > 1e-6) {
393             // tmp: verbose output
394             cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
395                  << scientific << setw(10) << RsddX[j][k][i] << endl;
396           }
397           if (abs(LsuuX[j][k][i]) > 1e-6) {
398             // tmp: verbose output
399             cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
400                  << scientific << setw(10) << LsuuX[j][k][i] << endl;
401           }
402           if (abs(RsuuX[j][k][i]) > 1e-6) {
403             // tmp: verbose output
404             cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
405                  << scientific << setw(10) << RsuuX[j][k][i] << endl;
406           }
407         }
408       }
409       
410     }
411     
412   }
413   
414   // Construct ~chi+ couplings 
415   for (int i=1;i<=nChar;i++) {
416     
417     // Ui1, Ui2, Vi1, Vi2
418     complex ui1,ui2,vi1,vi2;    
419     ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
420     ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
421     vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
422     vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
423     
424     // ~chi+ [i] ~chi- [j] Z : loop over [j]
425     for (int j=1; j<=nChar; j++) {
426       
427       // Chargino mixing
428       complex uj1, uj2, vj1, vj2;      
429       uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
430       uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
431       vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
432       vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
433       
434       // ~chi+ [i] ~chi- [j] Z : couplings
435       OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2) 
436         + ( (i == j) ? sin2W : 0.0);
437       ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2 
438         + ( (i == j) ? sin2W : 0.0);
439       
440       if (DEBUG) {
441         // tmp: verbose output
442         cout << " OL'   [" << i << "][" << j << "] = " 
443              << scientific << setw(10) << OLp[i][j]
444              << " OR'   [" << i << "][" << j << "] = " 
445              << scientific << setw(10) << ORp[i][j] << endl;
446       }
447     }
448     
449     // Loop over quark [l] flavour
450     for (int l=1;l<=3;l++) {
451       
452       // Set quark [l] masses 
453       // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
454       double mul = particleDataPtr->m0(2*l);
455       double mdl = particleDataPtr->m0(2*l-1);
456       if (l == 1) { mul=0.0 ; mdl=0.0; }
457       if (l == 2) { mdl=0.0 ; mul=0.0; } 
458       
459       // Compute running mass from Yukawas and vevs if possible.
460       if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
461         double yll=slhaPtr->yd(l,l);
462         double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
463         if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
464       }
465       if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
466         double yll=slhaPtr->yu(l,l);
467         double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
468         if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
469       }
470       
471       // Loop over squark [j] flavour
472       for (int j=1;j<=6;j++) {
473
474         // Loop over off-diagonal quark [k] generation
475         for (int k=1;k<=3;k++) {
476
477           // Set quark [k] masses 
478           // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
479           double muk = particleDataPtr->m0(2*k);
480           double mdk = particleDataPtr->m0(2*k-1);
481           if (k == 1) { muk=0.0 ; mdk=0.0; }
482           if (k == 2) { mdk=0.0 ; muk=0.0; } 
483       
484           // Compute running mass from Yukawas and vevs if possible.
485           if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
486             double ykk=slhaPtr->yd(k,k);
487             double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
488             if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
489           }
490           if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
491             double ykk=slhaPtr->yu(k,k);
492             double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
493             if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
494           }       
495
496           // CKM matrix (use Pythia one if no SLHA)
497           // (NB: could also try input one if no running one found, but
498           // would then need to compute from Wolfenstein)
499           complex Vlk=coupSMPtr->VCKMgen(l,k);
500           complex Vkl=coupSMPtr->VCKMgen(k,l);
501           if (slhaPtr->vckm.exists()) {
502             Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
503             Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
504           }
505       
506           // Squark mixing
507           complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
508           complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
509           complex Rujk  = complex(Ru(j,k),  imRu(j,k)  );
510           complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
511
512           // sqrt(2)
513           double rt2 = sqrt(2.0);
514           
515           // ~d[j] u[l] ~chi+[i]
516           LsduX[j][l][i] += 0.5*(ui1*conj(Rdjk) 
517                                  - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
518           RsduX[j][l][i] -= 0.5*mul*conj(vi2)*Vlk*Rdjk/rt2/mW/sinb; 
519
520           // ~u[j] d[l] ~chi+[i]
521           LsudX[j][l][i] += 0.5*(vi1*conj(Rujk)
522                                  - muk*vi2*conj(Rujk3)/rt2/mW/sinb)*conj(Vkl);
523           RsudX[j][l][i] -= 0.5*mdl*conj(ui2)*conj(Vkl)*conj(Rujk)/rt2/mW/cosb;
524
525         }
526
527         if (DEBUG) {
528           if (max(abs(LsduX[j][l][i]),abs(RsduX[j][l][i])) > 1e-6) {
529             // tmp: verbose output
530             cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
531                  << scientific << setw(10) << LsduX[j][l][i];
532             cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
533                  << scientific << setw(10) << RsduX[j][l][i] << endl;
534           }
535           if (max(abs(LsudX[j][l][i]),abs(RsudX[j][l][i])) > 1e-6) {
536             // tmp: verbose output
537             cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
538                  << scientific << setw(10) << LsudX[j][l][i];
539             cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
540                  << scientific << setw(10) << RsudX[j][l][i] << endl;
541           }
542         }
543         
544       }
545       
546       
547     }
548     
549   }
550
551   // SLHA2 compatibility
552   // Check whether scalar particle masses are ordered
553   bool orderedQ = true;
554   bool orderedL = true;
555   for (int j=1;j<=5;j++) {
556     if (particleDataPtr->m0(idSlep(j+1)) < particleDataPtr->m0(idSlep(j))) 
557       orderedL  = false;
558     if (particleDataPtr->m0(idSup(j+1)) < particleDataPtr->m0(idSup(j)))
559       orderedQ  = false;
560     if (particleDataPtr->m0(idSdown(j+1)) <particleDataPtr->m0(idSdown(j)))
561       orderedQ  = false;
562   }
563
564   // If ordered, change sparticle labels to mass-ordered enumeration
565   for (int i=1;i<=6;i++) {
566     ostringstream indx;
567     indx << i;
568     string uName = "~u_"+indx.str();
569     string dName = "~d_"+indx.str();
570     string lName = "~e_"+indx.str();
571     ParticleDataEntry* pdePtr;
572     if (orderedQ) {
573       pdePtr = particleDataPtr->particleDataEntryPtr(idSup(i));
574       pdePtr->setNames(uName,uName+"bar");
575       pdePtr = particleDataPtr->particleDataEntryPtr(idSdown(i));
576       pdePtr->setNames(dName,dName+"bar");
577     }
578     if (orderedL) {
579       pdePtr = particleDataPtr->particleDataEntryPtr(idSlep(i));
580       pdePtr->setNames(lName,lName+"bar");
581     }
582   }
583                         
584   // Let everyone know we are ready
585   isInit = true;
586   
587 }
588
589 //--------------------------------------------------------------------------
590
591 // Return neutralino flavour codes.
592
593 int CoupSUSY::idNeut(int idChi) {
594   int id = 0;
595   if      (idChi == 1) id = 1000022; 
596   else if (idChi == 2) id = 1000023; 
597   else if (idChi == 3) id = 1000025; 
598   else if (idChi == 4) id = 1000035; 
599   else if (idChi == 5) id = 1000045; 
600   return id;
601 }
602
603 //--------------------------------------------------------------------------
604
605 // Return chargino flavour codes.
606
607 int CoupSUSY::idChar(int idChi) {
608   int id = 0;
609   if      (idChi ==  1) id =  1000024; 
610   else if (idChi == -1) id = -1000024;     
611   else if (idChi ==  2) id =  1000037; 
612   else if (idChi == -2) id = -1000037; 
613   return id;
614 }
615
616 //--------------------------------------------------------------------------
617
618 // Return sup flavour codes.
619
620 int CoupSUSY::idSup(int iSup) {
621   int id = 0;
622   int sgn = ( iSup > 0 ) ? 1 : -1;
623   iSup = abs(iSup);
624   if      (iSup ==  1) id =  1000002; 
625   else if (iSup ==  2) id =  1000004;     
626   else if (iSup ==  3) id =  1000006; 
627   else if (iSup ==  4) id =  2000002; 
628   else if (iSup ==  5) id =  2000004;     
629   else if (iSup ==  6) id =  2000006; 
630   return sgn*id;
631 }
632
633 //--------------------------------------------------------------------------
634
635 // Return sdown flavour codes
636
637 int CoupSUSY::idSdown(int iSdown) {
638   int id = 0;
639   int sgn = ( iSdown > 0 ) ? 1 : -1;
640   iSdown = abs(iSdown);
641   if      (iSdown ==  1) id =  1000001; 
642   else if (iSdown ==  2) id =  1000003;     
643   else if (iSdown ==  3) id =  1000005; 
644   else if (iSdown ==  4) id =  2000001; 
645   else if (iSdown ==  5) id =  2000003;     
646   else if (iSdown ==  6) id =  2000005; 
647   return sgn*id;
648 }
649
650 //--------------------------------------------------------------------------
651
652 // Function to return slepton flavour codes
653
654 int CoupSUSY::idSlep(int iSlep) {
655   int id = 0;
656   int sgn = ( iSlep > 0 ) ? 1 : -1;
657   iSlep = abs(iSlep);
658   if      (iSlep ==  1) id =  1000011; 
659   else if (iSlep ==  2) id =  1000013;     
660   else if (iSlep ==  3) id =  1000015; 
661   else if (iSlep ==  4) id =  2000011; 
662   else if (iSlep ==  5) id =  2000013;     
663   else if (iSlep ==  6) id =  2000015; 
664   return sgn*id;
665 }
666
667 //--------------------------------------------------------------------------
668
669 // Return a particle name, given the PDG code.
670
671 string CoupSUSY::getName(int pdgCode) {    
672
673   // Absolute value and corresponding SM code
674   int codeA = abs(pdgCode);
675   int idSM  = codeA % 1000000;
676
677   // Name
678   string name;
679
680   // Flag to indicate whether SLHA1 or SLHA2
681   bool slha1 = false;
682
683   // SM particles
684   if (codeA == idSM) {
685     // Neutrinos
686     if (!slhaPtr->upmns.exists()) slha1=true;
687     if (codeA == 12) name = (slha1) ? "nu_e" : "nu_1";
688     if (codeA == 14) name = (slha1) ? "nu_mu" : "nu_2";
689     if (codeA == 16) name = (slha1) ? "nu_tau" : "nu_3";
690   }
691
692   // Squarks
693   else if ( idSM <= 6) {
694     // up squarks
695     if (idSM % 2 == 0) {
696       // If SLHA1, return old PDG names
697       if (slhaPtr->stopmix.exists()) slha1 = true;
698       if (codeA == 1000002) name = (slha1) ? "~u_L" : "~u_1";
699       if (codeA == 1000004) name = (slha1) ? "~c_L" : "~u_2";
700       if (codeA == 1000006) name = (slha1) ? "~t_1" : "~u_3";
701       if (codeA == 2000002) name = (slha1) ? "~u_R" : "~u_4";
702       if (codeA == 2000004) name = (slha1) ? "~c_R" : "~u_5";
703       if (codeA == 2000006) name = (slha1) ? "~t_2" : "~u_6";
704     } 
705     // down squarks
706     else {
707       // If SLHA1, return old PDG names
708       if (slhaPtr->sbotmix.exists()) slha1 = true;
709       if (codeA == 1000001) name = (slha1) ? "~d_L" : "~d_1";
710       if (codeA == 1000003) name = (slha1) ? "~s_L" : "~d_2";
711       if (codeA == 1000005) name = (slha1) ? "~b_1" : "~d_3";
712       if (codeA == 2000001) name = (slha1) ? "~d_R" : "~d_4";
713       if (codeA == 2000003) name = (slha1) ? "~s_R" : "~d_5";
714       if (codeA == 2000005) name = (slha1) ? "~b_2" : "~d_6";
715     }
716     if (pdgCode < 0) name += "bar";
717   } 
718
719   // Sleptons
720   else if ( idSM <= 19 ) {
721
722     // Sneutrinos
723    if (idSM % 2 == 0) {
724       // If SLHA1, return old PDG names
725       if (slhaPtr->staumix.exists()) slha1 = true;
726       if (codeA == 1000012) name = (slha1) ? "~nu_eL" : "~nu_1";
727       if (codeA == 1000014) name = (slha1) ? "~nu_muL" : "~nu_2";
728       if (codeA == 1000016) name = (slha1) ? "~nu_tauL" : "~nu_3";
729       if (codeA == 2000012) name = (slha1) ? "~nu_eR" : "~nu_4";
730       if (codeA == 2000014) name = (slha1) ? "~nu_muR" : "~nu_5";
731       if (codeA == 2000016) name = (slha1) ? "~nu_tauR" : "~nu_6";
732       if (pdgCode < 0) name += "bar";
733     }
734     // charged sleptons
735     else {
736       // If SLHA1, return old PDG names
737       if (slhaPtr->staumix.exists()) slha1 = true;
738       if (codeA == 1000011) name = (slha1) ? "~e_L" : "~e_1";
739       if (codeA == 1000013) name = (slha1) ? "~mu_L" : "~e_2";
740       if (codeA == 1000015) name = (slha1) ? "~tau_1" : "~e_3";
741       if (codeA == 2000011) name = (slha1) ? "~e_R" : "~e_4";
742       if (codeA == 2000013) name = (slha1) ? "~mu_R" : "~e_5";
743       if (codeA == 2000015) name = (slha1) ? "~tau_2" : "~e_6";
744       if (pdgCode < 0) name += "-";
745       else name += "+";
746     }
747   }
748
749   else if (codeA == 1000021) name = "~g";
750   else if (codeA == 1000022) name = "~chi_10";  
751   else if (codeA == 1000023) name = "~chi_20";
752   else if (codeA == 1000024) name = (pdgCode > 0) ? "~chi_1+" : "~chi_1-";
753   else if (codeA == 1000025) name = "~chi_30";
754   else if (codeA == 1000035) name = "~chi_40";
755   else if (codeA == 1000037) name = (pdgCode > 0) ? "~chi_2+" : "~chi_2-";
756
757   return name;
758
759 }
760   
761 //==========================================================================
762
763 } // end namespace Pythia8
764
765