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) 2020 Torbjorn Sjostrand.
3 // Main authors of this file: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL v2 or later, 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( bool& useSLHAcouplings,
20  stringstream& particleDataBuffer) {
21 
22  // By default no SLHA couplings.
23  useSLHAcouplings = false;
24 
25  // Check if SUSY couplings need to be read in
26  if( !initSLHA())
27  infoPtr->errorMsg("Error in SLHAinterface::init: "
28  "Could not read SLHA file");
29 
30  // Reset any particle-related user settings.
31  string line;
32  string warnPref = "Warning in SLHAinterface::init: ";
33  while (getline(particleDataBuffer, line)
34  && settingsPtr->flag("SLHA:allowUserOverride")) {
35  bool pass = particleDataPtr->readString(line, true);
36  if (!pass) infoPtr->errorMsg(warnPref + "Unable to process line "
37  + line);
38  else infoPtr->errorMsg(warnPref + "Overwriting SLHA by " + line);
39  }
40 
41  // SLHA sets isSUSY flag to tell us if there was an SLHA SUSY spectrum
42  if (coupSUSYPtr->isSUSY) {
43  // Initialize the derived SUSY couplings class.
44  coupSUSYPtr->initSUSY(&slha, infoPtr);
45  useSLHAcouplings = true;
46  }
47 
48 }
49 
50 //--------------------------------------------------------------------------
51 
52 // Initialize SUSY Les Houches Accord data.
53 
54 bool SLHAinterface::initSLHA() {
55 
56  // Error and warning prefixes for this method
57  string errPref = "Error in SLHAinterface::initSLHA: ";
58  string warnPref = "Warning in SLHAinterface::initSLHA: ";
59  string infoPref = "Info from SLHAinterface::initSLHA: ";
60 
61  // Initial and settings values.
62  int ifailLHE = 1;
63  int ifailSpc = 1;
64  int readFrom = settingsPtr->mode("SLHA:readFrom");
65  string lhefFile = settingsPtr->word("Beams:LHEF");
66  string lhefHeader = settingsPtr->word("Beams:LHEFheader");
67  string slhaFile = settingsPtr->word("SLHA:file");
68  int verboseSLHA = settingsPtr->mode("SLHA:verbose");
69  bool slhaUseDec = settingsPtr->flag("SLHA:useDecayTable");
70  bool noSLHAFile = ( slhaFile == "none" || slhaFile == "void"
71  || slhaFile == "" || slhaFile == " " );
72 
73  // Set internal data members
74  meMode = settingsPtr->mode("SLHA:meMode");
75 
76  // No SUSY by default
77  coupSUSYPtr->isSUSY = false;
78 
79  // Option with no SLHA read-in at all.
80  if (readFrom == 0) return true;
81 
82  // First check LHEF header (if reading from LHEF)
83  if (readFrom == 1) {
84  // Check if there is a <slha> tag in the LHEF header
85  // Note: if the <slha> tag is NOT inside the <header>, it will be ignored.
86  string slhaInHeader( infoPtr->header("slha") );
87  if (slhaInHeader == "" && noSLHAFile) return true;
88  // If there is an <slha> tag, read file.
89  if (lhefHeader != "void")
90  ifailLHE = slha.readFile(lhefHeader, verboseSLHA, slhaUseDec );
91  else if (lhefFile != "void")
92  ifailLHE = slha.readFile(lhefFile, verboseSLHA, slhaUseDec );
93  else if (noSLHAFile) {
94  istringstream slhaInHeaderStream(slhaInHeader);
95  ifailLHE = slha.readFile(slhaInHeaderStream, verboseSLHA, slhaUseDec );
96  }
97  }
98 
99  // If LHEF read successful, everything needed should already be ready
100  if (ifailLHE == 0) {
101  ifailSpc = 0;
102  coupSUSYPtr->isSUSY = true;
103  // If no LHEF file or no SLHA info in header, read from SLHA:file
104  } else {
105  lhefFile = "void";
106  if ( noSLHAFile ) return true;
107  ifailSpc = slha.readFile(slhaFile,verboseSLHA, slhaUseDec);
108  }
109 
110  // In case of problems, print error and fail init.
111  if (ifailSpc != 0) {
112  infoPtr->errorMsg(errPref + "problem reading SLHA file", slhaFile);
113  return false;
114  } else {
115  coupSUSYPtr->isSUSY = true;
116  }
117 
118  // Check spectrum for consistency. Switch off SUSY if necessary.
119  ifailSpc = slha.checkSpectrum();
120 
121  // ifail >= 1 : no MODSEL found -> don't switch on SUSY
122  if (ifailSpc == 1) {
123  // no SUSY, but MASS ok
124  coupSUSYPtr->isSUSY = false;
125  infoPtr->errorMsg(infoPref +
126  "No MODSEL found, keeping internal SUSY switched off");
127  } else if (ifailSpc >= 2) {
128  // no SUSY, but problems
129  infoPtr->errorMsg(warnPref +
130  "No MASS or MODSEL blocks found.");
131  coupSUSYPtr->isSUSY = false;
132  }
133  // ifail = 0 : MODSEL found, spectrum OK
134  else if (ifailSpc == 0) {
135  // Print spectrum. Done.
136  slha.listSpectrum(0);
137  }
138  else if (ifailSpc < 0) {
139  infoPtr->errorMsg(warnPref + "Problem with SLHA spectrum.",
140  "\n Only using masses and switching off SUSY.");
141  settingsPtr->flag("SUSY:all", false);
142  coupSUSYPtr->isSUSY = false;
143  slha.listSpectrum(ifailSpc);
144  }
145 
146  // SLHA1 : SLHA2 compatibility
147  // Check whether scalar particle masses are ordered
148  bool isOrderedQ = true;
149  bool isOrderedL = true;
150  int idSdown[6]={1000001,1000003,1000005,2000001,2000003,2000005};
151  int idSup[6]={1000002,1000004,1000006,2000002,2000004,2000006};
152  int idSlep[6]={1000011,1000013,1000015,2000011,2000013,2000015};
153  for (int j=0;j<=4;j++) {
154  if (slha.mass(idSlep[j+1]) < slha.mass(idSlep[j]))
155  isOrderedL = false;
156  if (slha.mass(idSup[j+1]) < slha.mass(idSup[j]))
157  isOrderedQ = false;
158  if (slha.mass(idSdown[j+1]) < slha.mass(idSdown[j]))
159  isOrderedQ = false;
160  }
161 
162  // If ordered, change sparticle labels to mass-ordered enumeration
163  for (int i=1;i<=6;i++) {
164  ostringstream indx;
165  indx << i;
166  string uName = "~u_"+indx.str();
167  string dName = "~d_"+indx.str();
168  string lName = "~e_"+indx.str();
169  if (isOrderedQ) {
170  particleDataPtr->names(idSup[i-1],uName,uName+"bar");
171  particleDataPtr->names(idSdown[i-1],dName,dName+"bar");
172  }
173  if (isOrderedL) particleDataPtr->names(idSlep[i-1],lName+"-",lName+"+");
174  }
175 
176  // NMSSM spectrum (modify existing Higgs names and add particles)
177  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(3) >= 1 ) {
178  // Modify Higgs names
179  particleDataPtr->name(25,"H_1");
180  particleDataPtr->name(35,"H_2");
181  particleDataPtr->name(45,"H_3");
182  particleDataPtr->name(36,"A_1");
183  particleDataPtr->name(46,"A_2");
184  particleDataPtr->name(1000045,"~chi_50");
185  }
186 
187  // SLHA2 spectrum with flavour mixing (modify squark and/or slepton names)
188  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(6) >= 1 ) {
189  // Squark flavour violation
190  if ( (slha.modsel(6) == 1 || slha.modsel(6) >= 3)
191  && slha.usqmix.exists() && slha.dsqmix.exists() ) {
192  // Modify squark names
193  particleDataPtr->names(1000001,"~d_1","~d_1bar");
194  particleDataPtr->names(1000002,"~u_1","~u_1bar");
195  particleDataPtr->names(1000003,"~d_2","~d_2bar");
196  particleDataPtr->names(1000004,"~u_2","~u_2bar");
197  particleDataPtr->names(1000005,"~d_3","~d_3bar");
198  particleDataPtr->names(1000006,"~u_3","~u_3bar");
199  particleDataPtr->names(2000001,"~d_4","~d_4bar");
200  particleDataPtr->names(2000002,"~u_4","~u_4bar");
201  particleDataPtr->names(2000003,"~d_5","~d_5bar");
202  particleDataPtr->names(2000004,"~u_5","~u_5bar");
203  particleDataPtr->names(2000005,"~d_6","~d_6bar");
204  particleDataPtr->names(2000006,"~u_6","~u_6bar");
205  }
206  // Slepton flavour violation
207  if ( (slha.modsel(6) == 2 || slha.modsel(6) >= 3)
208  && slha.selmix.exists()) {
209  // Modify slepton names
210  particleDataPtr->names(1000011,"~e_1-","~e_1+");
211  particleDataPtr->names(1000013,"~e_2-","~e_2+");
212  particleDataPtr->names(1000015,"~e_3-","~e_3+");
213  particleDataPtr->names(2000011,"~e_4-","~e_4+");
214  particleDataPtr->names(2000013,"~e_5-","~e_5+");
215  particleDataPtr->names(2000015,"~e_6-","~e_6+");
216  }
217  // Neutrino flavour violation
218  if ( (slha.modsel(6) == 2 || slha.modsel(6) >= 3)
219  && slha.upmns.exists()) {
220  // Modify neutrino names (note that SM processes may not use UPMNS)
221  particleDataPtr->names(12,"nu_1","nu_1bar");
222  particleDataPtr->names(14,"nu_2","nu_2bar");
223  particleDataPtr->names(16,"nu_3","nu_3bar");
224  }
225  // Sneutrino flavour violation
226  if ( (slha.modsel(6) == 2 || slha.modsel(6) >= 3)
227  && slha.snumix.exists()) {
228  // Modify sneutrino names
229  particleDataPtr->names(1000012,"~nu_1","~nu_1bar");
230  particleDataPtr->names(1000014,"~nu_2","~nu_2bar");
231  particleDataPtr->names(1000016,"~nu_3","~nu_3bar");
232  }
233  // Optionally allow for separate scalar and pseudoscalar sneutrinos
234  if ( slha.snsmix.exists() && slha.snamix.exists() ) {
235  // Scalar sneutrinos
236  particleDataPtr->names(1000012,"~nu_S1","~nu_S1bar");
237  particleDataPtr->names(1000014,"~nu_S2","~nu_S2bar");
238  particleDataPtr->names(1000016,"~nu_S3","~nu_S3bar");
239  // Add the pseudoscalar sneutrinos
240  particleDataPtr->addParticle(1000017, "~nu_A1", "~nu_A1bar",1, 0., 0);
241  particleDataPtr->addParticle(1000018, "~nu_A2", "~nu_A2bar",1, 0., 0);
242  particleDataPtr->addParticle(1000019, "~nu_A3", "~nu_A3bar",1, 0., 0);
243  }
244  }
245 
246  // SLHA2 spectrum with RPV
247  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(4) >= 1 ) {
248  if ( slha.rvnmix.exists() ) {
249  // Neutralinos -> neutrinos
250  // Maintain R-conserving names since mass-ordering unlikely to change.
251  particleDataPtr->names(12,"nu_1","nu_1bar");
252  particleDataPtr->names(14,"nu_2","nu_2bar");
253  particleDataPtr->names(16,"nu_3","nu_3bar");
254  particleDataPtr->name(1000022,"~chi_10");
255  particleDataPtr->name(1000023,"~chi_20");
256  particleDataPtr->name(1000025,"~chi_30");
257  particleDataPtr->name(1000035,"~chi_40");
258  }
259  if ( slha.rvumix.exists() && slha.rvvmix.exists() ) {
260  // Charginos -> charged leptons (note sign convention)
261  // Maintain R-conserving names since mass-ordering unlikely to change.
262  particleDataPtr->names(11,"e-","e+");
263  particleDataPtr->names(13,"mu-","mu+");
264  particleDataPtr->names(15,"tau-","tau+");
265  particleDataPtr->names(1000024,"~chi_1+","~chi_1-");
266  particleDataPtr->names(1000037,"~chi_2+","~chi_2-");
267  }
268  if ( slha.rvhmix.exists() ) {
269  // Sneutrinos -> higgses (general mass-ordered names)
270  particleDataPtr->name(25,"H_10");
271  particleDataPtr->name(35,"H_20");
272  particleDataPtr->names(1000012,"H_30","H_30");
273  particleDataPtr->names(1000014,"H_40","H_40");
274  particleDataPtr->names(1000016,"H_50","H_50");
275  }
276  if ( slha.rvamix.exists() ) {
277  // Sneutrinos -> higgses (general mass-ordered names)
278  particleDataPtr->name(36,"A_10");
279  particleDataPtr->names(1000017,"A_20","A_20");
280  particleDataPtr->names(1000018,"A_30","A_30");
281  particleDataPtr->names(1000019,"A_40","A_40");
282  }
283  if ( slha.rvlmix.exists() ) {
284  // sleptons -> charged higgses (note sign convention)
285  particleDataPtr->names(37,"H_1+","H_1-");
286  particleDataPtr->names(1000011,"H_2-","H_2+");
287  particleDataPtr->names(1000013,"H_3-","H_3+");
288  particleDataPtr->names(1000015,"H_4-","H_4+");
289  particleDataPtr->names(2000011,"H_5-","H_5+");
290  particleDataPtr->names(2000013,"H_6-","H_6+");
291  particleDataPtr->names(2000015,"H_7-","H_7+");
292  }
293  }
294 
295  // SLHA2 spectrum with CPV
296  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.modsel(5) >= 1 ) {
297  // no scalar/pseudoscalar distinction
298  particleDataPtr->name(25,"H_10");
299  particleDataPtr->name(35,"H_20");
300  particleDataPtr->name(36,"H_30");
301  }
302 
303 
304  // Import qnumbers
305  vector<int> isQnumbers;
306  bool foundLowCode = false;
307  if ( (ifailSpc == 1 || ifailSpc == 0) && slha.qnumbers.size() > 0) {
308  for (int iQnum=0; iQnum < int(slha.qnumbers.size()); iQnum++) {
309  // Always use positive id codes
310  int id = abs(slha.qnumbers[iQnum](0));
311  ostringstream idCode;
312  idCode << id;
313  if (particleDataPtr->isParticle(id)) {
314  infoPtr->errorMsg(warnPref + "ignoring QNUMBERS", "for id = "
315  + idCode.str() + " (already exists)", true);
316  } else {
317  // Note: qnumbers entries stored internally as doubles to allow for
318  // millicharged particles. Round to nearest int when handing to PYTHIA.
319  // (So charge ~ 0 treated as charge = 0; good enough for showers.)
320  int qEM3 = lrint(slha.qnumbers[iQnum](1));
321  int nSpins = lrint(slha.qnumbers[iQnum](2));
322  int colRep = lrint(slha.qnumbers[iQnum](3));
323  int hasAnti = lrint(slha.qnumbers[iQnum](4));
324  // Translate colRep to PYTHIA colType
325  int colType = 0;
326  if (colRep == 3) colType = 1;
327  else if (colRep == -3) colType = -1;
328  else if (colRep == 8) colType = 2;
329  else if (colRep == 6) colType = 3;
330  else if (colRep == -6) colType = -3;
331  // Default name: PDG code
332  string name, antiName;
333  ostringstream idStream;
334  idStream<<id;
335  name = idStream.str();
336  antiName = "-"+name;
337  if (iQnum < int(slha.qnumbersName.size())) {
338  name = slha.qnumbersName[iQnum];
339  antiName = slha.qnumbersAntiName[iQnum];
340  if (antiName == "") {
341  if (name.find("+") != string::npos) {
342  antiName = name;
343  antiName.replace(antiName.find("+"),1,"-");
344  } else if (name.find("-") != string::npos) {
345  antiName = name;
346  antiName.replace(antiName.find("-"),1,"+");
347  } else {
348  antiName = name+"bar";
349  }
350  }
351  }
352  if ( hasAnti == 0) {
353  antiName = "";
354  particleDataPtr->addParticle(id, name, nSpins, qEM3, colType);
355  } else {
356  particleDataPtr->addParticle(id, name, antiName, nSpins, qEM3,
357  colType);
358  }
359  // Store list of particle codes added by QNUMBERS
360  isQnumbers.push_back(id);
361  if (id < 1000000) foundLowCode = true;
362  }
363  }
364  }
365  // Inform user that BSM particles should ideally be assigned id codes > 1M
366  if (foundLowCode)
367  infoPtr->errorMsg(warnPref
368  + "using QNUMBERS for id codes < 1000000 may clash with SM.");
369 
370  // Import mass spectrum.
371  double minMassSM = settingsPtr->parm("SLHA:minMassSM");
372  map<int,bool> idModified;
373  if (ifailSpc == 1 || ifailSpc == 0) {
374 
375  // Start at beginning of mass array
376  int id = slha.mass.first();
377  vector<int> ignoreMassKeepSM;
378  vector<int> ignoreMassM0;
379  vector<int> importMass;
380  // Loop through to update particle data.
381  for (int i = 1; i <= slha.mass.size() ; i++) {
382  // Step to next mass entry
383  if (i>1) id = slha.mass.next();
384  double mass = abs(slha.mass(id));
385 
386  // Check if this ID was added by qnumbers
387  bool isInternal = true;
388  for (unsigned int iq = 0; iq<isQnumbers.size(); ++iq)
389  if (id == isQnumbers[iq]) isInternal = false;
390 
391  // Ignore masses for known SM particles or particles with
392  // default masses < minMassSM; overwrite masses for rest.
393  // SM particles: (idRes < 25 || idRes > 80 && idRes < 1000000);
394  // Top and Higgs (25) can be overwritten if mass > minMassSM
395  // Extra higgses (26 - 40) & DM (51 - 60) not SM
396  // minMassSM : mass above which SM masses may be overwritten
397  bool isSM = id < 26 || ( id > 80 && id < 1000000);
398  if (isSM && isInternal) {
399  if (particleDataPtr->m0(id) < minMassSM) {
400  ignoreMassM0.push_back(id);
401  }
402  } else {
403  ParticleDataEntry* tmpPtr = particleDataPtr->findParticle(id);
404  if( tmpPtr == NULL ) {
405  ostringstream idCode;
406  idCode << id;
407  infoPtr->errorMsg(infoPref + " attempting to set properties",
408  "for unknown id = {" + idCode.str() + "}", true);
409  continue;
410  }
411  particleDataPtr->m0(id,mass);
412  idModified[id] = true;
413  importMass.push_back(id);
414  // If the mMin and mMax cutoffs on Breit-Wigner tails were not already
415  // set by user, set default bounds to at most m0 +- m0/2.
416  // Treat these values as new defaults: do not set hasChanged flags.
417  // Note: tighter bounds may apply if a width is given later; see below.
418  if (!particleDataPtr->hasChangedMMin(id))
419  particleDataPtr->findParticle(id)->setMMinNoChange( mass/2. );
420  if (!particleDataPtr->hasChangedMMax(id))
421  particleDataPtr->findParticle(id)->setMMaxNoChange( 3.*mass/2. );
422  }
423  };
424  // Give summary of any imported/ignored MASS entries, and state reason
425  if (importMass.size() >= 1) {
426  string idImport;
427  for (unsigned int i=0; i<importMass.size(); ++i) {
428  ostringstream idCode;
429  idCode << importMass[i];
430  if (i != 0) idImport +=",";
431  idImport += idCode.str();
432  }
433  infoPtr->errorMsg(infoPref + "importing MASS entries","for id = {"
434  + idImport + "}", true);
435  }
436  if (ignoreMassM0.size() >= 1) {
437  string idIgnore;
438  for (unsigned int i=0; i<ignoreMassM0.size(); ++i) {
439  ostringstream idCode;
440  idCode << ignoreMassM0[i];
441  if (i != 0) idIgnore +=",";
442  idIgnore += idCode.str();
443  }
444  infoPtr->errorMsg(warnPref + "ignoring MASS entries", "for id = {"
445  + idIgnore + "}" + " (m0 < SLHA:minMassSM)", true);
446  }
447  }
448 
449  // Update decay data.
450  vector<int> ignoreDecayKeepSM;
451  vector<int> ignoreDecayM0;
452  vector<int> ignoreDecayBR;
453  vector<int> importDecay;
454  for (int iTable=0; iTable < int(slha.decays.size()); iTable++) {
455 
456  // Pointer to this SLHA table
457  LHdecayTable* slhaTable=&(slha.decays[iTable]);
458 
459  // Extract ID and create pointer to corresponding particle data object
460  int idRes = slhaTable->getId();
461  ostringstream idCode;
462  idCode << idRes;
463  ParticleDataEntry* particlePtr
464  = particleDataPtr->particleDataEntryPtr(idRes);
465 
466  // Check if this ID was added by qnumbers
467  bool isInternal = true;
468  for (unsigned int iq = 0; iq<isQnumbers.size(); ++iq)
469  if (idRes == isQnumbers[iq]) isInternal = false;
470 
471 
472  // Ignore decay channels for known SM particles or particles with
473  // default masses < minMassSM; overwrite masses for rest.
474  // Let extra Higgses & Dark Matter sector be non-SM
475  bool isSM = idRes < 26 || ( idRes > 80 && idRes < 1000000);
476  if (isSM && isInternal) {
477  if(particleDataPtr->m0(idRes) < minMassSM) {
478  ignoreDecayM0.push_back(idRes);
479  continue;
480  }
481  }
482 
483  // Extract and store total width (absolute value, neg -> switch off)
484  double widRes = abs(slhaTable->getWidth());
485  double pythiaMinWidth = settingsPtr->parm("ResonanceWidths:minWidth");
486  if (widRes > 0. && widRes < pythiaMinWidth) {
487  infoPtr->errorMsg(warnPref + "forcing width = 0 ","for id = "
488  + idCode.str() + " (width < ResonanceWidths:minWidth)" , true);
489  widRes = 0.0;
490  }
491  particlePtr->setMWidth(widRes);
492  // If the mMin and mMax cutoffs on Breit-Wigner tails were not already
493  // set by user, set default values for them to 5*width though at most m0/2.
494  // Treat these values as defaults, ie do not set hasChanged flags.
495  // (After all channels have been read in, we also check that mMin is
496  // high enough to allow at least one channel to be on shell; see below.)
497  if (!particlePtr->hasChangedMMin()) {
498  double m0 = particlePtr->m0();
499  double mMin = m0 - min(5*widRes , m0/2.);
500  particlePtr->setMMinNoChange(mMin);
501  }
502  if (!particlePtr->hasChangedMMax()) {
503  double m0 = particlePtr->m0();
504  double mMax = m0 + min(5*widRes , m0/2.);
505  particlePtr->setMMaxNoChange(mMax);
506  }
507 
508  // Set lifetime for displaced vertex calculations (convert GeV^-1 to mm).
509  if (widRes > 0.) {
510  double decayLength = HBARC * FM2MM / widRes;
511  particlePtr->setTau0(decayLength);
512 
513  // Reset decay table of the particle. Allow decays, treat as resonance.
514  if (slhaTable->size() > 0) {
515  particlePtr->clearChannels();
516  particleDataPtr->mayDecay(idRes,true);
517  particleDataPtr->isResonance(idRes,true);
518  // Let user know we are using these tables.
519  importDecay.push_back(idRes);
520  } else
521  ignoreDecayBR.push_back(idRes);
522  }
523 
524  // Reset to stable if width <= 0.0
525  else {
526  particlePtr->clearChannels();
527  particleDataPtr->mayDecay(idRes,false);
528  particleDataPtr->isResonance(idRes,false);
529  }
530 
531  // Loop over SLHA channels, import into Pythia, treating channels
532  // with negative branching fractions as having the equivalent positive
533  // branching fraction, but being switched off for this run
534  for (int iChannel=0 ; iChannel<slhaTable->size(); iChannel++) {
535  LHdecayChannel slhaChannel = slhaTable->getChannel(iChannel);
536  double brat = slhaChannel.getBrat();
537  vector<int> idDa = slhaChannel.getIdDa();
538  if (idDa.size() >= 9) {
539  infoPtr->errorMsg(errPref + "max number of DECAY products is 8 ",
540  "for id = "+idCode.str(), true);
541  } else if (idDa.size() <= 1) {
542  infoPtr->errorMsg(errPref + "min number of DECAY products is 2 ",
543  "for id = "+idCode.str(), true);
544  }
545  else {
546  int onMode = 1;
547  if (brat < 0.0) onMode = 0;
548  int meModeNow = meMode;
549 
550  // Check phase space, including margin ~ sqrt(sum(widths^2))
551  double massSum(0.);
552  double widSqSum = pow2(widRes);
553  int nDa = idDa.size();
554  for (int jDa=0; jDa<nDa; ++jDa) {
555  massSum += particleDataPtr->m0( idDa[jDa] );
556  widSqSum += pow2(particleDataPtr->mWidth( idDa[jDa] ));
557  }
558  double deltaM = particleDataPtr->m0(idRes) - massSum;
559  // Negative mass difference: intrinsically off shell
560  if (onMode == 1 && brat > 0.0 && deltaM < 0.) {
561  // String containing decay name
562  ostringstream errCode;
563  errCode << idRes <<" ->";
564  for (int jDa=0; jDa<nDa; ++jDa) errCode<<" "<<idDa[jDa];
565  // Could mass fluctuations at all give the needed deltaM ?
566  if (abs(deltaM) > 100. * sqrt(widSqSum)) {
567  infoPtr->errorMsg(warnPref + "switched off DECAY mode",
568  ": " + errCode.str()+" (too far off shell)",true);
569  onMode = 0;
570  }
571  // If ~ OK within fluctuations
572  else {
573  // Ignore user-selected meMode
574  if (meModeNow != 100) {
575  infoPtr->errorMsg(warnPref + "adding off shell DECAY mode",
576  ": "+errCode.str()+" (forced meMode = 100)",true);
577  meModeNow = 100;
578  } else {
579  infoPtr->errorMsg(warnPref + "adding off shell DECAY mode",
580  errCode.str(), true);
581  }
582  }
583  }
584 
585  // Add channel
586  int id0 = idDa[0];
587  int id1 = idDa[1];
588  int id2 = (idDa.size() >= 3) ? idDa[2] : 0;
589  int id3 = (idDa.size() >= 4) ? idDa[3] : 0;
590  int id4 = (idDa.size() >= 5) ? idDa[4] : 0;
591  int id5 = (idDa.size() >= 6) ? idDa[5] : 0;
592  int id6 = (idDa.size() >= 7) ? idDa[6] : 0;
593  int id7 = (idDa.size() >= 8) ? idDa[7] : 0;
594  particlePtr->addChannel(onMode,abs(brat),meModeNow,
595  id0,id1,id2,id3,id4,id5,id6,id7);
596 
597  }
598  }
599 
600  // Add to list of particles that have been modified
601  idModified[idRes]=true;
602 
603  }
604 
605  // Give summary of imported/ignored DECAY tables, and state reason
606  if (importDecay.size() >= 1) {
607  string idImport;
608  for (unsigned int i=0; i<importDecay.size(); ++i) {
609  ostringstream idCode;
610  idCode << importDecay[i];
611  if (i != 0) idImport +=",";
612  idImport += idCode.str();
613  }
614  infoPtr->errorMsg(infoPref + "importing DECAY tables","for id = {"
615  + idImport + "}", true);
616  }
617  if (ignoreDecayM0.size() >= 1) {
618  string idIgnore;
619  for (unsigned int i=0; i<ignoreDecayM0.size(); ++i) {
620  ostringstream idCode;
621  idCode << ignoreDecayM0[i];
622  if (i != 0) idIgnore +=",";
623  idIgnore += idCode.str();
624  }
625  infoPtr->errorMsg(warnPref + "ignoring DECAY tables", "for id = {"
626  + idIgnore + "}" + " (m0 < SLHA:minMassSM)", true);
627  }
628  if (ignoreDecayBR.size() >= 1) {
629  string idIgnore;
630  for (unsigned int i=0; i<ignoreDecayBR.size(); ++i) {
631  ostringstream idCode;
632  idCode << ignoreDecayBR[i];
633  if (i != 0) idIgnore +=",";
634  idIgnore += idCode.str();
635  }
636  infoPtr->errorMsg(warnPref + "ignoring empty DECAY tables", "for id = {"
637  + idIgnore + "}" + " (total width provided but no Branching Ratios)",
638  true);
639  }
640 
641  // Sanity check of all decay tables with modified MASS or DECAY info
642  map<int,bool>::iterator it;
643  for (it=idModified.begin(); it!=idModified.end(); ++it) {
644  int id = it->first;
645  if (idModified[id] == false) continue;
646  ostringstream idCode;
647  idCode << id;
648  ParticleDataEntry* particlePtr
649  = particleDataPtr->particleDataEntryPtr(id);
650  double m0 = particlePtr->m0();
651  double wid = particlePtr->mWidth();
652  // Always set massless particles stable
653  if (m0 <= 0.0 && (wid > 0.0 || particlePtr->mayDecay())) {
654  infoPtr->errorMsg(warnPref + "massless particle forced stable",
655  " id = " + idCode.str(), true);
656  particlePtr->clearChannels();
657  particlePtr->setMWidth(0.0);
658  particlePtr->setMayDecay(false);
659  particleDataPtr->isResonance(id,false);
660  continue;
661  }
662  // Declare zero-width particles to be stable (for now)
663  if (wid == 0.0 && particlePtr->mayDecay()) {
664  particlePtr->setMayDecay(false);
665  continue;
666  }
667  // Check at least one on-shell channel is available
668  double mSumMin = 10. * m0;
669  int nChannels = particlePtr->sizeChannels();
670  if (nChannels >= 1) {
671  for (int iChannel=0; iChannel<nChannels; ++iChannel) {
672  DecayChannel channel = particlePtr->channel(iChannel);
673  if (channel.onMode() <= 0) continue;
674  int nProd = channel.multiplicity();
675  double mSum = 0.;
676  for (int iDa = 0; iDa < nProd; ++iDa) {
677  int idDa = channel.product(iDa);
678  mSum += particleDataPtr->m0(idDa);
679  }
680  mSumMin = min(mSumMin, mSum);
681  }
682  // Require at least one on-shell channel
683  if (mSumMin > m0 && particlePtr->id() != 25) {
684  infoPtr->errorMsg(warnPref + "particle forced stable"," id = "
685  + idCode.str() + " (no on-shell decay channels)", true);
686  particlePtr->setMWidth(0.0);
687  particlePtr->setMayDecay(false);
688  continue;
689  }
690  else if (mSumMin > m0 && particlePtr->id() == 25) {
691  infoPtr->errorMsg(warnPref
692  + "allowing particle with no on-shell decays ",
693  " id = " + idCode.str() , true);
694  }
695  else {
696  // mMin: lower cutoff on Breit-Wigner; see above.
697  // Increase minimum if needed to ensure at least one channel on shell
698  double mMin = max(mSumMin, particlePtr->mMin());
699  particlePtr->setMMin(mMin);
700  }
701  }
702  }
703 
704  return true;
705 
706 }
707 
708 //--------------------------------------------------------------------------
709 
710 // Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
711 // E.g., to make sure that there are no important unfilled entries
712 
713 void SLHAinterface::pythia2slha() {
714 
715  // Initialize block SMINPUTS.
716  string blockName = "sminputs";
717  double mZ = particleDataPtr->m0(23);
718  slha.set(blockName,1,1.0/coupSMPtr->alphaEM(pow2(mZ)));
719  slha.set(blockName,2,coupSMPtr->GF());
720  slha.set(blockName,3,coupSMPtr->alphaS(pow2(mZ)));
721  slha.set(blockName,4,mZ);
722  // b mass (should be running mass, here pole for time being)
723  slha.set(blockName,5,particleDataPtr->m0(5));
724  slha.set(blockName,6,particleDataPtr->m0(6));
725  slha.set(blockName,7,particleDataPtr->m0(15));
726  slha.set(blockName,8,particleDataPtr->m0(16));
727  slha.set(blockName,11,particleDataPtr->m0(11));
728  slha.set(blockName,12,particleDataPtr->m0(12));
729  slha.set(blockName,13,particleDataPtr->m0(13));
730  slha.set(blockName,14,particleDataPtr->m0(14));
731  // Force 3 lightest quarks massless
732  slha.set(blockName,21,double(0.0));
733  slha.set(blockName,22,double(0.0));
734  slha.set(blockName,23,double(0.0));
735  // c mass (should be running mass, here pole for time being)
736  slha.set(blockName,24,particleDataPtr->m0(4));
737 
738  // Initialize block MASS.
739  blockName = "mass";
740  int id = 1;
741  int count = 0;
742  while (particleDataPtr->nextId(id) > id) {
743  slha.set(blockName,id,particleDataPtr->m0(id));
744  id = particleDataPtr->nextId(id);
745  ++count;
746  if (count > 10000) {
747  infoPtr->errorMsg("Error in SLHAinterface::pythia2slha(): "
748  "encountered infinite loop when saving mass block");
749  break;
750  }
751  }
752 
753 }
754 
755 //==========================================================================
756 
757 } // end namespace Pythia8
Generated by   doxygen 1.8.5