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