StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SLHAinterface.cc
1 // SLHAinterface.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 #include "Pythia8/SLHAinterface.h"
8 
9 namespace Pythia8 {
10 
11 //==========================================================================
12 
13 // The SLHAinterface class.
14 
15 //--------------------------------------------------------------------------
16 
17 // Initialize and switch to SUSY couplings if reading SLHA spectrum.
18 
19 void SLHAinterface::init( Settings& settings, Rndm* rndmPtr,
20  Couplings* couplingsPtrIn, ParticleData* particleDataPtr,
21  bool& useSLHAcouplings) {
22 
23  // Initialize SLHA couplingsPtr to PYTHIA one by default
24  couplingsPtr = couplingsPtrIn;
25  useSLHAcouplings = false;
26 
27  // Check if SUSY couplings need to be read in
28  if( !initSLHA(settings, particleDataPtr))
29  infoPtr->errorMsg("Error in SLHAinterface::init: "
30  "Could not read SLHA file");
31 
32  // SLHA sets isSUSY flag to tell us if there was an SLHA SUSY spectrum
33  if (couplingsPtr->isSUSY) {
34  // Initialize the derived SUSY couplings class (SM first, then SUSY)
35  coupSUSY.init( settings, rndmPtr);
36  coupSUSY.initSUSY(&slha, infoPtr, particleDataPtr, &settings);
37  // Switch couplingsPtr to point to the derived class
38  // and tell PYTHIA to use it
39  couplingsPtr = (Couplings *) &coupSUSY;
40  useSLHAcouplings = true;
41  }
42 
43 }
44 
45 //--------------------------------------------------------------------------
46 
47 // Initialize SUSY Les Houches Accord data.
48 
49 bool SLHAinterface::initSLHA(Settings& settings,
50  ParticleData* particleDataPtr) {
51 
52  // Error and warning prefixes for this method
53  string errPref = "Error in SLHAinterface::initSLHA: ";
54  string warnPref = "Warning in SLHAinterface::initSLHA: ";
55  string infoPref = "Info from SLHAinterface::initSLHA: ";
56 
57  // Initial and settings values.
58  int ifailLHE = 1;
59  int ifailSpc = 1;
60  int readFrom = settings.mode("SLHA:readFrom");
61  string lhefFile = settings.word("Beams:LHEF");
62  string lhefHeader = settings.word("Beams:LHEFheader");
63  string slhaFile = settings.word("SLHA:file");
64  int verboseSLHA = settings.mode("SLHA:verbose");
65  bool slhaUseDec = settings.flag("SLHA:useDecayTable");
66 
67  // Set internal data members
68  meMode = settings.mode("SLHA:meMode");
69 
70  // No SUSY by default
71  couplingsPtr->isSUSY = false;
72 
73  // Option with no SLHA read-in at all.
74  if (readFrom == 0) return true;
75 
76  // First check LHEF header (if reading from LHEF)
77  if (readFrom == 1) {
78  if (lhefHeader != "void")
79  ifailLHE = slha.readFile(lhefHeader, verboseSLHA, slhaUseDec );
80  else if (lhefFile != "void")
81  ifailLHE = slha.readFile(lhefFile, verboseSLHA, slhaUseDec );
82  }
83 
84  // If LHEF read successful, everything needed should already be ready
85  if (ifailLHE == 0) {
86  ifailSpc = 0;
87  couplingsPtr->isSUSY = true;
88  // If no LHEF file or no SLHA info in header, read from SLHA:file
89  } else {
90  lhefFile = "void";
91  if ( settings.word("SLHA:file") == "none"
92  || settings.word("SLHA:file") == "void"
93  || settings.word("SLHA:file") == ""
94  || settings.word("SLHA:file") == " ") return true;
95  ifailSpc = slha.readFile(slhaFile,verboseSLHA, slhaUseDec);
96  }
97 
98  // In case of problems, print error and fail init.
99  if (ifailSpc != 0) {
100  infoPtr->errorMsg(errPref + "problem reading SLHA file", slhaFile);
101  return false;
102  } else {
103  couplingsPtr->isSUSY = true;
104  }
105 
106  // Check spectrum for consistency. Switch off SUSY if necessary.
107  ifailSpc = slha.checkSpectrum();
108 
109  // ifail >= 1 : no MODSEL found -> don't switch on SUSY
110  if (ifailSpc == 1) {
111  // no SUSY, but MASS ok
112  couplingsPtr->isSUSY = false;
113  infoPtr->errorMsg(infoPref +
114  "No MODSEL found, keeping internal SUSY switched off");
115  } else if (ifailSpc >= 2) {
116  // no SUSY, but problems
117  infoPtr->errorMsg(warnPref + "Problem with SLHA MASS or QNUMBERS.");
118  couplingsPtr->isSUSY = false;
119  }
120  // ifail = 0 : MODSEL found, spectrum OK
121  else if (ifailSpc == 0) {
122  // Print spectrum. Done.
123  slha.printSpectrum(0);
124  }
125  else if (ifailSpc < 0) {
126  infoPtr->errorMsg(warnPref + "Problem with SLHA spectrum.",
127  "\n Only using masses and switching off SUSY.");
128  settings.flag("SUSY:all", false);
129  couplingsPtr->isSUSY = false;
130  slha.printSpectrum(ifailSpc);
131  }
132 
133  // NMSSM spectrum (modify existing Higgs names and add particles)
134  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(3) >= 1 ) {
135  // Modify Higgs names
136  particleDataPtr->name(25,"H_1");
137  particleDataPtr->name(35,"H_2");
138  particleDataPtr->name(45,"H_3");
139  particleDataPtr->name(36,"A_1");
140  particleDataPtr->name(46,"A_2");
141  particleDataPtr->name(1000045,"~chi_50");
142  }
143 
144  // SLHA2 spectrum with flavour mixing (modify squark and/or slepton names)
145  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(6) >= 1 ) {
146  // Squark flavour violation
147  if ( (slha.modsel(6) == 1 || slha.modsel(6) >= 3)
148  && slha.usqmix.exists() && slha.dsqmix.exists() ) {
149  // Modify squark names
150  particleDataPtr->names(1000001,"~d1","~d1bar");
151  particleDataPtr->names(1000002,"~u1","~u1bar");
152  particleDataPtr->names(1000003,"~d2","~d2bar");
153  particleDataPtr->names(1000004,"~u2","~u2bar");
154  particleDataPtr->names(1000005,"~d3","~d3bar");
155  particleDataPtr->names(1000006,"~u3","~u3bar");
156  particleDataPtr->names(2000001,"~d4","~d4bar");
157  particleDataPtr->names(2000002,"~u4","~u4bar");
158  particleDataPtr->names(2000003,"~d5","~d5bar");
159  particleDataPtr->names(2000004,"~u5","~u5bar");
160  particleDataPtr->names(2000005,"~d6","~d6bar");
161  particleDataPtr->names(2000006,"~u6","~u6bar");
162  }
163  // Slepton flavour violation
164  if ( (slha.modsel(6) == 2 || slha.modsel(6) >= 3)
165  && slha.selmix.exists()) {
166  // Modify slepton names
167  particleDataPtr->names(1000011,"~e1-","~e1+");
168  particleDataPtr->names(1000013,"~e2-","~e2+");
169  particleDataPtr->names(1000015,"~e3-","~e3+");
170  particleDataPtr->names(2000011,"~e4-","~e4+");
171  particleDataPtr->names(2000013,"~e5-","~e5+");
172  particleDataPtr->names(2000015,"~e6-","~e6+");
173  }
174  // Sneutrino flavour violation
175  if ( (slha.modsel(6) == 2 || slha.modsel(6) >= 3)
176  && slha.snumix.exists()) {
177  // Modify sneutrino names
178  particleDataPtr->names(1000012,"~nu_1","~nu_1bar");
179  particleDataPtr->names(1000014,"~nu_2","~nu_2bar");
180  particleDataPtr->names(1000016,"~nu_3","~nu_3bar");
181  }
182  // Optionally allow for separate scalar and pseudoscalar sneutrinos
183  if ( slha.snsmix.exists() && slha.snamix.exists() ) {
184  // Scalar sneutrinos
185  particleDataPtr->names(1000012,"~nu_S1","~nu_S1bar");
186  particleDataPtr->names(1000014,"~nu_S2","~nu_S2bar");
187  particleDataPtr->names(1000016,"~nu_S3","~nu_S3bar");
188  // Add the pseudoscalar sneutrinos
189  particleDataPtr->addParticle(1000017, "~nu_A1", "~nu_A1bar",1, 0., 0);
190  particleDataPtr->addParticle(1000018, "~nu_A2", "~nu_A2bar",1, 0., 0);
191  particleDataPtr->addParticle(1000019, "~nu_A3", "~nu_A3bar",1, 0., 0);
192  }
193  }
194 
195  // SLHA2 spectrum with RPV
196  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(4) >= 1 ) {
197  if ( slha.rvnmix.exists() ) {
198  // Neutralinos -> neutrinos
199  particleDataPtr->names(12,"nu_1","nu_1bar");
200  particleDataPtr->names(14,"nu_2","nu_2bar");
201  particleDataPtr->names(16,"nu_3","nu_3bar");
202  particleDataPtr->names(1000022,"nu_4","nu_4bar");
203  particleDataPtr->names(1000023,"nu_5","nu_5bar");
204  particleDataPtr->names(1000025,"nu_6","nu_6bar");
205  particleDataPtr->names(1000035,"nu_7","nu_7bar");
206  }
207  if ( slha.rvumix.exists() && slha.rvvmix.exists() ) {
208  // Charginos -> charged leptons (note sign convention)
209  particleDataPtr->names(11,"e_1-","e_1+");
210  particleDataPtr->names(13,"e_2-","e_2+");
211  particleDataPtr->names(15,"e_3-","e_3+");
212  particleDataPtr->names(1000024,"e_4+","e_4-");
213  particleDataPtr->names(1000037,"e_5+","e_5-");
214  }
215  if ( slha.rvhmix.exists() ) {
216  // Sneutrinos -> higgses
217  particleDataPtr->name(25,"H_1");
218  particleDataPtr->name(35,"H_2");
219  particleDataPtr->name(1000012,"H_3");
220  particleDataPtr->name(1000014,"H_4");
221  particleDataPtr->name(1000016,"H_5");
222  }
223  if ( slha.rvamix.exists() ) {
224  // Sneutrinos -> higgses
225  particleDataPtr->name(36,"A_1");
226  particleDataPtr->name(1000017,"A_2");
227  particleDataPtr->name(1000018,"A_3");
228  particleDataPtr->name(1000019,"A_4");
229  }
230  if ( slha.rvlmix.exists() ) {
231  // sleptons -> charged higgses (note sign convention)
232  particleDataPtr->names(37,"H_1+","H_1-");
233  particleDataPtr->names(1000011,"H_2-","H_2+");
234  particleDataPtr->names(1000013,"H_3-","H_3+");
235  particleDataPtr->names(1000015,"H_4-","H_4+");
236  particleDataPtr->names(2000011,"H_5-","H_5+");
237  particleDataPtr->names(2000013,"H_6-","H_6+");
238  particleDataPtr->names(2000015,"H_7-","H_7+");
239  }
240  }
241 
242  // SLHA2 spectrum with CPV
243  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(5) >= 1 ) {
244  // no scalar/pseudoscalar distinction
245  particleDataPtr->name(25,"H_1");
246  particleDataPtr->name(35,"H_2");
247  particleDataPtr->name(36,"H_3");
248  }
249 
250  // Import qnumbers
251  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.qnumbers.size() > 0) {
252  for (int iQnum=0; iQnum < int(slha.qnumbers.size()); iQnum++) {
253  // Always use positive id codes
254  int id = abs(slha.qnumbers[iQnum](0));
255  ostringstream idCode;
256  idCode << id;
257  if (particleDataPtr->isParticle(id)) {
258  infoPtr->errorMsg(warnPref + "ignoring QNUMBERS", "for id = "
259  + idCode.str() + " (already exists)", true);
260  } else {
261  int qEM3 = slha.qnumbers[iQnum](1);
262  int nSpins = slha.qnumbers[iQnum](2);
263  int colRep = slha.qnumbers[iQnum](3);
264  int hasAnti = slha.qnumbers[iQnum](4);
265  // Translate colRep to PYTHIA colType
266  int colType = 0;
267  if (colRep == 3) colType = 1;
268  else if (colRep == -3) colType = -1;
269  else if (colRep == 8) colType = 2;
270  else if (colRep == 6) colType = 3;
271  else if (colRep == -6) colType = -3;
272  // Default name: PDG code
273  string name, antiName;
274  ostringstream idStream;
275  idStream<<id;
276  name = idStream.str();
277  antiName = "-"+name;
278  if (iQnum < int(slha.qnumbersName.size())) {
279  name = slha.qnumbersName[iQnum];
280  antiName = slha.qnumbersAntiName[iQnum];
281  if (antiName == "") {
282  if (name.find("+") != string::npos) {
283  antiName = name;
284  antiName.replace(antiName.find("+"),1,"-");
285  } else if (name.find("-") != string::npos) {
286  antiName = name;
287  antiName.replace(antiName.find("-"),1,"+");
288  } else {
289  antiName = name+"bar";
290  }
291  }
292  }
293  if ( hasAnti == 0) {
294  antiName = "";
295  particleDataPtr->addParticle(id, name, nSpins, qEM3, colType);
296  } else {
297  particleDataPtr->addParticle(id, name, antiName, nSpins, qEM3,
298  colType);
299  }
300  }
301  }
302  }
303 
304  // Import mass spectrum.
305  bool keepSM = settings.flag("SLHA:keepSM");
306  double minMassSM = settings.parm("SLHA:minMassSM");
307  bool allowUserOverride = settings.flag("SLHA:allowUserOverride");
308  vector<int> idModified;
309  if (ifailSpc == 1 || ifailSpc == 0) {
310 
311  // Start at beginning of mass array
312  int id = slha.mass.first();
313  // Loop through to update particle data.
314  for (int i = 1; i <= slha.mass.size() ; i++) {
315  // Step to next mass entry
316  if (i>1) id = slha.mass.next();
317  ostringstream idCode;
318  idCode << id;
319  double mass = abs(slha.mass(id));
320 
321  // Ignore masses for known SM particles or particles with
322  // default masses < minMassSM; overwrite masses for rest.
323  if (keepSM && (id < 25 || (id > 80 && id < 1000000))) ;
324  else if (id < 1000000 && particleDataPtr->m0(id) < minMassSM) {
325  cout<<" id = "<<id<<" m0 = "<<particleDataPtr->m0(id)<<" minMassSM = "
326  <<minMassSM<<endl;
327  infoPtr->errorMsg(warnPref + "ignoring MASS entry", "for id = "
328  + idCode.str() + " (m0 < SLHA:minMassSM)", true);
329  }
330 
331  // Also ignore SLHA mass values if user has already set
332  // a different value and is allowed to override them.
333  else if (allowUserOverride && particleDataPtr->hasChanged(id)) {
334  ostringstream mValue;
335  mValue << particleDataPtr->m0(id);
336  infoPtr->errorMsg(warnPref + "keeping user mass",
337  "for id = " + idCode.str() + ", m0 = " + mValue.str(), true);
338  idModified.push_back(id);
339  }
340  else {
341  particleDataPtr->m0(id,mass);
342  idModified.push_back(id);
343  }
344  };
345 
346  }
347 
348  // Update decay data.
349  for (int iTable=0; iTable < int(slha.decays.size()); iTable++) {
350 
351  // Pointer to this SLHA table
352  LHdecayTable* slhaTable=&(slha.decays[iTable]);
353 
354  // Extract ID and create pointer to corresponding particle data object
355  int idRes = slhaTable->getId();
356  ostringstream idCode;
357  idCode << idRes;
358  ParticleDataEntry* particlePtr
359  = particleDataPtr->particleDataEntryPtr(idRes);
360 
361  // Ignore decay channels for known SM particles or particles with
362  // default masses < minMassSM; overwrite masses for rest.
363  if (keepSM && (idRes < 25 || (idRes > 80 && idRes < 1000000))) continue;
364  else if (idRes < 1000000 && particleDataPtr->m0(idRes) < minMassSM
365  && !particleDataPtr->hasChanged(idRes) ) {
366  infoPtr->errorMsg(warnPref + "ignoring DECAY table", "for id = "
367  + idCode.str() + " (m0 < SLHA:minMassSM)", true);
368  continue;
369  }
370 
371  // Verbose output. Let user know we are using these tables.
372  if (verboseSLHA <= 1)
373  infoPtr->errorMsg(infoPref + "importing SLHA decay table(s)","");
374  else
375  infoPtr->errorMsg(infoPref + "importing SLHA decay table","for id = "
376  +idCode.str(),true);
377 
378  // Extract and store total width (absolute value, neg -> switch off)
379  double widRes = abs(slhaTable->getWidth());
380  double pythiaMinWidth = settings.parm("ResonanceWidths:minWidth");
381  if (widRes > 0. && widRes < pythiaMinWidth) {
382  infoPtr->errorMsg(warnPref + "forcing width = 0 ","for id = "
383  + idCode.str() + " (width < ResonanceWidths:minWidth)" , true);
384  widRes = 0.0;
385  }
386  particlePtr->setMWidth(widRes);
387 
388  // Set lifetime in mm for displaced vertex calculations
389  // (convert GeV^-1 to mm)
390  if (widRes > 0.) {
391  double decayLength = 1.97e-13/widRes;
392  particlePtr->setTau0(decayLength);
393 
394  // Reset decay table of the particle. Allow decays, treat as resonance.
395  if (slhaTable->size() > 0) {
396  particlePtr->clearChannels();
397  particleDataPtr->mayDecay(idRes,true);
398  particleDataPtr->isResonance(idRes,true);
399  } else {
400  infoPtr->errorMsg(warnPref + "empty DECAY table ","for id = "
401  + idCode.str() + " (total width provided but no branching"
402  + " fractions)", true);
403  }
404  }
405  // Reset to stable if width <= 0.0
406  else {
407  particlePtr->clearChannels();
408  particleDataPtr->mayDecay(idRes,false);
409  particleDataPtr->isResonance(idRes,false);
410  }
411 
412  // Set initial minimum mass.
413  double brWTsum = 0.;
414  double massWTsum = 0.;
415 
416  // Loop over SLHA channels, import into Pythia, treating channels
417  // with negative branching fractions as having the equivalent positive
418  // branching fraction, but being switched off for this run
419  for (int iChannel=0 ; iChannel<slhaTable->size(); iChannel++) {
420  LHdecayChannel slhaChannel = slhaTable->getChannel(iChannel);
421  double brat = slhaChannel.getBrat();
422  vector<int> idDa = slhaChannel.getIdDa();
423  if (idDa.size() >= 9) {
424  infoPtr->errorMsg(errPref + "max number of DECAY products is 8");
425  } else if (idDa.size() <= 1) {
426  infoPtr->errorMsg(errPref + "min number of DECAY products is 2");
427  }
428  else {
429  int onMode = 1;
430  if (brat < 0.0) onMode = 0;
431  int meModeNow = meMode;
432 
433  // Check phase space, including margin ~ sqrt(sum(widths^2))
434  double massSum(0.);
435  double widSqSum = pow2(widRes);
436  int nDa = idDa.size();
437  for (int jDa=0; jDa<nDa; ++jDa) {
438  massSum += particleDataPtr->m0( idDa[jDa] );
439  widSqSum += pow2(particleDataPtr->mWidth( idDa[jDa] ));
440  }
441  double deltaM = particleDataPtr->m0(idRes) - massSum;
442  // Negative mass difference: intrinsically off shell
443  if (onMode == 1 && brat > 0.0 && deltaM < 0.) {
444  // String containing decay name
445  ostringstream errCode;
446  errCode << idRes <<" ->";
447  for (int jDa=0; jDa<nDa; ++jDa) errCode<<" "<<idDa[jDa];
448  // Could mass fluctuations at all give the needed deltaM ?
449  if (abs(deltaM) > 100. * sqrt(widSqSum)) {
450  infoPtr->errorMsg(warnPref + "switched off DECAY mode",
451  ": " + errCode.str()+" (too far off shell)",true);
452  onMode = 0;
453  }
454  // If ~ OK within fluctuations
455  else {
456  // Ignore user-selected meMode
457  if (meModeNow != 100) {
458  infoPtr->errorMsg(warnPref + "adding off shell DECAY mode",
459  ": "+errCode.str()+" (forced meMode = 100)",true);
460  meModeNow = 100;
461  } else {
462  infoPtr->errorMsg(warnPref + "adding off shell DECAY mode",
463  errCode.str(), true);
464  }
465  }
466  }
467  // Branching-ratio-weighted average mass in decay.
468  brWTsum += abs(brat);
469  massWTsum += abs(brat) * massSum;
470 
471  // Add channel
472  int id0 = idDa[0];
473  int id1 = idDa[1];
474  int id2 = (idDa.size() >= 3) ? idDa[2] : 0;
475  int id3 = (idDa.size() >= 4) ? idDa[3] : 0;
476  int id4 = (idDa.size() >= 5) ? idDa[4] : 0;
477  int id5 = (idDa.size() >= 6) ? idDa[5] : 0;
478  int id6 = (idDa.size() >= 7) ? idDa[6] : 0;
479  int id7 = (idDa.size() >= 8) ? idDa[7] : 0;
480  particlePtr->addChannel(onMode,abs(brat),meModeNow,
481  id0,id1,id2,id3,id4,id5,id6,id7);
482 
483  }
484  }
485 
486  // Set minimal mass, but always below nominal one.
487  if (slhaTable->size() > 0) {
488  double massAvg = massWTsum / brWTsum;
489  double massMin = min( massAvg, particlePtr->m0()) ;
490  particlePtr->setMMin(massMin);
491  }
492 
493  // Add to list of particles that have been modified
494  idModified.push_back(idRes);
495 
496  }
497 
498  // Sanity check of all decay tables with modified MASS or DECAY info
499  for (int iMod = 0; iMod < int(idModified.size()); ++iMod) {
500  int id = idModified[iMod];
501  ostringstream idCode;
502  idCode << id;
503  ParticleDataEntry* particlePtr
504  = particleDataPtr->particleDataEntryPtr(id);
505  double m0 = particlePtr->m0();
506  double wid = particlePtr->mWidth();
507  // Always set massless particles stable
508  if (m0 <= 0.0 && (wid > 0.0 || particlePtr->mayDecay())) {
509  infoPtr->errorMsg(warnPref + "massless particle forced stable"," id = "
510  + idCode.str(), true);
511  particlePtr->clearChannels();
512  particlePtr->setMWidth(0.0);
513  particlePtr->setMayDecay(false);
514  particleDataPtr->isResonance(id,false);
515  continue;
516  }
517  // Declare zero-width particles to be stable (for now)
518  if (wid == 0.0 && particlePtr->mayDecay()) {
519  particlePtr->setMayDecay(false);
520  continue;
521  }
522  // Check at least one on-shell channel is available
523  double mSumMin = 10. * m0;
524  int nChannels = particlePtr->sizeChannels();
525  if (nChannels >= 1) {
526  for (int iChannel=0; iChannel<nChannels; ++iChannel) {
527  DecayChannel channel = particlePtr->channel(iChannel);
528  if (channel.onMode() <= 0) continue;
529  int nProd = channel.multiplicity();
530  double mSum = 0.;
531  for (int iDa = 0; iDa < nProd; ++iDa) {
532  int idDa = channel.product(iDa);
533  mSum += particleDataPtr->m0(idDa);
534  }
535  mSumMin = min(mSumMin, mSum);
536  }
537  // Require at least one on-shell channel
538  if (mSumMin > m0) {
539  infoPtr->errorMsg(warnPref + "particle forced stable"," id = "
540  + idCode.str() + " (no on-shell decay channels)", true);
541  particlePtr->setMWidth(0.0);
542  particlePtr->setMayDecay(false);
543  continue;
544  }
545  else {
546  // mMin: lower cutoff on Breit-Wigner: default is mMin = m0 - 5*Gamma
547  // (User is allowed to specify a lower value if desired.)
548  // Increase minimum if needed to ensure at least one channel on shell
549  double mMin = min(particlePtr->mMin(), max(0.0,m0 - 5.*wid));
550  mMin = max(mSumMin,mMin);
551  particlePtr->setMMin(mMin);
552  }
553  }
554  }
555 
556  return true;
557 
558 }
559 
560 //--------------------------------------------------------------------------
561 
562 // Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
563 // E.g., to make sure that there are no important unfilled entries
564 
565 void SLHAinterface::pythia2slha(ParticleData* particleDataPtr) {
566 
567  // Initialize block SMINPUTS.
568  string blockName = "sminputs";
569  double mZ = particleDataPtr->m0(23);
570  slha.set(blockName,1,1.0/couplingsPtr->alphaEM(pow2(mZ)));
571  slha.set(blockName,2,couplingsPtr->GF());
572  slha.set(blockName,3,couplingsPtr->alphaS(pow2(mZ)));
573  slha.set(blockName,4,mZ);
574  // b mass (should be running mass, here pole for time being)
575  slha.set(blockName,5,particleDataPtr->m0(5));
576  slha.set(blockName,6,particleDataPtr->m0(6));
577  slha.set(blockName,7,particleDataPtr->m0(15));
578  slha.set(blockName,8,particleDataPtr->m0(16));
579  slha.set(blockName,11,particleDataPtr->m0(11));
580  slha.set(blockName,12,particleDataPtr->m0(12));
581  slha.set(blockName,13,particleDataPtr->m0(13));
582  slha.set(blockName,14,particleDataPtr->m0(14));
583  // Force 3 lightest quarks massless
584  slha.set(blockName,21,double(0.0));
585  slha.set(blockName,22,double(0.0));
586  slha.set(blockName,23,double(0.0));
587  // c mass (should be running mass, here pole for time being)
588  slha.set(blockName,24,particleDataPtr->m0(4));
589 
590  // Initialize block MASS.
591  blockName = "mass";
592  int id = 1;
593  int count = 0;
594  while (particleDataPtr->nextId(id) > id) {
595  slha.set(blockName,id,particleDataPtr->m0(id));
596  id = particleDataPtr->nextId(id);
597  ++count;
598  if (count > 10000) {
599  infoPtr->errorMsg("Error in SLHAinterface::pythia2slha(): "
600  "encountered infinite loop when saving mass block");
601  break;
602  }
603  }
604 
605 }
606 
607 //==========================================================================
608 
609 } // end namespace Pythia8
610 
611 
612 
Generated by   doxygen 1.8.5