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