- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / SusyLesHouches.cxx
1 // SusyLesHouches.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 Peter Skands, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6 #include "SusyLesHouches.h"
7
8 //==========================================================================
9
10 // Main routine to read in SLHA and LHEF+SLHA files
11
12 int SusyLesHouches::readFile(string slhaFileIn, int verboseIn) {
13
14   // Copy inputs to local
15   slhaFile = slhaFileIn;
16   verbose  = verboseIn;
17
18   // Check that input file is OK.
19   int iFailFile=0;
20   const char* cstring = slhaFile.c_str();
21   ifstream file(cstring);  
22   if (!file.good()) {
23     message(2,"readFile",slhaFile+" not found",0);
24     return -1;
25     slhaRead=false;
26   }  
27
28   if (verbose >= 3) message(0,"readFile","parsing "+slhaFile,0);
29
30   // Array of particles read in.
31   vector<int> idRead;
32
33   //Initial values for read-in variables.  
34   slhaRead=true;
35   lhefRead=false;
36   lhefSlha=false;
37   bool foundSlhaTag = false;
38   bool xmlComment   = false;
39   bool decayPrinted = false;
40   string line="";
41   string blockIn="";
42   string decay="";
43   string comment="";
44   string blockName="";
45   string nameNow="";
46   int idNow=0;
47   double width=0.0;
48
49   //Initialize line counter
50   int iLine=0;
51
52   // Read in one line at a time.
53   while ( getline(file, line) ) {
54     iLine++;
55
56     //Rewrite string in lowercase
57     for (unsigned int i=0;i<line.length();i++) line[i]=tolower(line[i]);
58
59     // Remove extra blanks 
60     while (line.find("  ") != string::npos) line.erase( line.find("  "), 1);
61
62     //Detect whether read-in is from a Les Houches Event File (LHEF).
63     if (line.find("<leshouches") != string::npos 
64         || line.find("<slha") != string::npos) {
65       lhefRead=true;
66     }
67
68     // If LHEF
69     if (lhefRead) {      
70       //Ignore XML comments (only works for whole lines so far)
71       if (line.find("-->") != string::npos) {
72         xmlComment = false;
73       }
74       else if (xmlComment) continue;
75       else if (line.find("<!--") != string::npos) {
76         xmlComment = true;
77       } 
78       //Detect when <slha> tag reached.
79       if (line.find("<slha") != string::npos) {
80         lhefSlha     = true;
81         foundSlhaTag = true;
82         //Print header if not already done
83         if (! headerPrinted) printHeader();
84       }
85       //Stop looking when </header> or <init> tag reached
86       if (line.find("</header>") != string::npos ||
87           line.find("<init") != string::npos) {      
88         if (!foundSlhaTag) return 101;
89         break;
90       }
91       //If <slha> tag not yet reached, skip
92       if (!lhefSlha) continue;
93     }
94
95     //Ignore comment lines with # as first character
96     if (line.find("#") == 0) continue;
97
98     //Move comment to separate string
99     if (line.find("#") != string::npos) {
100       comment=line.substr(line.find("#")+1,line.length()-line.find("#")-1);
101       line.erase(line.find("#"),line.length()-line.find("#"));
102     }
103
104     // Remove blanks before and after an = sign.
105     while (line.find(" =") != string::npos) line.erase( line.find(" ="), 1);
106     while (line.find("= ") != string::npos) line.erase( line.find("= ")+1, 1);
107
108     //New block. 
109     if (line.find("block") <= 1) { 
110
111       //Print header if not already done
112       if (! headerPrinted) printHeader();
113
114       blockIn=line ; 
115       decay="";
116       int nameBegin=6 ;
117       int nameEnd=blockIn.find(" ",7);
118       blockName=line.substr(nameBegin,nameEnd-nameBegin);
119       
120       //Find Q=... for DRbar running blocks
121       if (blockIn.find("q=") != string::npos) {
122         int qbegin=blockIn.find("q=")+2;
123         istringstream qstream(blockIn.substr(qbegin,blockIn.length()));
124         double q=0.0;
125         qstream >> q;
126         if (qstream) {
127           // SLHA1 running blocks
128           if (blockName=="hmix") hmix.setq(q);
129           if (blockName=="yu") yu.setq(q);
130           if (blockName=="yd") yd.setq(q);
131           if (blockName=="ye") ye.setq(q);
132           if (blockName=="au") au.setq(q);
133           if (blockName=="ad") ad.setq(q);
134           if (blockName=="ae") ae.setq(q);
135           if (blockName=="msoft") msoft.setq(q);
136           if (blockName=="gauge") gauge.setq(q);
137           // SLHA2 running blocks
138           if (blockName=="vckm") vckm.setq(q);
139           if (blockName=="upmns") upmns.setq(q);
140           if (blockName=="msq2") msq2.setq(q);
141           if (blockName=="msu2") msu2.setq(q);
142           if (blockName=="msd2") msd2.setq(q);
143           if (blockName=="msl2") msl2.setq(q);
144           if (blockName=="mse2") mse2.setq(q);
145           if (blockName=="tu") tu.setq(q);
146           if (blockName=="td") td.setq(q);
147           if (blockName=="te") te.setq(q);
148           if (blockName=="rvlamlle") rvlamlle.setq(q);
149           if (blockName=="rvlamlqd") rvlamlqd.setq(q);
150           if (blockName=="rvlamudd") rvlamudd.setq(q);
151           if (blockName=="rvtlle") rvtlle.setq(q);
152           if (blockName=="rvtlqd") rvtlqd.setq(q);
153           if (blockName=="rvtudd") rvtudd.setq(q);
154           if (blockName=="rvkappa") rvkappa.setq(q);
155           if (blockName=="rvd") rvd.setq(q);
156           if (blockName=="rvm2lh1") rvm2lh1.setq(q);
157           if (blockName=="rvsnvev") rvsnvev.setq(q);          
158           if (blockName=="imau") imau.setq(q);
159           if (blockName=="imad") imad.setq(q);
160           if (blockName=="imae") imae.setq(q);
161           if (blockName=="imhmix") imhmix.setq(q);
162           if (blockName=="immsoft") immsoft.setq(q);
163           if (blockName=="imtu") imtu.setq(q);
164           if (blockName=="imtd") imtd.setq(q);
165           if (blockName=="imte") imte.setq(q);
166           if (blockName=="imvckm") imvckm.setq(q);
167           if (blockName=="imupmns") imupmns.setq(q);
168           if (blockName=="immsq2") immsq2.setq(q);
169           if (blockName=="immsu2") immsu2.setq(q);
170           if (blockName=="immsd2") immsd2.setq(q);
171           if (blockName=="immsl2") immsl2.setq(q);
172           if (blockName=="immse2") immse2.setq(q);
173         };
174       };
175       
176       //Skip to next line.
177       continue ; 
178       
179     } 
180
181     //New decay table
182     else if (line.find("decay") <= 1) {
183
184       // Print header if not already done
185       if (! headerPrinted) printHeader();
186
187       // If previous had zero length, print now
188       if (decay != "" && ! decayPrinted) {
189         if (verbose >= 2) message(0,"readFile","reading  WIDTH for "+nameNow
190                 +" (but no decay channels found)",0);
191       }
192
193       //Set decay block name
194       decay=line;
195       blockIn="";
196       int nameBegin=6 ;
197       int nameEnd=decay.find(" ",7);
198       nameNow=decay.substr(nameBegin,nameEnd-nameBegin);
199       
200       //Extract PDG code and width
201       istringstream dstream(nameNow);
202       dstream >> idNow;
203       if (dstream) {
204         string widthName=decay.substr(nameEnd+1,decay.length());
205         istringstream wstream(widthName);
206         wstream >> width;
207         if (wstream) {
208           // Set 
209           decayTable tmp(idNow,width);
210           decays.push_back(decayTable(idNow,width));          
211           decayIndices[idNow]=decays.size()-1;
212           //Set PDG code and width
213           if (width <= 0.0) {
214             string endComment="";
215             if (width < -1e-6) {
216               endComment="(forced width < 0 to zero)";
217             }
218             if (verbose >= 2)
219               message(0,"readFile","reading  stable particle "+nameNow
220                       +" "+endComment,0);
221             width=0.0;
222             decayPrinted = true;
223             decays[decayIndices[idNow]].setWidth(width);
224           } else {
225             decayPrinted = false;
226           }
227         } else {
228           if (verbose >= 2) 
229             message(0,"readFile","ignoring DECAY table for "+nameNow
230                     +" (read failed)",iLine);
231           decayPrinted = true;
232           width=0.0;
233           decay="";
234           continue;
235         }
236       }
237       else {
238         message(0,"readFile",
239                     "PDG Code unreadable. Ignoring this DECAY block",iLine);
240         decayPrinted = true;
241         decay="";
242         continue;
243       }
244
245       //Skip to next line
246       continue ;
247     }
248
249     //Switch off SLHA read-in via LHEF if outside <slha> tag.
250     else if (line.find("</slha>") != string::npos) {
251       lhefSlha=false;
252       blockIn="";
253       decay="";
254       continue;
255     }
256
257     //Skip not currently reading block data lines.
258     if (blockIn != "") {
259
260       // Replace an equal sign by a blank to make parsing simpler.
261       while (line.find("=") != string::npos) {
262         int firstEqual = line.find_first_of("=");
263         line.replace(firstEqual, 1, " ");   
264       };
265     
266       //Parse data lines within given block
267       //Constructed explicitly so that each block can have its own types and
268       //own rules defined. For extra user blocks, just add more recognized 
269       //blockNames at the end and implement user defined rules accordingly.
270       //string comment = line.substr(line.find("#"),line.length());    
271       int ifail=-2;
272       istringstream linestream(line);
273       //MODEL
274       if (blockName == "modsel") {
275         int i;
276         linestream >> i; 
277         if (linestream) {
278           if (i == 12) {ifail=modsel12.set(0,linestream);} 
279           else if (i == 21) {ifail=modsel21.set(0,linestream);}
280           else {ifail=modsel.set(i,linestream);};}
281         else {
282           ifail = -1;}
283       };
284       if (blockName == "minpar") ifail=minpar.set(linestream); 
285       if (blockName == "sminputs") ifail=sminputs.set(linestream);
286       if (blockName == "extpar") ifail=extpar.set(linestream);
287       if (blockName == "qextpar") ifail=qextpar.set(linestream);
288       //FLV
289       if (blockName == "vckmin") ifail=vckmin.set(linestream);
290       if (blockName == "upmnsin") ifail=upmnsin.set(linestream);
291       if (blockName == "msq2in") ifail=msq2in.set(linestream);
292       if (blockName == "msu2in") ifail=msu2in.set(linestream);
293       if (blockName == "msd2in") ifail=msd2in.set(linestream);
294       if (blockName == "msl2in") ifail=msl2in.set(linestream);
295       if (blockName == "mse2in") ifail=mse2in.set(linestream);
296       if (blockName == "tuin") ifail=tuin.set(linestream);
297       if (blockName == "tdin") ifail=tdin.set(linestream);
298       if (blockName == "tein") ifail=tein.set(linestream);
299       //RPV
300       if (blockName == "rvlamllein") ifail=rvlamllein.set(linestream);
301       if (blockName == "rvlamlqdin") ifail=rvlamlqdin.set(linestream);
302       if (blockName == "rvlamuddin") ifail=rvlamuddin.set(linestream);
303       if (blockName == "rvtllein") ifail=rvtllein.set(linestream);
304       if (blockName == "rvtlqdin") ifail=rvtlqdin.set(linestream);
305       if (blockName == "rvtuddin") ifail=rvtuddin.set(linestream);
306       if (blockName == "rvkappain") ifail=rvkappain.set(linestream);
307       if (blockName == "rvdin") ifail=rvdin.set(linestream);
308       if (blockName == "rvm2lh1in") ifail=rvm2lh1in.set(linestream);
309       if (blockName == "rvsnvevin") ifail=rvsnvevin.set(linestream);
310       //CPV 
311       if (blockName == "imminpar") ifail=imminpar.set(linestream);
312       if (blockName == "imextpar") ifail=imextpar.set(linestream);
313       //CPV +FLV
314       if (blockName == "immsq2in") ifail=immsq2in.set(linestream);
315       if (blockName == "immsu2in") ifail=immsu2in.set(linestream);
316       if (blockName == "immsd2in") ifail=immsd2in.set(linestream);
317       if (blockName == "immsl2in") ifail=immsl2in.set(linestream);
318       if (blockName == "immse2in") ifail=immse2in.set(linestream);
319       if (blockName == "imtuin") ifail=imtuin.set(linestream);
320       if (blockName == "imtdin") ifail=imtdin.set(linestream);
321       if (blockName == "imtein") ifail=imtein.set(linestream);
322       //Info:
323       if (blockName == "spinfo" || blockName=="dcinfo") {
324         int i;
325         string entry;
326         linestream >> i >> entry;
327         string blockStr="RGE";
328         if (blockName=="dcinfo") blockStr="DCY";
329
330         if (linestream) {
331           if ( i == 3 ) {
332             string warning=line.substr(line.find("3")+1,line.length());
333             message(1,"readFile","(from "+blockStr+" program): "+warning,0);
334             if (blockName == "spinfo") spinfo3.set(warning);
335             else dcinfo3.set(warning);
336           } else if ( i == 4 ) {
337             string error=line.substr(line.find("4")+1,line.length());
338             message(2,"readFile","(from "+blockStr+" program): "+error,0);
339             if (blockName == "spinfo") spinfo4.set(error);
340             else dcinfo4.set(error);
341           } else {
342             //Rewrite string in uppercase
343             for (unsigned int j=0;j<entry.length();j++) 
344               entry[j]=toupper(entry[j]);
345             ifail=(blockName=="spinfo") ? spinfo.set(i,entry)
346               : dcinfo.set(i,entry);
347           };
348         } else {
349           ifail=-1;
350         };
351       };
352       //SPECTRUM
353       //Pole masses
354       if (blockName == "mass") ifail=mass.set(linestream);
355
356       //Mixing
357       if (blockName == "alpha") ifail=alpha.set(linestream);
358       if (blockName == "stopmix") ifail=stopmix.set(linestream);
359       if (blockName == "sbotmix") ifail=sbotmix.set(linestream);
360       if (blockName == "staumix") ifail=staumix.set(linestream);
361       if (blockName == "nmix") ifail=nmix.set(linestream);
362       if (blockName == "umix") ifail=umix.set(linestream);
363       if (blockName == "vmix") ifail=vmix.set(linestream);
364       //FLV
365       if (blockName == "usqmix") ifail=usqmix.set(linestream);
366       if (blockName == "dsqmix") ifail=dsqmix.set(linestream);
367       if (blockName == "selmix") ifail=selmix.set(linestream);
368       if (blockName == "snumix") ifail=snumix.set(linestream);
369       if (blockName == "snsmix") ifail=snsmix.set(linestream);
370       if (blockName == "snamix") ifail=snamix.set(linestream);
371       //RPV
372       if (blockName == "rvnmix") ifail=rvnmix.set(linestream);
373       if (blockName == "rvumix") ifail=rvumix.set(linestream);
374       if (blockName == "rvvmix") ifail=rvvmix.set(linestream);
375       if (blockName == "rvhmix") ifail=rvhmix.set(linestream);
376       if (blockName == "rvamix") ifail=rvamix.set(linestream);
377       if (blockName == "rvlmix") ifail=rvlmix.set(linestream);
378       //CPV
379       if (blockName == "cvhmix") ifail=cvhmix.set(linestream);
380       if (blockName == "imcvhmix") ifail=imcvhmix.set(linestream);
381       //CPV + FLV
382       if (blockName == "imusqmix") ifail=imusqmix.set(linestream);
383       if (blockName == "imdsqmix") ifail=imdsqmix.set(linestream);
384       if (blockName == "imselmix") ifail=imselmix.set(linestream);
385       if (blockName == "imsnumix") ifail=imsnumix.set(linestream);
386       if (blockName == "imnmix") ifail=imnmix.set(linestream);
387       if (blockName == "imumix") ifail=imumix.set(linestream);
388       if (blockName == "imvmix") ifail=imvmix.set(linestream);
389       //NMSSM
390       if (blockName == "nmhmix") ifail=nmhmix.set(linestream);
391       if (blockName == "nmamix") ifail=nmamix.set(linestream);
392       if (blockName == "nmnmix") ifail=nmnmix.set(linestream);
393       
394       //DRbar Lagrangian parameters
395       if (blockName == "gauge") ifail=gauge.set(linestream);      
396       if (blockName == "yu") ifail=yu.set(linestream);
397       if (blockName == "yd") ifail=yd.set(linestream);
398       if (blockName == "ye") ifail=ye.set(linestream);
399       if (blockName == "au") ifail=au.set(linestream);
400       if (blockName == "ad") ifail=ad.set(linestream);
401       if (blockName == "ae") ifail=ae.set(linestream);
402       if (blockName == "hmix") ifail=hmix.set(linestream);
403       if (blockName == "msoft") ifail=msoft.set(linestream);
404       //FLV
405       if (blockName == "vckm") ifail=vckm.set(linestream);
406       if (blockName == "upmns") ifail=upmns.set(linestream);
407       if (blockName == "msq2") ifail=msq2.set(linestream);
408       if (blockName == "msu2") ifail=msu2.set(linestream);
409       if (blockName == "msd2") ifail=msd2.set(linestream);
410       if (blockName == "msl2") ifail=msl2.set(linestream);
411       if (blockName == "mse2") ifail=mse2.set(linestream);
412       if (blockName == "tu") ifail=tu.set(linestream);
413       if (blockName == "td") ifail=td.set(linestream);
414       if (blockName == "te") ifail=te.set(linestream);
415       //RPV
416       if (blockName == "rvlamlle") ifail=rvlamlle.set(linestream);
417       if (blockName == "rvlamlqd") ifail=rvlamlqd.set(linestream);
418       if (blockName == "rvlamudd") ifail=rvlamudd.set(linestream);
419       if (blockName == "rvtlle") ifail=rvtlle.set(linestream);
420       if (blockName == "rvtlqd") ifail=rvtlqd.set(linestream);
421       if (blockName == "rvtudd") ifail=rvtudd.set(linestream);
422       if (blockName == "rvkappa") ifail=rvkappa.set(linestream);
423       if (blockName == "rvd") ifail=rvd.set(linestream);
424       if (blockName == "rvm2lh1") ifail=rvm2lh1.set(linestream);
425       if (blockName == "rvsnvev") ifail=rvsnvev.set(linestream);
426       //CPV
427       if (blockName == "imau") ifail=imau.set(linestream);
428       if (blockName == "imad") ifail=imad.set(linestream);
429       if (blockName == "imae") ifail=imae.set(linestream);
430       if (blockName == "imhmix") ifail=imhmix.set(linestream);
431       if (blockName == "immsoft") ifail=immsoft.set(linestream);
432       //CPV+FLV
433       if (blockName == "imvckm") ifail=imvckm.set(linestream);
434       if (blockName == "imupmns") ifail=imupmns.set(linestream);
435       if (blockName == "immsq2") ifail=immsq2.set(linestream);
436       if (blockName == "immsu2") ifail=immsu2.set(linestream);
437       if (blockName == "immsd2") ifail=immsd2.set(linestream);
438       if (blockName == "immsl2") ifail=immsl2.set(linestream);
439       if (blockName == "immse2") ifail=immse2.set(linestream);
440       if (blockName == "imtu") ifail=imtu.set(linestream);
441       if (blockName == "imtd") ifail=imtd.set(linestream);
442       if (blockName == "imte") ifail=imte.set(linestream);
443       //NMSSM
444       if (blockName == "nmssmrun") ifail=nmssmrun.set(linestream);      
445
446       //Diagnostics
447       if (ifail != 0) { 
448         if (ifail == -2) {
449           message(1,"readFile","Ignoring unknown block: "+blockName,iLine);
450           blockIn="";
451         };
452         if (ifail == -1) {
453           message(2,"readFile","Error on line.",iLine);        
454         };
455         if (ifail == 1) {
456           message(0,"readFile",blockName+" existing entry overwritten.",iLine);
457         };
458       };
459       
460     } 
461
462     // Decay table read-in
463     else if (decay != "") {
464       if (! decayPrinted) {
465         if (verbose >= 2) 
466           message(0,"readFile","reading  DECAY table for "+nameNow,0);
467         decayPrinted = true;
468       }
469       //      cout << " found decay mode "<<line<<endl;
470       double brat;
471       bool ok=true;
472       int nDa = 0;
473       vector<int> idDa;
474       istringstream linestream(line);
475       linestream >> brat;
476       if (! linestream) ok = false;
477       if (ok) linestream >> nDa;
478       if (! linestream) ok = false;
479       else {
480         for (int i=0; i<nDa; i++) {
481           int idThis;
482           linestream >> idThis;
483           if (! linestream) {
484             ok = false;
485             break;
486           }
487           idDa.push_back(idThis);
488         }
489       }
490
491       // Stop reading decay channels if not consistent.
492       if (!ok || nDa < 2) {
493         message(0,"readFile","decay channel read failed",iLine);
494          
495       // Append decay channel.
496       } else {
497         decays[decayIndices[idNow]].addChannel(brat,nDa,idDa);
498       }
499     }
500   };
501
502   //Print footer 
503   printFooter();
504
505   //Return 0 if read-in successful 
506   if ( lhefRead && !foundSlhaTag) { 
507     return 102; 
508   }
509   else return iFailFile;
510     
511 }
512
513 //--------------------------------------------------------------------------
514
515 // Print a header with information on version, last date of change, etc.
516
517 void SusyLesHouches::printHeader() {
518   if (verbose == 0) return;
519   setprecision(3);
520   if (! headerPrinted) {
521     cout <<" *--------------------  SusyLesHouches v0.05 SUSY/BSM Interface  ---------------------*\n";
522     message(0,"","Last Change 12 May 2009 - P. Z. Skands",0);
523     headerPrinted=true;
524   }
525 }
526
527 //--------------------------------------------------------------------------
528
529 // Print a footer
530
531 void SusyLesHouches::printFooter() {
532   if (verbose == 0) return;
533   if (! footerPrinted) {
534     //    cout << " *"<<endl;
535     cout <<" *------------------------------------------------------------------------------------*\n";
536     footerPrinted=true;
537     //    headerPrinted=false; 
538   }
539 }
540
541 //--------------------------------------------------------------------------
542
543 // Print the current spectrum on stdout.
544 // Not yet fully implemented.
545
546 void SusyLesHouches::printSpectrum() {
547
548   // Exit if output switched off
549   if (verbose <= 0) return;
550
551   // Print header if not already done
552   if (! headerPrinted) printHeader();
553   message(0,"","");
554
555   // Print Calculator and File name
556   if (slhaRead) {
557     message(0,"","  Spectrum Calculator was:   "+spinfo(1)+"   version: "+spinfo(2));
558     if (lhefRead) message(0,"","  Read <slha> spectrum from: "+slhaFile);
559     else message(0,"","  Read SLHA spectrum from: "+slhaFile);
560   }
561
562   // gluino
563   message(0,"","");
564   cout<<" |  ~g                  m"<<endl; 
565   cout<<setprecision(3)<<" |     1000021 "<<setw(10)<<
566       ( (mass(2000003) > 1e7) ? scientific : fixed)<<mass(1000021)<<endl;
567
568   // d squarks
569   message(0,"","");
570   cout<<" |  ~d                  m     ~dL     ~sL     ~bL     ~dR     ~sR     ~bR"<<endl;
571
572   cout<<setprecision(3) <<" |     1000001 "<<setw(10)<<
573     ( (mass(1000001) > 1e7) ? scientific : fixed)<<mass(1000001)<<fixed<<"  ";
574   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(1,icur)<<"  ";
575
576   cout<<endl<<" |     1000003 "<<setw(10)<<
577     ( (mass(1000003) > 1e7) ? scientific : fixed)<<mass(1000003)<<fixed<<"  ";
578   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(2,icur)<<"  ";
579
580   cout<<endl<<" |     1000005 "<<setw(10)<<
581     ( (mass(1000005) > 1e7) ? scientific : fixed)<<mass(1000005)<<fixed<<"  ";
582   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(3,icur)<<"  ";
583
584   cout<<endl<<" |     2000001 "<<setw(10)<<
585     ( (mass(2000001) > 1e7) ? scientific : fixed)<<mass(2000001)<<fixed<<"  ";
586   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(4,icur)<<"  ";
587
588   cout<<endl<<" |     2000003 "<<setw(10)<<
589     ( (mass(2000003) > 1e7) ? scientific : fixed)<<mass(2000003)<<fixed<<"  ";
590   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(5,icur)<<"  ";
591
592   cout<<endl<<" |     2000005 "<<setw(10)<<
593     ( (mass(2000005) > 1e7) ? scientific : fixed)<<mass(2000005)<<fixed<<"  ";
594   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<dsqmix(6,icur)<<"  ";
595
596   cout<<endl;
597   
598   // u squarks
599   message(0,"","");
600   cout<<" |  ~u                  m     ~uL     ~cL     ~tL     ~uR     ~cR     ~tR"<<endl; 
601
602   cout<<setprecision(3)<<" |     1000002 "<<setw(10)<<
603     ( (mass(1000002) > 1e7) ? scientific : fixed)<<mass(1000002)<<fixed<<"  ";
604   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(1,icur)<<"  ";
605
606   cout<<endl<<" |     1000004 "<<setw(10)<<
607     ( (mass(1000004) > 1e7) ? scientific : fixed)<<mass(1000004)<<fixed<<"  ";
608   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(2,icur)<<"  ";
609
610   cout<<endl<<" |     1000006 "<<setw(10)<<
611     ( (mass(1000006) > 1e7) ? scientific : fixed)<<mass(1000006)<<fixed<<"  "; 
612   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(3,icur)<<"  ";
613
614   cout<<endl<<" |     2000002 "<<setw(10)<<
615     ( (mass(2000002) > 1e7) ? scientific : fixed)<<mass(2000002)<<fixed<<"  "; 
616   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(4,icur)<<"  ";
617
618   cout<<endl<<" |     2000004 "<<setw(10)<<
619     ( (mass(2000004) > 1e7) ? scientific : fixed)<<mass(2000004)<<fixed<<"  " ;
620   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(5,icur)<<"  ";
621
622   cout<<endl<<" |     2000006 "<<setw(10)<<
623     ( (mass(2000006) > 1e7) ? scientific : fixed)<<mass(2000006)<<fixed<<"  ";
624   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<usqmix(6,icur)<<"  ";
625
626   cout<<endl;
627
628   // sleptons  
629   // (NB: should be an if/then/else for RPV case here)
630   message(0,"","");
631   cout<<" |  ~e                  m     ~eL    ~muL   ~tauL     ~eR    ~muR   ~tauR"<<endl; 
632
633   cout<<setprecision(3)<<" |     1000011 "<<setw(10)<<
634     ( (mass(1000011) > 1e7) ? scientific : fixed)<<mass(1000011)<<fixed<<"  ";
635   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(1,icur)<<"  ";
636
637   cout<<endl<<" |     1000013 "<<setw(10)<<
638     ( (mass(1000013) > 1e7) ? scientific : fixed)<<mass(1000013)<<fixed<<"  ";
639   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(2,icur)<<"  ";
640
641   cout<<endl<<" |     1000015 "<<setw(10)<<
642     ( (mass(1000015) > 1e7) ? scientific : fixed)<<mass(1000015)<<fixed<<"  "; 
643   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(3,icur)<<"  ";
644
645   cout<<endl<<" |     2000011 "<<setw(10)<<
646     ( (mass(2000011) > 1e7) ? scientific : fixed)<<mass(2000011)<<fixed<<"  "; 
647   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(4,icur)<<"  ";
648
649   cout<<endl<<" |     2000013 "<<setw(10)<<
650     ( (mass(2000013) > 1e7) ? scientific : fixed)<<mass(2000013)<<fixed<<"  " ;
651   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(5,icur)<<"  ";
652
653   cout<<endl<<" |     2000015 "<<setw(10)<<
654     ( (mass(2000015) > 1e7) ? scientific : fixed)<<mass(2000015)<<fixed<<"  ";
655   for (int icur=1;icur<=6;icur++) cout<<setw(6)<<selmix(6,icur)<<"  ";
656
657   cout<<endl;
658
659   // sneutrinos
660   // (NB: should be an if/then/else for RPV case here)
661   message(0,"","");
662   cout<<" |  ~nu                 m";
663   if (snumix.exists()) cout<<"   ~nu_e  ~nu_mu ~nu_tau";
664   cout<<endl; 
665
666   cout<<setprecision(3)<<" |     1000012 "<<setw(10)<<
667     ( (mass(1000012) > 1e7) ? scientific : fixed)<<mass(1000012)<<fixed<<"  ";
668   if (snumix.exists()) for (int icur=1;icur<=3;icur++) 
669     cout<<setw(6)<<snumix(1,icur)<<"  ";
670
671   cout<<endl<<" |     1000014 "<<setw(10)<<
672     ( (mass(1000014) > 1e7) ? scientific : fixed)<<mass(1000014)<<fixed<<"  ";
673   if (snumix.exists()) for (int icur=1;icur<=3;icur++) 
674     cout<<setw(6)<<snumix(2,icur)<<"  ";
675
676   cout<<endl<<" |     1000016 "<<setw(10)<<
677     ( (mass(1000016) > 1e7) ? scientific : fixed)<<mass(1000016)<<fixed<<"  "; 
678   if (snumix.exists()) for (int icur=1;icur<=3;icur++) 
679     cout<<setw(6)<<snumix(3,icur)<<"  ";
680
681   cout<<endl;
682
683   // neutralinos
684   // (NB: should be an if/then/else for RPV case here)
685   message(0,"","");
686   cout<<" |  ~chi0               m      ~B    ~W_3    ~H_1    ~H_2"<<endl; 
687
688   cout<<setprecision(3)<<" |     1000022 "<<setw(10)<<
689     ( (mass(1000022) > 1e7) ? scientific : fixed)<<mass(1000022)<<fixed<<"  ";
690   for (int icur=1;icur<=4;icur++) cout<<setw(6)<<nmix(1,icur)<<"  ";
691
692   cout<<endl<<" |     1000023 "<<setw(10)<<
693     ( (mass(1000023) > 1e7) ? scientific : fixed)<<mass(1000023)<<fixed<<"  ";
694   for (int icur=1;icur<=4;icur++) cout<<setw(6)<<nmix(2,icur)<<"  ";
695
696   cout<<endl<<" |     1000025 "<<setw(10)<<
697     ( (mass(1000025) > 1e7) ? scientific : fixed)<<mass(1000025)<<fixed<<"  "; 
698   for (int icur=1;icur<=4;icur++) cout<<setw(6)<<nmix(3,icur)<<"  ";
699
700   cout<<endl<<" |     1000035 "<<setw(10)<<
701     ( (mass(1000035) > 1e7) ? scientific : fixed)<<mass(1000035)<<fixed<<"  "; 
702   for (int icur=1;icur<=4;icur++) cout<<setw(6)<<nmix(4,icur)<<"  ";
703
704   cout<<endl;
705
706   // charginos
707   // (NB: should be an if/then/else for RPV case here)
708   message(0,"","");
709   cout<<" |  ~chi+               m   U:   ~W      ~H  |  V:   ~W      ~H"
710       <<endl; 
711
712   cout<<setprecision(3)<<" |     1000024 "<<setw(10)<<
713     ((mass(1000024) > 1e7) ? scientific : fixed)<<mass(1000024)<<fixed<<"    ";
714   for (int icur=1;icur<=2;icur++) cout<<setw(6)<<umix(1,icur)<<"  ";
715   cout<<"|   ";
716   for (int icur=1;icur<=2;icur++) cout<<setw(6)<<vmix(1,icur)<<"  ";
717
718   cout<<endl<<" |     1000037 "<<setw(10)<<
719     ((mass(1000037) > 1e7) ? scientific : fixed)<<mass(1000037)<<fixed<<"    ";
720   for (int icur=1;icur<=2;icur++) cout<<setw(6)<<umix(2,icur)<<"  ";
721   cout<<"|   " ;
722   for (int icur=1;icur<=2;icur++) cout<<setw(6)<<vmix(2,icur)<<"  ";
723
724   cout<<endl;
725
726   // Higgs bosons 
727   // (NB: should be an if/then/else for RPV case here)
728   // ...
729
730   // Print footer  
731   footerPrinted=false;
732   message(0,"","");
733   printFooter();
734 }
735
736 //--------------------------------------------------------------------------
737
738 // Check consistency of spectrum, unitarity of matrices, etc.
739
740 int SusyLesHouches::checkSpectrum() {
741
742   if (! headerPrinted) printHeader();
743   int ifail=0;
744   bool foundModsel = modsel.exists();
745   if (! foundModsel) {
746     if (mass.exists()) return 1;
747     else return 2;
748   }
749
750   // Step 1) Check MODSEL. Assign default values where applicable.
751   if (!modsel.exists(1)) {
752     message(1,"checkSpectrum","MODSEL(1) undefined. Assuming =0.",0);
753     modsel.set(1,0);
754     ifail=-1;
755   }
756   if (!modsel.exists(3)) modsel.set(3,0);
757   if (!modsel.exists(4)) modsel.set(4,0);
758   if (!modsel.exists(5)) modsel.set(5,0);
759   if (!modsel.exists(6)) modsel.set(6,0);
760   if (!modsel.exists(11)) modsel.set(11,1);
761   
762   // Step 2) Check for existence / duplication of blocks
763
764   //Global
765   if (!minpar.exists()) {
766       message(1,"checkSpectrum","MINPAR not found",0);
767       ifail=-1;    
768   }
769   if (!sminputs.exists()) {
770       message(1,"checkSpectrum","SMINPUTS not found",0);
771       ifail=-1;    
772   }
773   if (!mass.exists()) {
774       message(1,"checkSpectrum","MASS not found",0);
775       ifail=-1;  
776   }
777   if (!gauge.exists()) {
778       message(1,"checkSpectrum","GAUGE not found",0);
779       ifail=-1;    
780   }
781
782   //SLHA1
783   if (modsel(3) == 0 && modsel(4) == 0 && modsel(5) == 0 && modsel(6) == 0) {
784     // Check for required SLHA1 blocks
785     if (!staumix.exists()) {
786       message(1,"checkSpectrum","STAUMIX not found",0);
787       ifail=-1;
788     };  
789     if (!sbotmix.exists()) {
790       message(1,"checkSpectrum","SBOTMIX not found",0);
791       ifail=-1;
792     };  
793     if (!stopmix.exists()) {
794       message(1,"checkSpectrum","STOPMIX not found",0);
795       ifail=-1;
796     };  
797     if (!nmix.exists()) {
798       message(1,"checkSpectrum","NMIX not found",0);
799       ifail=-1;
800     };  
801     if (!umix.exists()) {
802       message(1,"checkSpectrum","UMIX not found",0);
803       ifail=-1;
804     };  
805     if (!vmix.exists()) {
806       message(1,"checkSpectrum","VMIX not found",0);
807       ifail=-1;
808     };  
809     if (!alpha.exists()) {
810       message(1,"checkSpectrum","ALPHA not found",0);
811       ifail=-1;    
812     }
813     if (!hmix.exists()) {
814       message(1,"checkSpectrum","HMIX not found",0);
815       ifail=-1;    
816     }
817     if (!msoft.exists()) {
818       message(1,"checkSpectrum","MSOFT not found",0);
819       ifail=-1;    
820     }
821   } 
822
823   //RPV (+ FLV)
824   else if (modsel(4) != 0) {
825     // Check for required SLHA2 blocks
826     if (!rvnmix.exists()) {
827       message(1,"checkSpectrum","MODSEL 4 != 0 but RVNMIX not found",0);
828       ifail=-1;
829     }
830     if (!rvumix.exists()) {
831       message(1,"checkSpectrum","MODSEL 4 != 0 but RVUMIX not found",0);
832       ifail=-1;
833     }
834     if (!rvvmix.exists()) {
835       message(1,"checkSpectrum","MODSEL 4 != 0 but RVVMIX not found",0);
836       ifail=-1;
837     }
838     if (!rvhmix.exists()) {
839       message(1,"checkSpectrum","MODSEL 4 != 0 but RVHMIX not found",0);
840       ifail=-1;
841     }
842     if (!rvamix.exists()) {
843       message(1,"checkSpectrum","MODSEL 4 != 0 but RVAMIX not found",0);
844       ifail=-1;
845     }
846     if (!rvlmix.exists()) {
847       message(1,"checkSpectrum","MODSEL 4 != 0 but RVLMIX not found",0);
848       ifail=-1;
849     }
850     if (!usqmix.exists()) {
851       message(1,"checkSpectrum","MODSEL 4 != 0 but USQMIX not found",0);
852       ifail=-1;
853     }
854     if (!dsqmix.exists()) {
855       message(1,"checkSpectrum","MODSEL 4 != 0 but DSQMIX not found",0);
856       ifail=-1;
857     }
858   }
859
860   // FLV but not RPV (see above for FLV+RPV, below for FLV regardless of RPV)
861   else if (modsel(6) != 0) {
862     // Quark FLV
863     if (modsel(6) != 2) {
864       if (!usqmix.exists()) {
865         message(1,"checkSpectrum","quark FLV on but USQMIX not found",0);
866         ifail=-1;
867       }
868       if (!dsqmix.exists()) {
869         message(1,"checkSpectrum","quark FLV on but DSQMIX not found",0);
870         ifail=-1;
871       }
872     }
873     // Lepton FLV
874     if (modsel(6) != 1) {
875       if (!upmns.exists()) {
876         message(1,"checkSpectrum","lepton FLV on but UPMNSIN not found",0);
877         ifail=-1;
878       }
879       if (!selmix.exists()) {
880         message(1,"checkSpectrum","lepton FLV on but SELMIX not found",0);
881         ifail=-1;
882       }
883       if (!snumix.exists() && !snsmix.exists()) {
884         message(1,"checkSpectrum","lepton FLV on but SNUMIX not found",0);
885         ifail=-1;
886       }
887     }
888   }
889
890   // CPV
891   if (modsel(5) != 0) {
892     if (!cvhmix.exists()) {
893       message(1,"checkSpectrum","MODSEL 5 != 0 but CVHMIX not found",0);
894       ifail=-1;
895     }
896   }
897
898   // FLV (regardless of whether RPV or not)
899   if (modsel(6) != 0) {
900     // Quark FLV
901     if (modsel(6) != 2) {
902       if (!vckmin.exists()) {
903         message(1,"checkSpectrum","quark FLV on but VCKMIN not found",0);
904         ifail=-1;
905       }
906       if (!msq2in.exists()) {
907         message(0,"checkSpectrum","note: quark FLV on but MSQ2IN not found",0);
908         ifail=min(ifail,0);
909       }
910       if (!msu2in.exists()) {
911         message(0,"checkSpectrum","note: quark FLV on but MSU2IN not found",0);
912         ifail=min(ifail,0);
913       }
914       if (!msd2in.exists()) {
915         message(0,"checkSpectrum","note: quark FLV on but MSD2IN not found",0);
916         ifail=min(ifail,0);
917       }
918       if (!tuin.exists()) {
919         message(0,"checkSpectrum","note: quark FLV on but TUIN not found",0);
920         ifail=min(ifail,0);
921       }
922       if (!tdin.exists()) {
923         message(0,"checkSpectrum","note: quark FLV on but TDIN not found",0);
924         ifail=min(ifail,0);
925       }
926     }
927     // Lepton FLV
928     if (modsel(6) != 1) {
929       if (!msl2in.exists()) {
930         message(0,"checkSpectrum","note: lepton FLV on but MSL2IN not found",0);
931         ifail=min(ifail,0);
932       }
933       if (!mse2in.exists()) {
934         message(0,"checkSpectrum","note: lepton FLV on but MSE2IN not found",0);
935         ifail=min(ifail,0);
936       }
937       if (!tein.exists()) {
938         message(0,"checkSpectrum","note: lepton FLV on but TEIN not found",0);
939         ifail=min(ifail,0);
940       }
941     }
942   }
943   
944   // Step 3) SLHA1 --> SLHA2 interoperability
945   //Note: the mass basis is NOT mass-ordered in SLHA1, so be careful!
946   //Here, the mass basis is hence by PDG code, not by mass-ordered value.
947
948   if (stopmix.exists() && ! usqmix.exists() ) {
949     //1000002 = ~uL, 1000004 = ~cL, 2000002 = ~uR, 2000004 = ~cR 
950     usqmix.set(1,1, 1.0);
951     usqmix.set(2,2, 1.0); 
952     usqmix.set(4,4, 1.0);
953     usqmix.set(5,5, 1.0);
954     //Fill (1000006,2000006) sector from stopmix
955     usqmix.set(3,3, stopmix(1,1));
956     usqmix.set(3,6, stopmix(1,2));
957     usqmix.set(6,3, stopmix(2,1));
958     usqmix.set(6,6, stopmix(2,2));    
959   };
960   if (sbotmix.exists() && ! dsqmix.exists() ) {
961     //1000001 = ~dL, 1000003 = ~sL, 2000001 = ~dR, 2000003 = ~sR 
962     dsqmix.set(1,1, 1.0);
963     dsqmix.set(2,2, 1.0); 
964     dsqmix.set(4,4, 1.0);
965     dsqmix.set(5,5, 1.0);
966     //Fill (1000005,2000005) sector from sbotmix
967     dsqmix.set(3,3, sbotmix(1,1));
968     dsqmix.set(3,6, sbotmix(1,2));
969     dsqmix.set(6,3, sbotmix(2,1));
970     dsqmix.set(6,6, sbotmix(2,2));
971   };
972   if (staumix.exists() && ! selmix.exists() ) {
973     //1000011 = ~eL, 1000013 = ~muL, 2000011 = ~eR, 2000013 = ~muR 
974     selmix.set(1,1, 1.0);
975     selmix.set(2,2, 1.0); 
976     selmix.set(4,4, 1.0);
977     selmix.set(5,5, 1.0);
978     //Fill (1000015,2000015) sector from staumix
979     selmix.set(3,3, staumix(1,1));
980     selmix.set(3,6, staumix(1,2));
981     selmix.set(6,3, staumix(2,1));
982     selmix.set(6,6, staumix(2,2));
983   };
984   if (! snumix.exists() && ! snsmix.exists()) {
985     //1000012 = ~nu_e, 1000014 = ~nu_mu, 1000016 = ~nu_tau
986     snumix.set(1,1, 1.0);
987     snumix.set(2,2, 1.0); 
988     snumix.set(3,3, 1.0);
989   };
990
991   // Step 4) Check unitarity/orthogonality of mixing matrices
992
993   //NMIX
994   if (nmix.exists()) {
995     for (int i=1;i<=4;i++) {
996       double cn1=0.0;
997       double cn2=0.0;
998       for (int j=1;j<=4;j++) {
999         cn1 += pow(nmix(i,j),2);
1000         cn2 += pow(nmix(j,i),2);
1001       }
1002       if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) { 
1003         ifail=2; 
1004         message(2,"checkSpectrum","inconsistent normalization of NMIX",0);
1005       }
1006     }
1007   }
1008
1009   //VMIX, UMIX
1010   if (vmix.exists() && umix.exists()) {
1011     for (int i=1;i<=2;i++) {
1012       double cu1=0.0;
1013       double cu2=0.0;
1014       double cv1=0.0;
1015       double cv2=0.0;
1016       for (int j=1;j<=2;j++) {
1017         cu1 += pow(umix(i,j),2);
1018         cu2 += pow(umix(j,i),2);
1019         cv1 += pow(vmix(i,j),2);
1020         cv2 += pow(vmix(j,i),2);
1021       }
1022       if (abs(1.0-cu1) > 1e-3 || abs(1.0-cu2) > 1e-3) { 
1023         ifail=2; 
1024         message(2,"checkSpectrum","inconsistent normalization of UMIX",0);
1025       }
1026       if (abs(1.0-cv1) > 1e-3 || abs(1.0-cv2) > 1e-3) { 
1027         ifail=2; 
1028         message(2,"checkSpectrum","inconsistent normalization of VMIX",0);
1029       }
1030     }
1031     
1032   }
1033
1034   //STOPMIX, SBOTMIX
1035   if (stopmix.exists() && sbotmix.exists()) {
1036     for (int i=1;i<=2;i++) {
1037       double ct1=0.0;
1038       double ct2=0.0;
1039       double cb1=0.0;
1040       double cb2=0.0;
1041       for (int j=1;j<=2;j++) {
1042         ct1 += pow(stopmix(i,j),2);
1043         ct2 += pow(stopmix(j,i),2);
1044         cb1 += pow(sbotmix(i,j),2);
1045         cb2 += pow(sbotmix(j,i),2);
1046       }
1047       if (abs(1.0-ct1) > 1e-3 || abs(1.0-ct2) > 1e-3) { 
1048         ifail=2; 
1049         message(2,"checkSpectrum","inconsistent normalization of STOPMIX",0);
1050       }
1051       if (abs(1.0-cb1) > 1e-3 || abs(1.0-cb2) > 1e-3) { 
1052         ifail=2; 
1053         message(2,"checkSpectrum","inconsistent normalization of SBOTMIX",0);
1054       }
1055     }    
1056   }
1057
1058   //STAUMIX
1059   if (staumix.exists()) {
1060     for (int i=1;i<=2;i++) {
1061       double ct1=0.0;
1062       double ct2=0.0;
1063       for (int j=1;j<=2;j++) {
1064         ct1 += pow(staumix(i,j),2);
1065         ct2 += pow(staumix(j,i),2);
1066       }
1067       if (abs(1.0-ct1) > 1e-3 || abs(1.0-ct2) > 1e-3) { 
1068         ifail=2; 
1069         message(2,"checkSpectrum","inconsistent normalization of STAUMIX",0);
1070       }
1071     }    
1072   }
1073
1074   //NMSSM:
1075   if (modsel(3) == 1) {
1076     //NMNMIX
1077     if ( nmnmix.exists() ) {
1078       for (int i=1;i<=5;i++) {
1079         double cn1=0.0;
1080         double cn2=0.0;
1081         for (int j=1;j<=4;j++) {
1082           cn1 += pow(nmnmix(i,j),2);
1083           cn2 += pow(nmnmix(j,i),2);
1084         }
1085         if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) { 
1086           ifail=max(ifail,2); 
1087           message(2,"checkSpectrum","inconsistent normalization of NMNMIX",0);
1088         }
1089       }
1090     }
1091     else {
1092       ifail=-1;
1093       message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMNMIX found",0);
1094     }
1095     //NMAMIX
1096     if ( nmamix.exists() ) {
1097       for (int i=1;i<=2;i++) {
1098         double cn1=0.0;
1099         for (int j=1;j<=3;j++) {
1100           cn1 += pow(nmamix(i,j),2);
1101         }
1102         if (abs(1.0-cn1) > 1e-3) { 
1103           ifail=max(ifail,2); 
1104           message(2,"checkSpectrum","inconsistent normalization of NMAMIX",0);
1105         }
1106       }
1107     }
1108     else {
1109       ifail=-1;
1110       message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMAMIX found",0);
1111     }
1112     //NMHMIX
1113     if ( nmhmix.exists() ) {
1114       for (int i=1;i<=3;i++) {
1115         double cn1=0.0;
1116         double cn2=0.0;
1117         for (int j=1;j<=3;j++) {
1118           cn1 += pow(nmhmix(i,j),2);
1119           cn2 += pow(nmhmix(j,i),2);
1120         }
1121         if (abs(1.0-cn1) > 1e-3 || abs(1.0-cn2) > 1e-3) { 
1122           ifail=max(ifail,2); 
1123           message(2,"checkSpectrum","inconsistent normalization of NMHMIX",0);
1124         }
1125       }
1126     }
1127     else {
1128       ifail=-1;
1129       message(1,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMHMIX found",0);
1130     }
1131     //NMSSMRUN
1132     if (! nmssmrun.exists() ) {
1133       ifail=-1;
1134       message(2,"checkSpectrum","MODSEL 3 = 1 (NMSSM) but no NMSSMRUN found",
1135               0);
1136     }
1137   }
1138   
1139   //Check for documentation
1140   if (slhaRead && ! spinfo.exists(1)) spinfo.set(1,"unknown");
1141   if (slhaRead && ! spinfo.exists(2)) spinfo.set(2,"unknown");
1142   if (! slhaRead && ! spinfo.exists(1)) {
1143     spinfo.set(1,"DEFAULT");
1144     spinfo.set(2,"n/a");
1145   }
1146
1147   //Give status 
1148   if (ifail >= 2) 
1149     message(0,"checkSpectrum","one or more serious problems were found");
1150
1151   //Print Footer
1152   printFooter();
1153
1154   //Return
1155   return ifail;
1156 }
1157
1158 //--------------------------------------------------------------------------
1159
1160 // Check consistency of decay tables
1161
1162 int SusyLesHouches::checkDecays() {
1163
1164   if (! headerPrinted) printHeader();
1165   int iFailDecays=0;
1166   
1167   // Loop over all particles read in
1168   for (int i = 0; i < int(decays.size()); ++i) { 
1169     
1170     // Shorthand
1171     decayTable decTab = decays[i];
1172     int idRes = decTab.getId();
1173     double width = decTab.getWidth();
1174     if (width <= 0.0 || decTab.size() == 0) continue;
1175     
1176     // Check sum of branching ratios and phase spaces
1177     double sum = 0.0;
1178     double absSum = 0.0;
1179     int decSize = decTab.size();
1180     for (int j = 0; j < decSize; ++j) {
1181       
1182       double brat = decTab.getBrat(j);
1183       
1184       // Check phase space
1185       if (abs(brat) > 0.0) {
1186         vector<int> idDa = decTab.getIdDa(j);
1187         double massSum=abs(mass(idRes));      
1188         for (int k=0; k<int(idDa.size()); ++k) {
1189           if (mass.exists(idDa[k])) massSum -= mass(abs(idDa[k]));
1190           // If no MASS information read, use lowish values for check
1191           else if (abs(idDa[k]) == 24) massSum -=  79.0;
1192           else if (abs(idDa[k]) == 23) massSum -=  91.0;
1193           else if (abs(idDa[k]) ==  6) massSum -= 165.0;
1194           else if (abs(idDa[k]) ==  5) massSum -=   4.0;
1195           else if (abs(idDa[k]) ==  4) massSum -=   1.0;
1196         }
1197         if (massSum < 0.0) {
1198           // String containing decay name
1199           ostringstream errCode;
1200           errCode << idRes <<" ->";
1201           for (int jDa=0; jDa<int(idDa.size()); ++jDa) errCode<<" "<<idDa[jDa];
1202           message(1,"checkDecays",errCode.str()
1203                   +": Phase Space Closed, but BR != 0");
1204           iFailDecays = 1;
1205         }
1206         
1207       }
1208
1209       // Sum up branching rations
1210       sum += brat;
1211       absSum += abs(brat);
1212       
1213     }
1214
1215     if (abs(1.0-absSum) > 1e-6) {
1216       message(1,"checkDecays","sum(BR) != 1");
1217       cout << " | offending particle: "<<idRes<<" sum(BR) = "<<absSum<<endl;      
1218       iFailDecays = 2;
1219     }
1220     
1221   }
1222   // End of loop over particles. Done.
1223   
1224   return iFailDecays;
1225
1226 }
1227
1228 //--------------------------------------------------------------------------
1229
1230 // Simple utility to print messages, warnings, and errors
1231
1232 void SusyLesHouches::message(int level, string place,string themessage,int line) {
1233   if (verbose == 0) return;
1234   //Send normal messages and warnings to stdout, errors to stderr.
1235   ostream* outstream = &cerr;
1236   if (level <= 1) outstream = &cout;
1237   // if (level == 2) { *outstream<<endl; }
1238   if (place != "") *outstream << " | (SLHA::"+place+") ";
1239   else *outstream << " | ";
1240   if (level == 1) *outstream<< "warning: "; 
1241   if (level == 2) { *outstream <<"ERROR: "; } 
1242   if (line != 0) *outstream<< "line "<<line<<" - ";
1243   *outstream << themessage << endl;
1244   //  if (level == 2) *outstream <<endl;
1245   footerPrinted=false;
1246   return;
1247 }
1248
1249 //==========================================================================
1250
1251
1252
1253