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