source: HiSusy/trunk/Pythia8/pythia8170/src/SusyCouplings.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 37.0 KB
Line 
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
12namespace 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
30void 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
1028int 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
1043int 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
1057int 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
1075int 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
1093int 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
1111string 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
1205int 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
1221int 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
Note: See TracBrowser for help on using the repository browser.