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