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