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