]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8170/src/SusyLesHouches.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / SusyLesHouches.cxx
CommitLineData
63ba5337 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