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