StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SusyLesHouches.cc
1 // SusyLesHouches.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/SusyLesHouches.h"
8 #include "Pythia8/Streams.h"
9 
10 namespace Pythia8 {
11 
12 //==========================================================================
13 
14 // The SusyLesHouches class.
15 
16 //--------------------------------------------------------------------------
17 
18 // Main routine to read in SLHA and LHEF+SLHA files
19 
20 int SusyLesHouches::readFile(string slhaFileIn, int verboseIn,
21  bool useDecayIn) {
22 
23  slhaFile = slhaFileIn;
24  // Check that input file is OK.
25  const char* cstring = slhaFile.c_str();
26  igzstream file(cstring);
27 
28  // Exit if input file not found. Else print file name.
29  if (!file.good()) {
30  message(2,"readFile",slhaFile+" not found",0);
31  slhaRead = false;
32  return -1;
33  }
34  if (verboseSav >= 3) {
35  message(0,"readFile","parsing "+slhaFile,0);
36  filePrinted = true;
37  }
38 
39  return readFile( file, verboseIn, useDecayIn );
40 }
41 
42 int SusyLesHouches::readFile(istream& is, int verboseIn,
43  bool useDecayIn) {
44 
45  int iFailFile=0;
46  // Copy inputs to local
47  verboseSav = verboseIn;
48  useDecay = useDecayIn;
49 
50  // Array of particles read in.
51  vector<int> idRead;
52 
53  // Array of block names read in.
54  vector<string> processedBlocks;
55 
56  //Initial values for read-in variables.
57  slhaRead=true;
58  lhefRead=false;
59  lhefSlha=false;
60  bool foundSlhaTag = false;
61  bool xmlComment = false;
62  bool decayPrinted = false;
63  string line="";
64  string blockIn="";
65  string decay="";
66  string comment="";
67  string blockName="";
68  string nameNow="";
69  int idNow=0;
70  double width=0.0;
71 
72  //Initialize line counter
73  int iLine=0;
74 
75  // Read in one line at a time.
76  while ( getline(is, line) ) {
77  iLine++;
78 
79  //Rewrite string in lowercase, removing initial and tralining blanks
80  //as well as garbage characters
81  toLowerRep(line);
82 
83  //Detect whether read-in is from a Les Houches Event File (LHEF).
84  if (line.find("<leshouches") != string::npos
85  || line.find("<slha") != string::npos) {
86  lhefRead=true;
87  }
88 
89  // If LHEF
90  if (lhefRead) {
91  //Ignore XML comments (only works for whole lines so far)
92  if (line.find("-->") != string::npos) {
93  xmlComment = false;
94  }
95  else if (xmlComment) continue;
96  else if (line.find("<!--") != string::npos) {
97  xmlComment = true;
98  }
99  //Detect when <slha> tag reached.
100  if (line.find("<slha") != string::npos) {
101  lhefSlha = true;
102  foundSlhaTag = true;
103  //Print header if not already done
104  if (! headerPrinted) listHeader();
105  }
106  //Stop looking when </header> or <init> tag reached
107  if (line.find("</header>") != string::npos ||
108  line.find("<init") != string::npos) {
109  if (!foundSlhaTag) return 101;
110  break;
111  }
112  //If <slha> tag not yet reached, skip
113  if (!lhefSlha) continue;
114  }
115 
116  //Ignore comment lines with # as first character
117  if (line.find("#") == 0) continue;
118 
119  //Ignore empty lines
120  if (line.size() == 0) continue;
121  if (line.size() == 1 && line.substr(0,1) == " ") continue;
122 
123  //Move comment to separate string
124  if (line.find("#") != string::npos) {
125  if (line.find("#") + 1 < line.length() ) {
126  int commentLength = line.length()-(line.find("#")+1);
127  comment = line.substr(line.find("#")+1,commentLength);
128  } else
129  comment = "";
130  line.erase(line.find("#"),line.length()-line.find("#")-1);
131  }
132 
133  // Remove blanks before and after an = sign. Also remove multiple blanks
134  while (line.find(" =") != string::npos) line.erase( line.find(" ="), 1);
135  while (line.find("= ") != string::npos) line.erase( line.find("= ")+1, 1);
136  while (line.find(" ") != string::npos) line.erase( line.find(" ")+1, 1);
137 
138  //New block.
139  if (line.find("block") <= 1) {
140 
141  //Print header if not already done
142  if (! headerPrinted) listHeader();
143 
144  blockIn=line ;
145  decay="";
146  int nameBegin=6 ;
147  int nameEnd=blockIn.find(" ",7);
148  blockName=blockIn.substr(nameBegin,nameEnd-nameBegin);
149 
150  // QNUMBERS blocks (cf. arXiv:0712.3311 [hep-ph])
151  if (blockIn.find("qnumbers") != string::npos) {
152  // Extract ID code for new particle
153  int pdgBegin=blockIn.find(" ",7)+1;
154  int pdgEnd=blockIn.find(" ",pdgBegin);
155  string pdgString = blockIn.substr(pdgBegin,pdgEnd-pdgBegin);
156  istringstream linestream(pdgString);
157  // Create and add new block with this code as zero'th entry
158  LHblock<double> newQnumbers;
159  newQnumbers.set(0,linestream);
160  qnumbers.push_back(newQnumbers);
161  // Default name: PDG code
162  string defName, defAntiName, newName, newAntiName;
163  ostringstream idStream;
164  idStream << newQnumbers(0);
165  defName = idStream.str();
166  defAntiName = "-"+defName;
167  newName = defName;
168  newAntiName = defAntiName;
169  // Attempt to extract names from comment string
170  if (comment.length() >= 1) {
171  int firstCommentBeg(0), firstCommentEnd(0);
172  if ( comment.find(" ") == 0) firstCommentBeg = 1;
173  if ( comment.find(" ",firstCommentBeg+1) == string::npos)
174  firstCommentEnd = comment.length();
175  else
176  firstCommentEnd = comment.find(" ",firstCommentBeg+1);
177  if (firstCommentEnd > firstCommentBeg)
178  newName = comment.substr(firstCommentBeg,
179  firstCommentEnd-firstCommentBeg);
180  // Now see if there is a separate name for antiparticle
181  int secondCommentBeg(firstCommentEnd+1), secondCommentEnd(0);
182  if (secondCommentBeg < int(comment.length())) {
183  if ( comment.find(" ",secondCommentBeg+1) == string::npos)
184  secondCommentEnd = comment.length();
185  else
186  secondCommentEnd = comment.find(" ",secondCommentBeg+1);
187  if (secondCommentEnd > secondCommentBeg)
188  newAntiName = comment.substr(secondCommentBeg,
189  secondCommentEnd-secondCommentBeg);
190  }
191  }
192  // If name given without specific antiname, set antiname to ""
193  if (newName != defName && newAntiName == defAntiName) newAntiName = "";
194  qnumbersName.push_back(newName);
195  qnumbersAntiName.push_back(newAntiName);
196  if (pdgString != newName) {
197  message(0,"readFile","storing QNUMBERS for id = "+pdgString+" "
198  +newName+" "+newAntiName,iLine);
199  } else {
200  message(0,"readFile","storing QNUMBERS for id = "+pdgString,iLine);
201  }
202  }
203 
204  // Non-qnumbers blocks
205  // Skip if several copies of same block
206  // (facility to use interpolation of different q= not implemented)
207  // only first copy of a given block type is kept
208  else {
209  bool exists = false;
210  for (int i=0; i<int(processedBlocks.size()); ++i)
211  if (blockName == processedBlocks[i]) exists = true;
212  if (exists) {
213  message(0,"readFile","skipping copy of block "+blockName,iLine);
214  blockIn = "";
215  continue;
216  }
217  processedBlocks.push_back(blockName);
218 
219  // Copy input file as generic blocks (containing strings)
220  // (more will be done with SLHA1 & 2 specific blocks below, this is
221  // just to make sure we have a complete copy of the input file,
222  // including also any unknown/user/generic blocks)
223  LHgenericBlock gBlock;
224  genericBlocks[blockName]=gBlock;
225  }
226 
227  //Find Q=... for DRbar running blocks
228  if (blockIn.find("q=") != string::npos) {
229  int qbegin=blockIn.find("q=")+2;
230  istringstream qstream(blockIn.substr(qbegin,blockIn.length()));
231  double q=0.0;
232  qstream >> q;
233  if (qstream) {
234  // SLHA1 running blocks
235  if (blockName=="hmix") hmix.setq(q);
236  if (blockName=="yu") yu.setq(q);
237  if (blockName=="yd") yd.setq(q);
238  if (blockName=="ye") ye.setq(q);
239  if (blockName=="au") au.setq(q);
240  if (blockName=="ad") ad.setq(q);
241  if (blockName=="ae") ae.setq(q);
242  if (blockName=="msoft") msoft.setq(q);
243  if (blockName=="gauge") gauge.setq(q);
244  // SLHA2 running blocks
245  if (blockName=="vckm") vckm.setq(q);
246  if (blockName=="upmns") upmns.setq(q);
247  if (blockName=="msq2") msq2.setq(q);
248  if (blockName=="msu2") msu2.setq(q);
249  if (blockName=="msd2") msd2.setq(q);
250  if (blockName=="msl2") msl2.setq(q);
251  if (blockName=="mse2") mse2.setq(q);
252  if (blockName=="tu") tu.setq(q);
253  if (blockName=="td") td.setq(q);
254  if (blockName=="te") te.setq(q);
255  if (blockName=="rvlamlle") rvlamlle.setq(q);
256  if (blockName=="rvlamlqd") rvlamlqd.setq(q);
257  if (blockName=="rvlamudd") rvlamudd.setq(q);
258  if (blockName=="rvtlle") rvtlle.setq(q);
259  if (blockName=="rvtlqd") rvtlqd.setq(q);
260  if (blockName=="rvtudd") rvtudd.setq(q);
261  if (blockName=="rvkappa") rvkappa.setq(q);
262  if (blockName=="rvd") rvd.setq(q);
263  if (blockName=="rvm2lh1") rvm2lh1.setq(q);
264  if (blockName=="rvsnvev") rvsnvev.setq(q);
265  if (blockName=="imau") imau.setq(q);
266  if (blockName=="imad") imad.setq(q);
267  if (blockName=="imae") imae.setq(q);
268  if (blockName=="imhmix") imhmix.setq(q);
269  if (blockName=="immsoft") immsoft.setq(q);
270  if (blockName=="imtu") imtu.setq(q);
271  if (blockName=="imtd") imtd.setq(q);
272  if (blockName=="imte") imte.setq(q);
273  if (blockName=="imvckm") imvckm.setq(q);
274  if (blockName=="imupmns") imupmns.setq(q);
275  if (blockName=="immsq2") immsq2.setq(q);
276  if (blockName=="immsu2") immsu2.setq(q);
277  if (blockName=="immsd2") immsd2.setq(q);
278  if (blockName=="immsl2") immsl2.setq(q);
279  if (blockName=="immse2") immse2.setq(q);
280  if (blockName=="nmssmrun") nmssmrun.setq(q);
281  };
282  };
283 
284  //Skip to next line.
285  continue ;
286 
287  }
288 
289  //New decay table
290  else if (line.find("decay") <= 1) {
291 
292  // Print header if not already done
293  if (! headerPrinted) listHeader();
294 
295  // If previous had zero length, print now
296  if (decay != "" && ! decayPrinted) {
297  if (verboseSav >= 2) message(0,"readFile","reading WIDTH for "+nameNow
298  +" (but no decay channels found)",0);
299  }
300 
301  //Set decay block name
302  decay=line;
303  blockIn="";
304  int nameBegin=6 ;
305  int nameEnd=decay.find(" ",7);
306  nameNow=decay.substr(nameBegin,nameEnd-nameBegin);
307 
308  //Extract PDG code and width
309  istringstream dstream(nameNow);
310  dstream >> idNow;
311 
312  //Ignore decay if decay table read-in switched off
313  if( !useDecay ) {
314  decay = "";
315  message(0,"readFile","ignoring DECAY table for "+nameNow
316  +" (DECAY read-in switched off)",iLine);
317  continue;
318  }
319 
320  if (dstream) {
321  string widthName=decay.substr(nameEnd+1,decay.length());
322  istringstream wstream(widthName);
323  wstream >> width;
324  if (wstream) {
325  // Set
326  decays.push_back(LHdecayTable(idNow,width));
327  decayIndices[idNow]=decays.size()-1;
328  //Set PDG code and width
329  if (width <= 0.0) {
330  string endComment="";
331  if (width < -1e-6) {
332  endComment="(forced width < 0 to zero)";
333  }
334  if (verboseSav >= 2)
335  message(0,"readFile","reading stable particle "+nameNow
336  +" "+endComment,0);
337  width=0.0;
338  decayPrinted = true;
339  decays[decayIndices[idNow]].setWidth(width);
340  } else {
341  decayPrinted = false;
342  }
343  } else {
344  if (verboseSav >= 2)
345  message(0,"readFile","ignoring DECAY table for "+nameNow
346  +" (read failed)",iLine);
347  decayPrinted = true;
348  width=0.0;
349  decay="";
350  continue;
351  }
352  }
353  else {
354  message(0,"readFile",
355  "PDG Code unreadable. Ignoring this DECAY block",iLine);
356  decayPrinted = true;
357  decay="";
358  continue;
359  }
360 
361  //Skip to next line
362  continue ;
363  }
364 
365  //Switch off SLHA read-in via LHEF if outside <slha> tag.
366  else if (line.find("</slha>") != string::npos) {
367  lhefSlha=false;
368  blockIn="";
369  decay="";
370  continue;
371  }
372 
373  //Skip not currently reading block data lines.
374  if (blockIn != "") {
375 
376  // Replace an equal sign by a blank to make parsing simpler.
377  while (line.find("=") != string::npos) {
378  int firstEqual = line.find_first_of("=");
379  line.replace(firstEqual, 1, " ");
380  };
381 
382  //Parse data lines within given block
383  //Constructed explicitly so that each block can have its own types and
384  //own rules defined. For extra user blocks, just add more recognized
385  //blockNames at the end and implement user defined rules accordingly.
386  //string comment = line.substr(line.find("#"),line.length());
387  int ifail=-2;
388  istringstream linestream(line);
389 
390  // Read line in QNUMBERS block, add to end of qnumbers vector
391  if (blockName == "qnumbers") {
392  int iEnd = qnumbers.size()-1;
393  if (iEnd >= 0) ifail = qnumbers[iEnd].set(linestream);
394  else ifail = -1;
395  }
396 
397  // MODEL
398  else if (blockName == "modsel") {
399  int i;
400  linestream >> i;
401  if (linestream) {
402  if (i == 12) {ifail=modsel12.set(0,linestream);}
403  else if (i == 21) {ifail=modsel21.set(0,linestream);}
404  else {ifail=modsel.set(i,linestream);};}
405  else {
406  ifail = -1;}
407  };
408  if (blockName == "minpar") ifail=minpar.set(linestream);
409  if (blockName == "sminputs") ifail=sminputs.set(linestream);
410  if (blockName == "extpar") ifail=extpar.set(linestream);
411  if (blockName == "qextpar") ifail=qextpar.set(linestream);
412  //FLV
413  if (blockName == "vckmin") ifail=vckmin.set(linestream);
414  if (blockName == "upmnsin") ifail=upmnsin.set(linestream);
415  if (blockName == "msq2in") ifail=msq2in.set(linestream);
416  if (blockName == "msu2in") ifail=msu2in.set(linestream);
417  if (blockName == "msd2in") ifail=msd2in.set(linestream);
418  if (blockName == "msl2in") ifail=msl2in.set(linestream);
419  if (blockName == "mse2in") ifail=mse2in.set(linestream);
420  if (blockName == "tuin") ifail=tuin.set(linestream);
421  if (blockName == "tdin") ifail=tdin.set(linestream);
422  if (blockName == "tein") ifail=tein.set(linestream);
423  //RPV
424  if (blockName == "rvlamllein") ifail=rvlamllein.set(linestream);
425  if (blockName == "rvlamlqdin") ifail=rvlamlqdin.set(linestream);
426  if (blockName == "rvlamuddin") ifail=rvlamuddin.set(linestream);
427  if (blockName == "rvtllein") ifail=rvtllein.set(linestream);
428  if (blockName == "rvtlqdin") ifail=rvtlqdin.set(linestream);
429  if (blockName == "rvtuddin") ifail=rvtuddin.set(linestream);
430  if (blockName == "rvkappain") ifail=rvkappain.set(linestream);
431  if (blockName == "rvdin") ifail=rvdin.set(linestream);
432  if (blockName == "rvm2lh1in") ifail=rvm2lh1in.set(linestream);
433  if (blockName == "rvsnvevin") ifail=rvsnvevin.set(linestream);
434  //CPV
435  if (blockName == "imminpar") ifail=imminpar.set(linestream);
436  if (blockName == "imextpar") ifail=imextpar.set(linestream);
437  //CPV +FLV
438  if (blockName == "immsq2in") ifail=immsq2in.set(linestream);
439  if (blockName == "immsu2in") ifail=immsu2in.set(linestream);
440  if (blockName == "immsd2in") ifail=immsd2in.set(linestream);
441  if (blockName == "immsl2in") ifail=immsl2in.set(linestream);
442  if (blockName == "immse2in") ifail=immse2in.set(linestream);
443  if (blockName == "imtuin") ifail=imtuin.set(linestream);
444  if (blockName == "imtdin") ifail=imtdin.set(linestream);
445  if (blockName == "imtein") ifail=imtein.set(linestream);
446  //Info:
447  if (blockName == "spinfo" || blockName=="dcinfo") {
448  int i;
449  string entry;
450  linestream >> i >> entry;
451  string blockStr="RGE";
452  if (blockName=="dcinfo") blockStr="DCY";
453 
454  if (linestream) {
455  if ( i == 3 ) {
456  string warning=line.substr(line.find("3")+1,line.length());
457  message(1,"readFile","(from "+blockStr+" program): "+warning,0);
458  if (blockName == "spinfo") spinfo3.set(warning);
459  else dcinfo3.set(warning);
460  } else if ( i == 4 ) {
461  string error=line.substr(line.find("4")+1,line.length());
462  message(2,"readFile","(from "+blockStr+" program): "+error,0);
463  if (blockName == "spinfo") spinfo4.set(error);
464  else dcinfo4.set(error);
465  } else {
466  //Rewrite string in uppercase
467  for (unsigned int j=0;j<entry.length();j++)
468  entry[j]=toupper(entry[j]);
469  ifail=(blockName=="spinfo") ? spinfo.set(i,entry)
470  : dcinfo.set(i,entry);
471  };
472  } else {
473  ifail=-1;
474  };
475  };
476  //SPECTRUM
477  //Pole masses
478  if (blockName == "mass") ifail=mass.set(linestream);
479 
480  //Mixing
481  if (blockName == "alpha") ifail=alpha.set(linestream,false);
482  if (blockName == "stopmix") ifail=stopmix.set(linestream);
483  if (blockName == "sbotmix") ifail=sbotmix.set(linestream);
484  if (blockName == "staumix") ifail=staumix.set(linestream);
485  if (blockName == "nmix") ifail=nmix.set(linestream);
486  if (blockName == "umix") ifail=umix.set(linestream);
487  if (blockName == "vmix") ifail=vmix.set(linestream);
488  //FLV
489  if (blockName == "usqmix") ifail=usqmix.set(linestream);
490  if (blockName == "dsqmix") ifail=dsqmix.set(linestream);
491  if (blockName == "selmix") ifail=selmix.set(linestream);
492  if (blockName == "snumix") ifail=snumix.set(linestream);
493  if (blockName == "snsmix") ifail=snsmix.set(linestream);
494  if (blockName == "snamix") ifail=snamix.set(linestream);
495  //RPV
496  if (blockName == "rvnmix") ifail=rvnmix.set(linestream);
497  if (blockName == "rvumix") ifail=rvumix.set(linestream);
498  if (blockName == "rvvmix") ifail=rvvmix.set(linestream);
499  if (blockName == "rvhmix") ifail=rvhmix.set(linestream);
500  if (blockName == "rvamix") ifail=rvamix.set(linestream);
501  if (blockName == "rvlmix") ifail=rvlmix.set(linestream);
502  //CPV
503  if (blockName == "cvhmix") ifail=cvhmix.set(linestream);
504  if (blockName == "imcvhmix") ifail=imcvhmix.set(linestream);
505  //CPV + FLV
506  if (blockName == "imusqmix") ifail=imusqmix.set(linestream);
507  if (blockName == "imdsqmix") ifail=imdsqmix.set(linestream);
508  if (blockName == "imselmix") ifail=imselmix.set(linestream);
509  if (blockName == "imsnumix") ifail=imsnumix.set(linestream);
510  if (blockName == "imnmix") ifail=imnmix.set(linestream);
511  if (blockName == "imumix") ifail=imumix.set(linestream);
512  if (blockName == "imvmix") ifail=imvmix.set(linestream);
513  //NMSSM
514  if (blockName == "nmhmix") ifail=nmhmix.set(linestream);
515  if (blockName == "nmamix") ifail=nmamix.set(linestream);
516  if (blockName == "nmnmix") ifail=nmnmix.set(linestream);
517 
518  //DRbar Lagrangian parameters
519  if (blockName == "gauge") ifail=gauge.set(linestream);
520  if (blockName == "yu") ifail=yu.set(linestream);
521  if (blockName == "yd") ifail=yd.set(linestream);
522  if (blockName == "ye") ifail=ye.set(linestream);
523  if (blockName == "au") ifail=au.set(linestream);
524  if (blockName == "ad") ifail=ad.set(linestream);
525  if (blockName == "ae") ifail=ae.set(linestream);
526  if (blockName == "hmix") ifail=hmix.set(linestream);
527  if (blockName == "msoft") ifail=msoft.set(linestream);
528  //FLV
529  if (blockName == "vckm") ifail=vckm.set(linestream);
530  if (blockName == "upmns") ifail=upmns.set(linestream);
531  if (blockName == "msq2") ifail=msq2.set(linestream);
532  if (blockName == "msu2") ifail=msu2.set(linestream);
533  if (blockName == "msd2") ifail=msd2.set(linestream);
534  if (blockName == "msl2") ifail=msl2.set(linestream);
535  if (blockName == "mse2") ifail=mse2.set(linestream);
536  if (blockName == "tu") ifail=tu.set(linestream);
537  if (blockName == "td") ifail=td.set(linestream);
538  if (blockName == "te") ifail=te.set(linestream);
539  //RPV
540  if (blockName == "rvlamlle") ifail=rvlamlle.set(linestream);
541  if (blockName == "rvlamlqd") ifail=rvlamlqd.set(linestream);
542  if (blockName == "rvlamudd") ifail=rvlamudd.set(linestream);
543  if (blockName == "rvtlle") ifail=rvtlle.set(linestream);
544  if (blockName == "rvtlqd") ifail=rvtlqd.set(linestream);
545  if (blockName == "rvtudd") ifail=rvtudd.set(linestream);
546  if (blockName == "rvkappa") ifail=rvkappa.set(linestream);
547  if (blockName == "rvd") ifail=rvd.set(linestream);
548  if (blockName == "rvm2lh1") ifail=rvm2lh1.set(linestream);
549  if (blockName == "rvsnvev") ifail=rvsnvev.set(linestream);
550  //CPV
551  if (blockName == "imau") ifail=imau.set(linestream);
552  if (blockName == "imad") ifail=imad.set(linestream);
553  if (blockName == "imae") ifail=imae.set(linestream);
554  if (blockName == "imhmix") ifail=imhmix.set(linestream);
555  if (blockName == "immsoft") ifail=immsoft.set(linestream);
556  //CPV+FLV
557  if (blockName == "imvckm") ifail=imvckm.set(linestream);
558  if (blockName == "imupmns") ifail=imupmns.set(linestream);
559  if (blockName == "immsq2") ifail=immsq2.set(linestream);
560  if (blockName == "immsu2") ifail=immsu2.set(linestream);
561  if (blockName == "immsd2") ifail=immsd2.set(linestream);
562  if (blockName == "immsl2") ifail=immsl2.set(linestream);
563  if (blockName == "immse2") ifail=immse2.set(linestream);
564  if (blockName == "imtu") ifail=imtu.set(linestream);
565  if (blockName == "imtd") ifail=imtd.set(linestream);
566  if (blockName == "imte") ifail=imte.set(linestream);
567  //NMSSM
568  if (blockName == "nmssmrun") ifail=nmssmrun.set(linestream);
569 
570  //Diagnostics
571  if (ifail != 0) {
572  if (ifail == -2 && !genericBlocks[blockName].exists() ) {
573  message(0,"readFile","storing non-SLHA(2) block: "+blockName,iLine);
574  };
575  if (ifail == -1) {
576  message(1,"readFile","read error or empty line",iLine);
577  };
578  if (ifail == 1) {
579  message(0,"readFile",blockName+" existing entry overwritten",iLine);
580  };
581  }
582 
583  // Add line to generic block (carbon copy of input structure)
584  // NB: do not save empty lines, defined as having length <= 1
585  if (line.size() >= 2) {
586  genericBlocks[blockName].set(line);
587  }
588 
589  }
590 
591  // Decay table read-in
592  else if (decay != "") {
593  if (! decayPrinted) {
594  if (verboseSav >= 2)
595  message(0,"readFile","reading DECAY table for "+nameNow,0);
596  decayPrinted = true;
597  }
598  double brat;
599  bool ok=true;
600  int nDa = 0;
601  vector<int> idDa;
602  istringstream linestream(line);
603  linestream >> brat;
604  if (! linestream) ok = false;
605  if (ok) linestream >> nDa;
606  if (! linestream) ok = false;
607  else {
608  for (int i=0; i<nDa; i++) {
609  int idThis;
610  linestream >> idThis;
611  if (! linestream) {
612  ok = false;
613  break;
614  }
615  idDa.push_back(idThis);
616  }
617  }
618 
619  // Stop reading decay channels if not consistent.
620  if (!ok || nDa < 2) {
621  message(1,"readFile","read error or empty line",iLine);
622 
623  // Append decay channel.
624  } else {
625  decays[decayIndices[idNow]].addChannel(brat,nDa,idDa);
626  }
627  }
628  };
629 
630  //Print footer
631  listFooter();
632 
633  //Return 0 if read-in successful
634  if ( lhefRead && !foundSlhaTag) {
635  return 102;
636  }
637  else return iFailFile;
638 
639 }
640 
641 //--------------------------------------------------------------------------
642 
643 // Print a header with information on version, last date of change, etc.
644 
645 void SusyLesHouches::listHeader() {
646  if (verboseSav == 0) return;
647  cout << setprecision(3);
648  if (! headerPrinted) {
649  cout << " *----------------------- SusyLesHouches SUSY/BSM"
650  << " Interface ------------------------*\n";
651  message(0,"","Last Change 12 Apr 2017 - P. Skands",0);
652  if (!filePrinted && slhaFile != "" && slhaFile != " ") {
653  message(0,"","Parsing: "+slhaFile,0);
654  filePrinted=true;
655  }
656  headerPrinted=true;
657  }
658 }
659 
660 //--------------------------------------------------------------------------
661 
662 // Print a footer
663 
664 void SusyLesHouches::listFooter() {
665  if (verboseSav == 0) return;
666  if (! footerPrinted) {
667  // cout << " *" << endl;
668  cout << " *-----------------------------------------------------"
669  << "-------------------------------*\n";
670  footerPrinted=true;
671  // headerPrinted=false;
672  }
673 }
674 
675 //--------------------------------------------------------------------------
676 
677 // Print the current spectrum on stdout.
678 // Not yet fully implemented.
679 
680 void SusyLesHouches::listSpectrum(int ifail) {
681 
682  // Exit if output switched off
683  if (verboseSav <= 0) return;
684 
685  // Print header if not already done
686  if (! headerPrinted) listHeader();
687  message(0,"","");
688 
689  // Print Calculator and File name
690  if (slhaRead) {
691  message(0,""," Spectrum Calculator was: "+spinfo(1)+" version: "
692  +spinfo(2));
693  if (lhefRead) message(0,""," Read <slha> spectrum from: "+slhaFile);
694  else message(0,""," Read SLHA spectrum from: "+slhaFile);
695  }
696 
697  // Failed?
698  if (ifail < 0) {
699  message(0,""," Check revealed problems. Only using masses.");
700  }
701 
702  // gluino
703  message(0,"","");
704  cout << " | ~g m" << endl;
705  cout << setprecision(3) << " | 1000021 " << setw(10) <<
706  ( (mass(2000003) > 1e7) ? scientific : fixed) << mass(1000021) << endl;
707 
708  // d squarks
709  message(0,"","");
710  cout << " | ~d m |~dL| |~sL| |~bL|"
711  << " |~dR| |~sR| |~bR|" << endl;
712 
713  cout << setprecision(3) << " | 1000001 " << setw(10)
714  << ( (mass(1000001) > 1e7) ? scientific : fixed) << mass(1000001)
715  << fixed << " ";
716  for (int icur=1;icur<=6;icur++) cout << setw(6)
717  << sqrt(pow2(dsqmix(1,icur))+pow2(imdsqmix(1,icur))) << " ";
718 
719  cout << endl << " | 1000003 " << setw(10)
720  << ( (mass(1000003) > 1e7) ? scientific : fixed) << mass(1000003)
721  << fixed << " ";
722  for (int icur=1;icur<=6;icur++) cout << setw(6)
723  << sqrt(pow2(dsqmix(2,icur))+pow2(imdsqmix(2,icur))) << " ";
724 
725  cout << endl << " | 1000005 " << setw(10)
726  << ( (mass(1000005) > 1e7) ? scientific : fixed) << mass(1000005)
727  << fixed << " ";
728  for (int icur=1;icur<=6;icur++) cout << setw(6)
729  << sqrt(pow2(dsqmix(3,icur))+pow2(imdsqmix(3,icur))) << " ";
730 
731  cout << endl << " | 2000001 " << setw(10)
732  << ( (mass(2000001) > 1e7) ? scientific : fixed) << mass(2000001)
733  << fixed << " ";
734  for (int icur=1;icur<=6;icur++) cout << setw(6)
735  << sqrt(pow2(dsqmix(4,icur))+pow2(imdsqmix(4,icur))) << " ";
736 
737  cout << endl << " | 2000003 " << setw(10)
738  << ( (mass(2000003) > 1e7) ? scientific : fixed) << mass(2000003)
739  << fixed << " ";
740  for (int icur=1;icur<=6;icur++) cout << setw(6)
741  << sqrt(pow2(dsqmix(5,icur))+pow2(imdsqmix(5,icur))) << " ";
742 
743  cout << endl << " | 2000005 " << setw(10)
744  << ( (mass(2000005) > 1e7) ? scientific : fixed) << mass(2000005)
745  << fixed << " ";
746  for (int icur=1;icur<=6;icur++) cout << setw(6)
747  << sqrt(pow2(dsqmix(6,icur))+pow2(imdsqmix(6,icur))) << " ";
748 
749  cout << endl;
750 
751  // u squarks
752  message(0,"","");
753  cout << " | ~u m |~uL| |~cL| |~tL|"
754  << " |~uR| |~cR| |~tR|" << endl;
755 
756  cout << setprecision(3) << " | 1000002 " << setw(10)
757  << ( (mass(1000002) > 1e7) ? scientific : fixed) << mass(1000002)
758  << fixed << " ";
759  for (int icur=1;icur<=6;icur++) cout <<setw(6)
760  << sqrt(pow2(usqmix(1,icur))+pow2(imusqmix(1,icur))) << " ";
761 
762  cout << endl << " | 1000004 " << setw(10)
763  << ( (mass(1000004) > 1e7) ? scientific : fixed) << mass(1000004)
764  << fixed << " ";
765  for (int icur=1;icur<=6;icur++) cout << setw(6)
766  << sqrt(pow2(usqmix(2,icur))+pow2(imusqmix(2,icur))) << " ";
767 
768  cout << endl << " | 1000006 " << setw(10)
769  << ( (mass(1000006) > 1e7) ? scientific : fixed) << mass(1000006)
770  << fixed << " ";
771  for (int icur=1;icur<=6;icur++) cout << setw(6)
772  << sqrt(pow2(usqmix(3,icur))+pow2(imusqmix(3,icur))) << " ";
773 
774  cout << endl << " | 2000002 " << setw(10)
775  << ( (mass(2000002) > 1e7) ? scientific : fixed) << mass(2000002)
776  << fixed << " ";
777  for (int icur=1;icur<=6;icur++) cout << setw(6)
778  << sqrt(pow2(usqmix(4,icur))+pow2(imusqmix(4,icur))) << " ";
779 
780  cout << endl << " | 2000004 " << setw(10)
781  << ( (mass(2000004) > 1e7) ? scientific : fixed) << mass(2000004)
782  << fixed << " " ;
783  for (int icur=1;icur<=6;icur++) cout << setw(6)
784  << sqrt(pow2(usqmix(5,icur))+pow2(imusqmix(5,icur))) << " ";
785 
786  cout << endl << " | 2000006 " << setw(10)
787  << ( (mass(2000006) > 1e7) ? scientific : fixed) << mass(2000006)
788  << fixed << " ";
789  for (int icur=1;icur<=6;icur++) cout << setw(6)
790  << sqrt(pow2(usqmix(6,icur))+pow2(imusqmix(6,icur))) << " ";
791 
792  cout << endl;
793 
794  // Charged scalars (sleptons)
795  message(0,"","");
796 
797  // R-conserving:
798  if (modsel(4) < 1) {
799  cout << " | ~e m |~eL| |~muL| |~tauL|"
800  << " |~eR| |~muR| |~tauR|" << endl;
801 
802  cout << setprecision(3) << " | 1000011 " << setw(10)
803  << ( (mass(1000011) > 1e7) ? scientific : fixed) << mass(1000011)
804  << fixed << " ";
805  for (int icur=1;icur<=6;icur++) cout << setw(6)
806  << sqrt(pow2(selmix(1,icur))+pow2(imselmix(1,icur))) << " ";
807 
808  cout << endl << " | 1000013 " << setw(10)
809  << ( (mass(1000013) > 1e7) ? scientific : fixed) << mass(1000013)
810  << fixed << " ";
811  for (int icur=1;icur<=6;icur++) cout << setw(6)
812  << sqrt(pow2(selmix(2,icur))+pow2(imselmix(2,icur))) << " ";
813 
814  cout << endl << " | 1000015 " << setw(10)
815  << ( (mass(1000015) > 1e7) ? scientific : fixed) << mass(1000015)
816  << fixed << " ";
817  for (int icur=1;icur<=6;icur++) cout << setw(6)
818  << sqrt(pow2(selmix(3,icur))+pow2(imselmix(3,icur))) << " ";
819 
820  cout << endl << " | 2000011 " << setw(10)
821  << ( (mass(2000011) > 1e7) ? scientific : fixed) << mass(2000011)
822  << fixed << " ";
823  for (int icur=1;icur<=6;icur++) cout << setw(6)
824  << sqrt(pow2(selmix(4,icur))+pow2(imselmix(4,icur))) << " ";
825 
826  cout << endl << " | 2000013 " << setw(10)
827  << ( (mass(2000013) > 1e7) ? scientific : fixed) << mass(2000013)
828  << fixed << " " ;
829  for (int icur=1;icur<=6;icur++) cout << setw(6)
830  << sqrt(pow2(selmix(5,icur))+pow2(imselmix(5,icur))) << " ";
831 
832  cout << endl << " | 2000015 " << setw(10)
833  << ( (mass(2000015) > 1e7) ? scientific : fixed) << mass(2000015)
834  << fixed << " ";
835  for (int icur=1;icur<=6;icur++) cout << setw(6)
836  << sqrt(pow2(selmix(6,icur))+pow2(imselmix(6,icur)))<< " ";
837  }
838 
839  // R-violating
840  else {
841  cout << " | H-/~e m H1- H2- ~eL ~muL"
842  << " ~tauL ~eR ~muR ~tauR" << endl;
843 
844  cout << setprecision(3) << " | -37 " << setw(10) <<
845  ( (mass(37) > 1e7) ? scientific : fixed) << mass(37) << fixed << " ";
846  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(1,icur) << " ";
847 
848  cout << endl << " | 1000011 " << setw(10)
849  << ( (mass(1000011) > 1e7) ? scientific : fixed) << mass(1000011)
850  << fixed << " ";
851  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(2,icur) << " ";
852 
853  cout << endl << " | 1000013 " << setw(10)
854  << ( (mass(1000013) > 1e7) ? scientific : fixed) << mass(1000013)
855  << fixed << " ";
856  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(3,icur) << " ";
857 
858  cout << endl << " | 1000015 " << setw(10)
859  << ( (mass(1000015) > 1e7) ? scientific : fixed) << mass(1000015)
860  << fixed << " ";
861  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(4,icur) << " ";
862 
863  cout << endl << " | 2000011 " << setw(10)
864  << ( (mass(2000011) > 1e7) ? scientific : fixed) << mass(2000011)
865  << fixed << " ";
866  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(5,icur) << " ";
867 
868  cout << endl << " | 2000013 " << setw(10)
869  << ( (mass(2000013) > 1e7) ? scientific : fixed) << mass(2000013)
870  << fixed << " " ;
871  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(6,icur) << " ";
872 
873  cout << endl << " | 2000015 " << setw(10)
874  << ( (mass(2000015) > 1e7) ? scientific : fixed) << mass(2000015)
875  << fixed << " ";
876  for (int icur=1;icur<=8;icur++) cout << setw(6) << rvlmix(7,icur) << " ";
877  }
878  cout << endl;
879 
880  // Neutral scalars (sneutrinos)
881  message(0,"","");
882 
883  // R-conserving:
884  if (modsel(4) < 1) {
885  cout << " | ~nu m";
886  if (snumix.exists()) cout << " |~nu_e| |~nu_mu||~nu_tau|";
887  cout << endl;
888 
889  cout << setprecision(3) << " | 1000012 " << setw(10)
890  << ( (mass(1000012) > 1e7) ? scientific : fixed) << mass(1000012)
891  << fixed << " ";
892  if (snumix.exists())
893  for (int icur=1;icur<=3;icur++) cout << setw(7)
894  << sqrt(pow2(snumix(1,icur))+pow2(imsnumix(1,icur))) << " ";
895 
896  cout << endl << " | 1000014 " << setw(10)
897  << ( (mass(1000014) > 1e7) ? scientific : fixed) << mass(1000014)
898  << fixed << " ";
899  if (snumix.exists())
900  for (int icur=1;icur<=3;icur++) cout << setw(7)
901  << sqrt(pow2(snumix(2,icur))+pow2(imsnumix(2,icur))) << " ";
902 
903  cout << endl << " | 1000016 " << setw(10)
904  << ( (mass(1000016) > 1e7) ? scientific : fixed) << mass(1000016)
905  << fixed << " ";
906  if (snumix.exists())
907  for (int icur=1;icur<=3;icur++) cout << setw(7)
908  << sqrt(pow2(snumix(3,icur))+pow2(imsnumix(3,icur))) << " ";
909  }
910 
911  // R-violating
912  else {
913  cout << " | H0/~nu m";
914  if (snumix.exists()) cout << " H0_1 H0_2 ~nu_e ~nu_mu ~nu_tau";
915  cout << endl;
916 
917  cout << setprecision(3) << " | 25 " << setw(10)
918  << ( (mass(25) > 1e7) ? scientific : fixed) << mass(25)
919  << fixed << " ";
920  if (rvhmix.exists()) for (int icur=1;icur<=5;icur++)
921  cout << setw(6) << rvhmix(1,icur) << " ";
922 
923  cout << endl << " | 35 " << setw(10)
924  << ( (mass(35) > 1e7) ? scientific : fixed) << mass(35)
925  << fixed << " ";
926  if (rvhmix.exists()) for (int icur=1;icur<=5;icur++)
927  cout << setw(6) << rvhmix(2,icur) << " ";
928 
929  cout << endl << " | 1000012 " << setw(10)
930  << ( (mass(1000012) > 1e7) ? scientific : fixed) << mass(1000012)
931  << fixed << " ";
932  if (rvhmix.exists()) for (int icur=1;icur<=5;icur++)
933  cout << setw(6) << rvhmix(3,icur) << " ";
934 
935  cout << endl << " | 1000014 " << setw(10)
936  << ( (mass(1000014) > 1e7) ? scientific : fixed) << mass(1000014)
937  << fixed << " ";
938  if (rvhmix.exists()) for (int icur=1;icur<=5;icur++)
939  cout << setw(6) << rvhmix(4,icur) << " ";
940 
941  cout << endl << " | 1000016 " << setw(10)
942  << ( (mass(1000016) > 1e7) ? scientific : fixed) << mass(1000016)
943  << fixed << " ";
944  if (rvhmix.exists()) for (int icur=1;icur<=5;icur++)
945  cout << setw(6) << rvhmix(5,icur) << " ";
946  }
947  cout << endl;
948 
949  // Neutral pseudoscalars (RPV only)
950  if (modsel(4) >= 1 && rvamix.exists()) {
951  message(0,"","");
952  cout << " | A0/~nu m A0_1 A0_2 ~nu_e ~nu_mu ~nu_tau"
953  << endl;
954 
955  cout << setprecision(3) << " | 36 " << setw(10)
956  << ( (mass(36) > 1e7) ? scientific : fixed) << mass(36)
957  << fixed << " ";
958  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvamix(1,icur) << " ";
959 
960  cout << endl << " | 1000017 " << setw(10)
961  << ( (mass(1000017) > 1e7) ? scientific : fixed) << mass(1000017)
962  << fixed << " ";
963  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvamix(2,icur) << " ";
964 
965  cout << endl << " | 1000018 " << setw(10)
966  << ( (mass(1000018) > 1e7) ? scientific : fixed) << mass(1000018)
967  << fixed << " ";
968  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvamix(3,icur) << " ";
969 
970  cout << endl << " | 1000019 " << setw(10)
971  << ( (mass(1000019) > 1e7) ? scientific : fixed) << mass(1000019)
972  << fixed << " ";
973  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvamix(4,icur) << " ";
974  cout << endl;
975 
976  }
977 
978  // Neutral fermions (neutralinos)
979  message(0,"","");
980 
981  // NMSSM
982  if (modsel(3) >= 1) {
983  cout << " | ~chi0 m |~B| |~W_3| |~H_1| |~H_2| |~S|"
984  << endl;
985 
986  cout << setprecision(3) << " | 1000022 " << setw(10)
987  << ( (mass(1000022) > 1e7) ? scientific : fixed) << mass(1000022)
988  << fixed << " ";
989  for (int icur=1;icur<=5;icur++)
990  cout <<setw(6)<< sqrt(pow2(nmnmix(1,icur))+pow2(imnmnmix(1,icur)))<<" ";
991 
992  cout << endl << " | 1000023 " << setw(10)
993  << ( (mass(1000023) > 1e7) ? scientific : fixed) << mass(1000023)
994  << fixed << " ";
995  for (int icur=1;icur<=5;icur++)
996  cout <<setw(6)<< sqrt(pow2(nmnmix(2,icur))+pow2(imnmnmix(2,icur)))<<" ";
997 
998  cout << endl << " | 1000025 " << setw(10)
999  << ( (mass(1000025) > 1e7) ? scientific : fixed) << mass(1000025)
1000  << fixed << " ";
1001  for (int icur=1;icur<=5;icur++) cout << setw(6)
1002  << sqrt(pow2(nmnmix(3,icur))+pow2(imnmnmix(3,icur))) << " ";
1003 
1004  cout << endl << " | 1000035 " << setw(10)
1005  << ( (mass(1000035) > 1e7) ? scientific : fixed) << mass(1000035)
1006  << fixed << " ";
1007  for (int icur=1;icur<=5;icur++) cout << setw(6)
1008  << sqrt(pow2(nmnmix(4,icur))+pow2(imnmnmix(4,icur))) << " ";
1009 
1010  cout << endl << " | 1000045 " << setw(10)
1011  << ( (mass(1000045) > 1e7) ? scientific : fixed) << mass(1000045)
1012  << fixed << " ";
1013  for (int icur=1;icur<=5;icur++) cout << setw(6)
1014  << sqrt(pow2(nmnmix(5,icur))+pow2(imnmnmix(5,icur))) << " ";
1015 
1016  }
1017 
1018  // R-Conserving MSSM
1019  else if (modsel(4) < 1) {
1020  cout << " | ~chi0 m |~B| |~W_3| |~H_1| |~H_2|"
1021  << endl;
1022 
1023  cout << setprecision(3) << " | 1000022 " << setw(10)
1024  << ( (mass(1000022) > 1e7) ? scientific : fixed) << mass(1000022)
1025  << fixed << " ";
1026  for (int icur=1;icur<=4;icur++)
1027  cout << setw(6) << sqrt(pow2(nmix(1,icur))+pow2(imnmix(1,icur))) << " ";
1028 
1029  cout << endl << " | 1000023 " << setw(10)
1030  << ( (mass(1000023) > 1e7) ? scientific : fixed) << mass(1000023)
1031  << fixed << " ";
1032  for (int icur=1;icur<=4;icur++)
1033  cout << setw(6) << sqrt(pow2(nmix(2,icur))+pow2(imnmix(2,icur))) << " ";
1034 
1035  cout << endl << " | 1000025 " << setw(10)
1036  << ( (mass(1000025) > 1e7) ? scientific : fixed) << mass(1000025)
1037  << fixed << " ";
1038  for (int icur=1;icur<=4;icur++)
1039  cout << setw(6) << sqrt(pow2(nmix(3,icur))+pow2(imnmix(3,icur))) << " ";
1040 
1041  cout << endl << " | 1000035 " << setw(10)
1042  << ( (mass(1000035) > 1e7) ? scientific : fixed) << mass(1000035)
1043  << fixed << " ";
1044  for (int icur=1;icur<=4;icur++)
1045  cout << setw(6) << sqrt(pow2(nmix(4,icur))+pow2(imnmix(4,icur))) << " ";
1046 
1047  }
1048 
1049  // R-violating MSSM
1050  else {
1051  cout << " | nu/~chi0 m nu_e nu_mu nu_tau ~B"
1052  << " ~W_3 ~H_1 ~H_2" << endl;
1053 
1054  cout << setprecision(3) << " | 12 " << setw(10)
1055  << ( (mass(12) > 1e7) ? scientific : fixed) << mass(12)
1056  << fixed << " ";
1057  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(1,icur) << " ";
1058 
1059  cout << endl << " | 14 " << setw(10)
1060  << ( (mass(14) > 1e7) ? scientific : fixed) << mass(14)
1061  << fixed << " ";
1062  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(2,icur) << " ";
1063 
1064  cout << endl << " | 16 " << setw(10) <<
1065  ( (mass(16) > 1e7) ? scientific : fixed) << mass(16) << fixed << " ";
1066  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(3,icur) << " ";
1067 
1068  cout << endl << " | 1000022 " << setw(10)
1069  << ( (mass(1000022) > 1e7) ? scientific : fixed) << mass(1000022)
1070  << fixed << " ";
1071  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(4,icur) << " ";
1072 
1073  cout << endl << " | 1000023 " << setw(10)
1074  << ( (mass(1000023) > 1e7) ? scientific : fixed) << mass(1000023)
1075  << fixed << " ";
1076  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(5,icur) << " ";
1077 
1078  cout << endl << " | 1000025 " << setw(10)
1079  << ( (mass(1000025) > 1e7) ? scientific : fixed) << mass(1000025)
1080  << fixed << " ";
1081  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(6,icur) << " ";
1082 
1083  cout << endl << " | 1000035 " << setw(10)
1084  << ( (mass(1000035) > 1e7) ? scientific : fixed) << mass(1000035)
1085  << fixed << " ";
1086  for (int icur=1;icur<=7;icur++) cout << setw(6) << rvnmix(7,icur) << " ";
1087  }
1088  cout << endl;
1089 
1090  // Charged fermions (charginos)
1091  message(0,"","");
1092 
1093  // R-conserving:
1094  if (modsel(4) < 1) {
1095  cout << " | ~chi+ m U: |~W| |~H| ; V: |~W| |~H|"
1096  << endl;
1097 
1098  cout << setprecision(3) << " | 1000024 " << setw(10)
1099  << ((mass(1000024) > 1e7) ? scientific : fixed) << mass(1000024)
1100  << fixed << " ";
1101  for (int icur=1;icur<=2;icur++)
1102  cout << setw(6) << sqrt(pow2(umix(1,icur))+pow2(imumix(1,icur))) << " ";
1103  cout << "; ";
1104  for (int icur=1;icur<=2;icur++)
1105  cout << setw(6) << sqrt(pow2(vmix(1,icur))+pow2(imvmix(1,icur))) << " ";
1106 
1107  cout << endl << " | 1000037 " << setw(10)
1108  << ((mass(1000037) > 1e7) ? scientific : fixed) << mass(1000037)
1109  << fixed << " ";
1110  for (int icur=1;icur<=2;icur++)
1111  cout << setw(6) << sqrt(pow2(umix(2,icur))+pow2(imumix(2,icur))) << " ";
1112  cout << "; " ;
1113  for (int icur=1;icur<=2;icur++)
1114  cout << setw(6) << sqrt(pow2(vmix(2,icur))+pow2(imvmix(2,icur))) << " ";
1115  }
1116 
1117  // R-violating
1118  else {
1119  cout << " | e+/~chi+ m U: eL+ muL+ tauL+ ~W+"
1120  << " ~H1+ | V: eR+ muR+ tauR+ ~W+ ~H2+" << endl;
1121 
1122  cout << setprecision(3) << " | -11 " << setw(10)
1123  << ((mass(11) > 1e7) ? scientific : fixed) << mass(11)
1124  << fixed << " ";
1125  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvumix(1,icur) << " ";
1126  cout << "| ";
1127  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvvmix(1,icur) << " ";
1128 
1129  cout << endl << " | -13 " << setw(10)
1130  << ((mass(13) > 1e7) ? scientific : fixed) << mass(13)
1131  << fixed << " ";
1132  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvumix(2,icur) << " ";
1133  cout << "| " ;
1134  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvvmix(2,icur) << " ";
1135 
1136  cout << endl << " | -15 " << setw(10)
1137  << ((mass(15) > 1e7) ? scientific : fixed) << mass(15)
1138  << fixed << " ";
1139  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvumix(3,icur) << " ";
1140  cout << "| " ;
1141  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvvmix(3,icur) << " ";
1142 
1143  cout << endl << " | 1000024 " << setw(10)
1144  << ((mass(1000024) > 1e7) ? scientific : fixed) << mass(1000024)
1145  << fixed << " ";
1146  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvumix(4,icur) << " ";
1147  cout << "| " ;
1148  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvvmix(4,icur) << " ";
1149 
1150  cout << endl << " | 1000037 " << setw(10)
1151  << ((mass(1000037) > 1e7) ? scientific : fixed) << mass(1000037)
1152  << fixed << " ";
1153  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvumix(5,icur) << " ";
1154  cout << "| " ;
1155  for (int icur=1;icur<=5;icur++) cout << setw(6) << rvvmix(5,icur) << " ";
1156  }
1157  cout << endl;
1158 
1159  // Higgs bosons
1160  message(0,"","");
1161 
1162  // NMSSM
1163  if (modsel(3) >= 1) {
1164  cout << " | H m";
1165  if (nmhmix.exists()) cout << " H_10 H_20 S0";
1166  cout << endl;
1167 
1168  cout << setprecision(3) << " | 25 " << setw(10)
1169  << ( (mass(25) > 1e7) ? scientific : fixed) << mass(25)
1170  << fixed << " ";
1171  if (nmhmix.exists()) for (int icur=1;icur<=3;icur++)
1172  cout << setw(6) << nmhmix(1,icur) << " ";
1173 
1174  cout << endl << " | 35 " << setw(10)
1175  << ( (mass(35) > 1e7) ? scientific : fixed) << mass(35)
1176  << fixed << " ";
1177  if (nmhmix.exists()) for (int icur=1;icur<=3;icur++)
1178  cout << setw(6) << nmhmix(2,icur) << " ";
1179 
1180  cout << endl << " | 45 " << setw(10)
1181  << ( (mass(45) > 1e7) ? scientific : fixed) << mass(45)
1182  << fixed << " ";
1183  if (nmhmix.exists()) for (int icur=1;icur<=3;icur++)
1184  cout << setw(6) << nmhmix(3,icur) << " ";
1185 
1186  cout << endl <<" |"<<endl;
1187  cout << " | A m";
1188  if (nmamix.exists()) cout << " H_10 H_20 S0";
1189  cout << endl;
1190 
1191  cout << setprecision(3) << " | 36 " << setw(10)
1192  << ( (mass(36) > 1e7) ? scientific : fixed) << mass(36)
1193  << fixed << " ";
1194  if (nmamix.exists()) for (int icur=1;icur<=3;icur++)
1195  cout << setw(6) << nmamix(1,icur) << " ";
1196 
1197  cout << endl << " | 46 " << setw(10)
1198  << ( (mass(46) > 1e7) ? scientific : fixed) << mass(46)
1199  << fixed << " ";
1200  if (nmamix.exists()) for (int icur=1;icur<=3;icur++)
1201  cout << setw(6) << nmamix(2,icur) << " ";
1202 
1203  cout << endl <<" |"<<endl;
1204  cout << " | H+ m"<< endl;
1205 
1206  cout << setprecision(3) << " | 37 " << setw(10)
1207  << ( (mass(37) > 1e7) ? scientific : fixed) << mass(37)<<endl;
1208 
1209  cout << endl<<" |"<<endl;
1210  }
1211  // R-conserving MSSM (R-violating case handled above, with sneutrinos)
1212  else if (modsel(4) < 1) {
1213  cout << " | Higgs m"<<endl;
1214  cout << setprecision(3) << " | 25 " << setw(10)
1215  << ( (mass(25) > 1e7) ? scientific : fixed) << mass(25)<<endl;
1216  cout << setprecision(3) << " | 35 " << setw(10)
1217  << ( (mass(35) > 1e7) ? scientific : fixed) << mass(35)<<endl;
1218  cout << setprecision(3) << " | 36 " << setw(10)
1219  << ( (mass(36) > 1e7) ? scientific : fixed) << mass(36)<<endl;
1220  cout << setprecision(3) << " | 37 " << setw(10)
1221  << ( (mass(37) > 1e7) ? scientific : fixed) << mass(37)<<endl;
1222  cout << " |"<<endl;
1223  cout << " | alpha ";
1224  if (alpha.exists()) cout << setw(8) << alpha();
1225  else cout << " absent";
1226  cout << endl<<" |"<<endl;
1227  }
1228  // Running Higgs parameters
1229  if (hmix.exists()) {
1230  cout << " | Hmix "<<endl;
1231  cout << " | mu ";
1232  if (hmix.exists(1)) cout << setw(8) << hmix(1)
1233  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1234  else cout << " absent";
1235  cout << "\n | tan(beta) ";
1236  if (hmix.exists(2)) cout << setw(8) << hmix(2)
1237  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1238  else cout << " absent";
1239  cout << "\n | v ";
1240  if (hmix.exists(3)) cout << setw(8) << hmix(3)
1241  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1242  else cout << " absent";
1243  cout << "\n | mA ";
1244  if (hmix.exists(4)) cout << setw(9)
1245  << ((abs(hmix(4)) > 1e5) ? scientific : fixed) << hmix(4)
1246  << " (DRbar running value at Q = " << fixed << hmix.q() << " GeV)";
1247  else cout << " absent";
1248  cout << "\n";
1249  }
1250 
1251  // Gauge
1252  message(0,"","");
1253  if (gauge.exists()) {
1254  cout << " | Gauge " << endl;
1255  cout << " | g' ";
1256  if (gauge.exists(1)) cout << setw(8) << gauge(1)
1257  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1258  else cout << " absent";
1259  cout << "\n | g ";
1260  if (gauge.exists(2)) cout << setw(8) << gauge(2)
1261  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1262  else cout << " absent";
1263  cout << "\n | g3 ";
1264  if (gauge.exists(3)) cout << setw(8) << gauge(3)
1265  << " (DRbar running value at Q = " << hmix.q() << " GeV)";
1266  else cout << " absent";
1267  cout << "\n";
1268  }
1269 
1270  // Print footer
1271  footerPrinted=false;
1272  message(0,"","");
1273  listFooter();
1274 }
1275 
1276 //--------------------------------------------------------------------------
1277 
1278 // Check consistency of spectrum, unitarity of matrices, etc.
1279 // Return codes:
1280 // 2 : no SUSY spectrum, and no MASS block either.
1281 // 1 : no SUSY spectrum; but there is a MASS block.
1282 // 0 : there is a SUSY spectrum and it looks ok.
1283 // -1 : there is a SUSY spectrum but it does not look ok.
1284 
1285 int SusyLesHouches::checkSpectrum() {
1286 
1287  if (! headerPrinted) listHeader();
1288  int ifail=0;
1289  bool foundModsel = modsel.exists();
1290  if (! foundModsel) {
1291  if (mass.exists()) return 1;
1292  else return 2;
1293  }
1294 
1295  // Step 1) Check MODSEL. Assign default values where applicable.
1296  if (!modsel.exists(1)) {
1297  message(1,"checkSpectrum","MODSEL(1) undefined. Assuming = 0",0);
1298  modsel.set(1,0);
1299  }
1300  if (!modsel.exists(3)) modsel.set(3,0);
1301  if (!modsel.exists(4)) modsel.set(4,0);
1302  if (!modsel.exists(5)) modsel.set(5,0);
1303  if (!modsel.exists(6)) modsel.set(6,0);
1304  if (!modsel.exists(11)) modsel.set(11,1);
1305 
1306  // Step 2) Check for existence / duplication of blocks
1307 
1308  //Global
1309  if (!minpar.exists()) {
1310  message(1,"checkSpectrum","MINPAR not found",0);
1311  }
1312  if (!sminputs.exists()) {
1313  message(1,"checkSpectrum","SMINPUTS not found",0);
1314  }
1315  if (!mass.exists()) {
1316  message(1,"checkSpectrum","MASS not found",0);
1317  }
1318  if (!gauge.exists()) {
1319  message(1,"checkSpectrum","GAUGE not found",0);
1320  }
1321 
1322  //SLHA1
1323  if (modsel(3) == 0 && modsel(4) == 0 && modsel(5) == 0 && modsel(6) == 0) {
1324  // Check for required SLHA1 blocks
1325  if (!staumix.exists() && !selmix.exists()) {
1326  message(1,"checkSpectrum","STAUMIX or SELMIX not found",0);
1327  };
1328  if (!sbotmix.exists() && !dsqmix.exists()) {
1329  message(1,"checkSpectrum","SBOTMIX or DSQMIX not found",0);
1330  };
1331  if (!stopmix.exists() && !usqmix.exists()) {
1332  message(1,"checkSpectrum","STOPMIX or USQMIX not found",0);
1333  };
1334  if (!nmix.exists()) {
1335  message(1,"checkSpectrum","NMIX not found",0);
1336  };
1337  if (!umix.exists()) {
1338  message(1,"checkSpectrum","UMIX not found",0);
1339  };
1340  if (!vmix.exists()) {
1341  message(1,"checkSpectrum","VMIX not found",0);
1342  };
1343  if (modsel(3) == 0 && !alpha.exists()) {
1344  message(1,"checkSpectrum","ALPHA not found",0);
1345  }
1346  if (!hmix.exists()) {
1347  message(1,"checkSpectrum","HMIX not found",0);
1348  }
1349  if (!msoft.exists()) {
1350  message(1,"checkSpectrum","MSOFT not found",0);
1351  }
1352  }
1353 
1354  //RPV (+ FLV)
1355  else if (modsel(4) != 0) {
1356  // Check for required SLHA2 blocks (or see if can be extracted from SLHA1)
1357  if (!rvnmix.exists()) {
1358  if (nmix.exists()) {
1359  message(1,"checkSpectrum",
1360  "MODSEL 4 != 0 but NMIX given instead of RVNMIX",0);
1361  for (int i=1; i<=4; i++) {
1362  if (i<=3) rvnmix.set(i,i,1.0);
1363  for (int j=1; j<=4; j++)
1364  rvnmix.set(i+3,j+3,nmix(i,j));
1365  }
1366  } else {
1367  message(1,"checkSpectrum","MODSEL 4 != 0 but RVNMIX not found",0);
1368  ifail=-1;
1369  }
1370  }
1371  if (!rvumix.exists()) {
1372  if (umix.exists()) {
1373  message(1,"checkSpectrum",
1374  "MODSEL 4 != 0 but UMIX given instead of RVUMIX",0);
1375  for (int i=1; i<=3; i++) rvumix.set(i,i,1.0);
1376  for (int i=1; i<=2; i++) {
1377  for (int j=1; j<=2; j++)
1378  rvumix.set(i+3,j+3,umix(i,j));
1379  }
1380  } else {
1381  message(1,"checkSpectrum","MODSEL 4 != 0 but RVUMIX not found",0);
1382  ifail=-1;
1383  }
1384  }
1385  if (!rvvmix.exists()) {
1386  if (vmix.exists()) {
1387  message(1,"checkSpectrum",
1388  "MODSEL 4 != 0 but VMIX given instead of RVVMIX",0);
1389  for (int i=1; i<=3; i++) rvvmix.set(i,i,1.0);
1390  for (int i=1; i<=2; i++) {
1391  for (int j=1; j<=2; j++)
1392  rvvmix.set(i+3,j+3,vmix(i,j));
1393  }
1394  } else {
1395  message(1,"checkSpectrum","MODSEL 4 != 0 but RVVMIX not found",0);
1396  ifail=-1;
1397  }
1398  }
1399  if (!rvhmix.exists()) {
1400  if (alpha.exists()) {
1401  message(1,"checkSpectrum",
1402  "MODSEL 4 != 0 but ALPHA given instead of RVHMIX",0);
1403  rvhmix.set(1,1,cos(alpha()));
1404  rvhmix.set(1,2,sin(alpha()));
1405  rvhmix.set(2,1,-sin(alpha()));
1406  rvhmix.set(2,2,cos(alpha()));
1407  rvhmix.set(3,3,1.0);
1408  rvhmix.set(4,4,1.0);
1409  rvhmix.set(5,5,1.0);
1410  } else {
1411  message(1,"checkSpectrum","MODSEL 4 != 0 but RVHMIX not found",0);
1412  ifail=-1;
1413  }
1414  }
1415  if (!rvamix.exists()) {
1416  message(1,"checkSpectrum","MODSEL 4 != 0 but RVAMIX not found",0);
1417  }
1418  if (!rvlmix.exists()) {
1419  if (selmix.exists()) {
1420  message(1,"checkSpectrum",
1421  "MODSEL 4 != 0 but SELMIX given instead of RVLMIX",0);
1422  for (int i=1; i<=6; i++) {
1423  for (int j=6; j<=6; j++)
1424  rvlmix.set(i+1,j+2,selmix(i,j));
1425  }
1426  } if (staumix.exists()) {
1427  message(1,"checkSpectrum",
1428  "MODSEL 4 != 0 but STAUMIX given instead of RVLMIX",0);
1429  rvlmix.set(2,3,1.0);
1430  rvlmix.set(3,4,1.0);
1431  rvlmix.set(4,5,staumix(1,1));
1432  rvlmix.set(4,8,staumix(1,2));
1433  rvlmix.set(5,6,1.0);
1434  rvlmix.set(6,7,1.0);
1435  rvlmix.set(7,5,staumix(2,1));
1436  rvlmix.set(7,8,staumix(2,2));
1437  } else {
1438  message(1,"checkSpectrum","MODSEL 4 != 0 but RVLMIX not found",0);
1439  ifail=-1;
1440  }
1441  }
1442  if (!usqmix.exists()) {
1443  if (stopmix.exists()) {
1444  message(1,"checkSpectrum",
1445  "MODSEL 4 != 0 but STOPMIX given instead of USQMIX",0);
1446  usqmix.set(1,1, 1.0);
1447  usqmix.set(2,2, 1.0);
1448  usqmix.set(4,4, 1.0);
1449  usqmix.set(5,5, 1.0);
1450  usqmix.set(3,3, stopmix(1,1));
1451  usqmix.set(3,6, stopmix(1,2));
1452  usqmix.set(6,3, stopmix(2,1));
1453  usqmix.set(6,6, stopmix(2,2));
1454  } else {
1455  message(1,"checkSpectrum","MODSEL 4 != 0 but USQMIX not found",0);
1456  ifail=-1;
1457  }
1458  }
1459  if (!dsqmix.exists()) {
1460  if (sbotmix.exists()) {
1461  message(1,"checkSpectrum",
1462  "MODSEL 4 != 0 but SBOTMIX given instead of DSQMIX",0);
1463  dsqmix.set(1,1, 1.0);
1464  dsqmix.set(2,2, 1.0);
1465  dsqmix.set(4,4, 1.0);
1466  dsqmix.set(5,5, 1.0);
1467  dsqmix.set(3,3, sbotmix(1,1));
1468  dsqmix.set(3,6, sbotmix(1,2));
1469  dsqmix.set(6,3, sbotmix(2,1));
1470  dsqmix.set(6,6, sbotmix(2,2));
1471  } else {
1472  message(1,"checkSpectrum","MODSEL 4 != 0 but DSQMIX not found",0);
1473  ifail=-1;
1474  }
1475  }
1476  }
1477 
1478  // FLV but not RPV (see above for FLV+RPV, below for FLV regardless of RPV)
1479  else if (modsel(6) != 0) {
1480  // Quark FLV
1481  if (modsel(6) != 2) {
1482  if (!usqmix.exists()) {
1483  message(1,"checkSpectrum","quark FLV on but USQMIX not found",0);
1484  ifail=-1;
1485  }
1486  if (!dsqmix.exists()) {
1487  message(1,"checkSpectrum","quark FLV on but DSQMIX not found",0);
1488  ifail=-1;
1489  }
1490  }
1491  // Lepton FLV
1492  if (modsel(6) != 1) {
1493  if (!upmns.exists()) {
1494  message(1,"checkSpectrum","lepton FLV on but UPMNSIN not found",0);
1495  ifail=-1;
1496  }
1497  if (!selmix.exists()) {
1498  message(1,"checkSpectrum","lepton FLV on but SELMIX not found",0);
1499  ifail=-1;
1500  }
1501  if (!snumix.exists() && !snsmix.exists()) {
1502  message(1,"checkSpectrum","lepton FLV on but SNUMIX not found",0);
1503  ifail=-1;
1504  }
1505  }
1506  }
1507 
1508  // CPV
1509  if (modsel(5) != 0) {
1510  if (!cvhmix.exists()) {
1511  message(1,"checkSpectrum","MODSEL 5 != 0 but CVHMIX not found",0);
1512  ifail=-1;
1513  }
1514  }
1515 
1516  // FLV (regardless of whether RPV or not)
1517  if (modsel(6) != 0) {
1518  // Quark FLV
1519  if (modsel(6) != 2) {
1520  if (!vckmin.exists()) {
1521  message(1,"checkSpectrum","quark FLV on but VCKMIN not found",0);
1522  ifail=-1;
1523  }
1524  // The following checks are not considered crucial since the entries
1525  // *may* not be needed to compute required derived couplings.
1526  if (!msq2in.exists()) {
1527  message(0,"checkSpectrum","note: quark FLV on but MSQ2IN not found",0);
1528  }
1529  if (!msu2in.exists()) {
1530  message(0,"checkSpectrum","note: quark FLV on but MSU2IN not found",0);
1531  }
1532  if (!msd2in.exists()) {
1533  message(0,"checkSpectrum","note: quark FLV on but MSD2IN not found",0);
1534  }
1535  if (!tuin.exists()) {
1536  message(0,"checkSpectrum","note: quark FLV on but TUIN not found",0);
1537  }
1538  if (!tdin.exists()) {
1539  message(0,"checkSpectrum","note: quark FLV on but TDIN not found",0);
1540  }
1541  }
1542  // Lepton FLV
1543  if (modsel(6) != 1) {
1544  if (!msl2in.exists()) {
1545  message(0,"checkSpectrum",
1546  "note: lepton FLV on but MSL2IN not found",0);
1547  }
1548  if (!mse2in.exists()) {
1549  message(0,"checkSpectrum",
1550  "note: lepton FLV on but MSE2IN not found",0);
1551  }
1552  if (!tein.exists()) {
1553  message(0,"checkSpectrum",
1554  "note: lepton FLV on but TEIN not found",0);
1555  }
1556  }
1557  }
1558 
1559  // Step 3) SLHA1 --> SLHA2 interoperability
1560  //Note: the mass basis is NOT mass-ordered in SLHA1, so be careful!
1561  //Here, the mass basis is hence by PDG code, not by mass-ordered value.
1562 
1563  if (stopmix.exists() && ! usqmix.exists() ) {
1564  //1000002 = ~uL, 1000004 = ~cL, 2000002 = ~uR, 2000004 = ~cR
1565  usqmix.set(1,1, 1.0);
1566  usqmix.set(2,2, 1.0);
1567  usqmix.set(4,4, 1.0);
1568  usqmix.set(5,5, 1.0);
1569  //Fill (1000006,2000006) sector from stopmix
1570  usqmix.set(3,3, stopmix(1,1));
1571  usqmix.set(3,6, stopmix(1,2));
1572  usqmix.set(6,3, stopmix(2,1));
1573  usqmix.set(6,6, stopmix(2,2));
1574  };
1575  if (sbotmix.exists() && ! dsqmix.exists() ) {
1576  //1000001 = ~dL, 1000003 = ~sL, 2000001 = ~dR, 2000003 = ~sR
1577  dsqmix.set(1,1, 1.0);
1578  dsqmix.set(2,2, 1.0);
1579  dsqmix.set(4,4, 1.0);
1580  dsqmix.set(5,5, 1.0);
1581  //Fill (1000005,2000005) sector from sbotmix
1582  dsqmix.set(3,3, sbotmix(1,1));
1583  dsqmix.set(3,6, sbotmix(1,2));
1584  dsqmix.set(6,3, sbotmix(2,1));
1585  dsqmix.set(6,6, sbotmix(2,2));
1586  };
1587  if (staumix.exists() && ! selmix.exists() ) {
1588  //1000011 = ~eL, 1000013 = ~muL, 2000011 = ~eR, 2000013 = ~muR
1589  selmix.set(1,1, 1.0);
1590  selmix.set(2,2, 1.0);
1591  selmix.set(4,4, 1.0);
1592  selmix.set(5,5, 1.0);
1593  //Fill (1000015,2000015) sector from staumix
1594  selmix.set(3,3, staumix(1,1));
1595  selmix.set(3,6, staumix(1,2));
1596  selmix.set(6,3, staumix(2,1));
1597  selmix.set(6,6, staumix(2,2));
1598  };
1599  if (! snumix.exists() && ! snsmix.exists()) {
1600  //1000012 = ~nu_e, 1000014 = ~nu_mu, 1000016 = ~nu_tau
1601  snumix.set(1,1, 1.0);
1602  snumix.set(2,2, 1.0);
1603  snumix.set(3,3, 1.0);
1604  };
1605 
1606  // Step 4) Check mass ordering and unitarity/orthogonality of mixing matrices
1607 
1608  // Check expected mass orderings
1609  if (mass.exists()) {
1610  // CP-even Higgs
1611  if (abs(mass(25)) > abs(mass(35))
1612  || (modsel(3) == 1 && abs(mass(35)) > abs(mass(45))) )
1613  message(0,"checkSpectrum","Note: Higgs sector is not mass-ordered",0);
1614  // CP-odd Higgs
1615  if (modsel(3) == 1 && abs(mass(36)) > abs(mass(46)))
1616  message(0,"checkSpectrum",
1617  "Note: CP-odd Higgs sector is not mass-ordered",0);
1618  // Neutralinos
1619  if (abs(mass(1000022)) > abs(mass(1000023))
1620  || abs(mass(1000023)) > abs(mass(1000025))
1621  || abs(mass(1000025)) > abs(mass(1000035))
1622  || (modsel(3) == 1 && abs(mass(1000035)) > abs(mass(1000045))) )
1623  message(0,"checkSpectrum","Note: Neutralino sector is not mass-ordered"
1624  ,0);
1625  // Charginos
1626  if (abs(mass(1000024)) > abs(mass(1000037)))
1627  message(0,"checkSpectrum","Note: Chargino sector is not mass-ordered",0);
1628  }
1629 
1630  //NMIX
1631  if (nmix.exists()) {
1632  for (int i=1;i<=4;i++) {
1633  double cn1=0.0;
1634  double cn2=0.0;
1635  for (int j=1;j<=4;j++) {
1636  cn1 += pow(nmix(i,j),2);
1637  cn2 += pow(nmix(j,i),2);
1638  if (imnmix.exists()) {
1639  cn1 += pow(imnmix(i,j),2);
1640  cn2 += pow(imnmix(j,i),2);
1641  }
1642  }
1643  if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) {
1644  ifail=-1;
1645  message(2,"checkSpectrum","NMIX is not unitary (wrong format?)",0);
1646  break;
1647  }
1648  }
1649  }
1650 
1651  //VMIX, UMIX
1652  if (vmix.exists() && umix.exists()) {
1653  // First check for non-standard "madgraph" convention
1654  // (2,2) entry not given explicitly
1655  for (int i=1;i<=2;i++) {
1656  double cu1=0.0;
1657  double cu2=0.0;
1658  double cv1=0.0;
1659  double cv2=0.0;
1660  for (int j=1;j<=2;j++) {
1661  cu1 += pow(umix(i,j),2);
1662  cu2 += pow(umix(j,i),2);
1663  cv1 += pow(vmix(i,j),2);
1664  cv2 += pow(vmix(j,i),2);
1665  if (imumix.exists()) {
1666  cu1 += pow(imumix(i,j),2);
1667  cu2 += pow(imumix(j,i),2);
1668  }
1669  if (imvmix.exists()) {
1670  cv1 += pow(imvmix(i,j),2);
1671  cv2 += pow(imvmix(j,i),2);
1672  }
1673  }
1674  if (abs(1.0-cu1) > 1e-3 || abs(1.0-cu2) > 1e-3) {
1675  cu1 += pow(umix(1,1),2);
1676  cu2 += pow(umix(1,1),2);
1677  if (abs(1.0-cu1) > 1e-3 || abs(1.0-cu2) > 1e-3) {
1678  ifail=-1;
1679  message(2,"checkSpectrum","UMIX is not unitary (wrong format?)",0);
1680  break;
1681  } else {
1682  // Fix madgraph non-standard convention problem
1683  message(1,"checkSpectrum","UMIX is not unitary (repaired)",0);
1684  umix.set(2,2,umix(1,1));
1685  }
1686  }
1687  if (abs(1.0-cv1) > 1e-3 || abs(1.0-cv2) > 1e-3) {
1688  cv1 += pow(vmix(1,1),2);
1689  cv2 += pow(vmix(1,1),2);
1690  if (abs(1.0-cv1) > 1e-3 || abs(1.0-cv2) > 1e-3) {
1691  ifail=-1;
1692  message(2,"checkSpectrum","VMIX is not unitary (wrong format?)",0);
1693  break;
1694  } else {
1695  // Fix madgraph non-standard convention problem
1696  message(1,"checkSpectrum","VMIX is not unitary (repaired)",0);
1697  vmix.set(2,2,vmix(1,1));
1698  }
1699  }
1700  }
1701 
1702  }
1703 
1704  //STOPMIX, SBOTMIX
1705  if (stopmix.exists() && sbotmix.exists()) {
1706  for (int i=1;i<=2;i++) {
1707  double ct1=0.0;
1708  double ct2=0.0;
1709  double cb1=0.0;
1710  double cb2=0.0;
1711  for (int j=1;j<=2;j++) {
1712  ct1 += pow(stopmix(i,j),2);
1713  ct2 += pow(stopmix(j,i),2);
1714  cb1 += pow(sbotmix(i,j),2);
1715  cb2 += pow(sbotmix(j,i),2);
1716  }
1717  if (abs(1.0-ct1) > 1e-3 || abs(1.0-ct2) > 1e-3) {
1718  ifail=-1;
1719  message(2,"checkSpectrum","STOPMIX is not unitary (wrong format?)",0);
1720  break;
1721  }
1722  if (abs(1.0-cb1) > 1e-3 || abs(1.0-cb2) > 1e-3) {
1723  ifail=-1;
1724  message(2,"checkSpectrum","SBOTMIX is not unitary (wrong format?)",0);
1725  break;
1726  }
1727  }
1728  }
1729 
1730  //STAUMIX
1731  if (staumix.exists()) {
1732  for (int i=1;i<=2;i++) {
1733  double ct1=0.0;
1734  double ct2=0.0;
1735  for (int j=1;j<=2;j++) {
1736  ct1 += pow(staumix(i,j),2);
1737  ct2 += pow(staumix(j,i),2);
1738  }
1739  if (abs(1.0-ct1) > 1e-3 || abs(1.0-ct2) > 1e-3) {
1740  ifail=-1;
1741  message(2,"checkSpectrum","STAUMIX is not unitary (wrong format?)",0);
1742  break;
1743  }
1744  }
1745  }
1746 
1747  //DSQMIX
1748  if (dsqmix.exists()) {
1749  for (int i=1;i<=6;i++) {
1750  double sr=0.0;
1751  double sc=0.0;
1752  for (int j=1;j<=6;j++) {
1753  sr += pow(dsqmix(i,j),2);
1754  sc += pow(dsqmix(j,i),2);
1755  if (imdsqmix.exists()) {
1756  sr += pow(imdsqmix(i,j),2);
1757  sc += pow(imdsqmix(j,i),2);
1758  }
1759  }
1760  if (abs(1.0-sr) > 1e-3 || abs(1.0-sc) > 1e-3) {
1761  ifail=-1;
1762  message(2,"checkSpectrum","DSQMIX is not unitary (wrong format?)",0);
1763  break;
1764  }
1765  }
1766  }
1767 
1768  //USQMIX
1769  if (usqmix.exists()) {
1770  for (int i=1;i<=6;i++) {
1771  double sr=0.0;
1772  double sc=0.0;
1773  for (int j=1;j<=6;j++) {
1774  sr += pow(usqmix(i,j),2);
1775  sc += pow(usqmix(j,i),2);
1776  if (imusqmix.exists()) {
1777  sr += pow(imusqmix(i,j),2);
1778  sc += pow(imusqmix(j,i),2);
1779  }
1780  }
1781  if (abs(1.0-sr) > 1e-3 || abs(1.0-sc) > 1e-3) {
1782  ifail=-1;
1783  message(2,"checkSpectrum","USQMIX is not unitary (wrong format?)",0);
1784  break;
1785  }
1786  }
1787  }
1788 
1789  //SELMIX
1790  if (selmix.exists()) {
1791  for (int i=1;i<=6;i++) {
1792  double sr=0.0;
1793  double sc=0.0;
1794  for (int j=1;j<=6;j++) {
1795  sr += pow(selmix(i,j),2);
1796  sc += pow(selmix(j,i),2);
1797  if (imselmix.exists()) {
1798  sr += pow(imselmix(i,j),2);
1799  sc += pow(imselmix(j,i),2);
1800  }
1801  }
1802  if (abs(1.0-sr) > 1e-3 || abs(1.0-sc) > 1e-3) {
1803  ifail=-1;
1804  message(2,"checkSpectrum","SELMIX is not unitary (wrong format?)",0);
1805  break;
1806  }
1807  }
1808  } //NMSSM:
1809  if (modsel(3) == 1) {
1810  //NMNMIX
1811  if ( nmnmix.exists() ) {
1812  for (int i=1;i<=5;i++) {
1813  double cn1=0.0;
1814  double cn2=0.0;
1815  for (int j=1;j<=5;j++) {
1816  cn1 += pow(nmnmix(i,j),2);
1817  cn2 += pow(nmnmix(j,i),2);
1818  if (imnmnmix.exists()) {
1819  cn1 += pow(imnmnmix(i,j),2);
1820  cn2 += pow(imnmnmix(j,i),2);
1821  }
1822  }
1823  if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) {
1824  ifail=-1;
1825  message(2,"checkSpectrum","NMNMIX is not unitary (wrong format?)",0);
1826  break;
1827  }
1828  }
1829  }
1830  else {
1831  ifail=-1;
1832  message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMNMIX found",0);
1833  }
1834  //NMAMIX
1835  if ( nmamix.exists() ) {
1836  for (int i=1;i<=2;i++) {
1837  double cn1=0.0;
1838  for (int j=1;j<=3;j++) {
1839  cn1 += pow(nmamix(i,j),2);
1840  }
1841  if (abs(1.0-cn1) > 1e-3) {
1842  ifail=-1;
1843  message(2,"checkSpectrum","NMAMIX is not unitary (wrong format?)",0);
1844  break;
1845  }
1846  }
1847  }
1848  else {
1849  ifail=-1;
1850  message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMAMIX found",0);
1851  }
1852  //NMHMIX
1853  if ( nmhmix.exists() ) {
1854  for (int i=1;i<=3;i++) {
1855  double cn1=0.0;
1856  double cn2=0.0;
1857  for (int j=1;j<=3;j++) {
1858  cn1 += pow(nmhmix(i,j),2);
1859  cn2 += pow(nmhmix(j,i),2);
1860  }
1861  if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) {
1862  ifail=-1;
1863  message(2,"checkSpectrum","NMHMIX is not unitary (wrong format?)",0);
1864  break;
1865  }
1866  }
1867  }
1868  else {
1869  ifail=-1;
1870  message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMHMIX found",0);
1871  }
1872  //NMSSMRUN
1873  if (! nmssmrun.exists() ) {
1874  ifail=-1;
1875  message(2,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMSSMRUN found",
1876  0);
1877  }
1878  }
1879 
1880  //Check for documentation
1881  if (slhaRead && ! spinfo.exists(1)) spinfo.set(1,"unknown");
1882  if (slhaRead && ! spinfo.exists(2)) spinfo.set(2,"unknown");
1883  if (! slhaRead && ! spinfo.exists(1)) {
1884  spinfo.set(1,"DEFAULT");
1885  spinfo.set(2,"n/a");
1886  }
1887 
1888  //Give status
1889  if (ifail >= 2)
1890  message(0,"checkSpectrum","one or more serious problems were found");
1891  else if (ifail == -1)
1892  message(0,"checkSpectrum","one or more inconsistencies in SUSY setup");
1893 
1894  //Print Footer
1895  listFooter();
1896 
1897  //Return
1898  return ifail;
1899 }
1900 
1901 //--------------------------------------------------------------------------
1902 
1903 // Simple utility to print messages, warnings, and errors
1904 
1905 void SusyLesHouches::message(int level, string place,string themessage,
1906  int line) {
1907  if (verboseSav == 0) return;
1908  // By default all output to cout, but lines below allow finer control.
1909  ostream* outstream = &cout;
1910  //Send normal messages and warnings to stdout, errors to stderr.
1911  //ostream* outstream = &cerr;
1912  //if (level <= 1) outstream = &cout;
1913  // if (level == 2) { *outstream << endl; }
1914  if (place != "") *outstream << " | (SLHA::"+place+") ";
1915  else *outstream << " | ";
1916  if (level == 1) *outstream << "Warning: ";
1917  if (level == 2) { *outstream << "ERROR: "; }
1918  if (line != 0) *outstream << "line " << line << " - ";
1919  *outstream << themessage << endl;
1920  // if (level == 2) *outstream << endl;
1921  footerPrinted=false;
1922  return;
1923 }
1924 
1925 //==========================================================================
1926 
1927 } // end namespace Pythia8
Generated by   doxygen 1.8.5