source: HiSusy/trunk/Pythia8/pythia8170/src/SusyLesHouches.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

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