- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / include / SusyLesHouches.h
1 // SusyLesHouches.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 Peter Skands, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6 // Header file for SUSY Les Houches Accord Interface
7 // Is independent of the rest of the Pythia implementation and thus could
8 // be re-used stand-alone or merged into other applications, subject to 
9 // the MCnet Guidelines mentioned above.
10
11 #ifndef SLHA_H
12 #define SLHA_H
13
14 // Stdlib header files for string and character manipulation.
15 #include <string>
16 #include <cctype>
17 // Stdlib header files for containers.
18 #include <vector>
19 #include <map>
20 // Stdlib header files for input/output.
21 #include <iostream>
22 #include <iomanip>
23 #include <fstream>
24 #include <sstream>
25 // Stdlib header files for mathematics.
26 #include <cmath>
27 #include <cstdlib>
28
29 // Stdlib namespace
30 using namespace std;
31
32 class SusyLesHouches {
33
34 public:
35
36   //Constructor, with and without filename.
37   SusyLesHouches(int verboseIn=1) : verbose(verboseIn), 
38     headerPrinted(false), footerPrinted(false),
39     slhaRead(false), lhefRead(false), lhefSlha(false) {};
40   SusyLesHouches(string filename, int verboseIn=1) : verbose(verboseIn), 
41     headerPrinted(false), footerPrinted(false),
42     slhaRead(true), lhefRead(false), lhefSlha(false) {readFile(filename);};
43
44   //***************************** SLHA FILE I/O *****************************//
45   // Read and write SLHA files 
46   int readFile(string slhaFileIn="slha.spc",int verboseIn=1); 
47   //int writeFile(string filename): write SLHA file on filename
48
49   //Output utilities
50   void printHeader();   // print Header
51   void printFooter();   // print Footer
52   void printSpectrum(); // print Spectrum
53   
54   // Check spectrum and decays
55   int checkSpectrum();
56   int checkDecays();
57
58   // File Name (can be either SLHA or LHEF)
59   string slhaFile;
60
61   // Class for SLHA data entry
62   class Entry {
63     
64   public:
65     //Constructor. 
66     Entry() : isIntP(false), isDoubleP(false), 
67       isStringP(false), n(0), d(0.0), s(""), commentP("") {}
68     
69     // Generic functions to inquire whether an int, double, or string
70     bool isInt(){return isIntP;}
71     bool isDouble(){return isDoubleP;}
72     bool isString(){return isStringP;}
73
74     // = Overloading: Set entry to int, double, or string
75     Entry& operator=(double& val)  {
76       d=val;isIntP=false;isDoubleP=true;isStringP=false;
77       return *this;      
78     };
79     Entry& operator=(int& val)  {
80       n=val;isIntP=true;isDoubleP=false;isStringP=false;
81       return *this;
82     };
83     Entry& operator=(string& val)  {
84       s=val;isIntP=false;isDoubleP=false;isStringP=true;
85       return *this;
86     };
87     
88     // Set and Get comment
89     void setComment(string comment) {commentP=comment;}
90     void getComment(string comment) {comment=commentP;}
91
92     // Generic functions to get value
93     bool get(int& val) {val=n; return isIntP;}
94     bool get(double& val) {val=d; return isDoubleP;}
95     bool get(string& val) {val=s; return isStringP;}
96
97   private:
98     bool isIntP, isDoubleP, isStringP;    
99     int n;
100     double d;
101     string s;
102     string commentP;
103     
104   };
105
106   //***************************** SLHA CLASSES *****************************//
107
108
109   //class block: the generic SLHA block (see below for matrices)
110   //Explicit typing required, e.g. block<double> minpar;
111   template <class T> class block {    
112
113   public: 
114
115     //Constructor. 
116     block<T>() : idnow(0) { } ;    
117
118     //Does block exist?
119     bool exists() { return int(entry.size()) == 0 ? false : true ; };
120     //Clear block
121     void clear() { entry.clear(); };
122
123     //set: set block entry values.
124     //Possible return values from set:
125     // 0: normal return. Entry did not previously exist and has been created.
126     // 1: normal return. Entry did previously exist and has been overwritten.
127     //-1: failure. 
128     int set(int iIn,T valIn) { 
129       int alreadyexisting=exists(iIn)?1:0;
130       entry[iIn]=valIn; 
131       return alreadyexisting;
132     };
133     // Read index and value from SLHA data line
134     int set(istringstream& linestream) {
135       linestream >> i >> val;
136       return linestream ? set(i,val) : -1;
137     };
138     // With i already given, read value from remaining SLHA data line
139     int set(int iIn,istringstream& linestream) {
140       linestream >> val;
141       return linestream ? set(iIn,val) : -1;
142     };
143     // Shorthand for entry[0]. Used e.g. for block ALPHA.
144     void set(T valIn) { entry[0]=valIn; };
145
146     // Does entry i already exist in this block?
147     bool exists(int iIn) {return entry.find(iIn) != entry.end() 
148       ? true : false;};
149
150     // Indexing with (). Output only.
151     T operator()(int iIn=0) {
152       if (exists(iIn)) {return entry[iIn];} else {T dummy(0); return dummy;};
153     };
154
155     // Size of map
156     int size() {return int(entry.size());};
157
158     // First and next key code
159     int first() { idnow = entry.begin()->first; return idnow; };
160     int next() { 
161       typename map<int,T>::iterator itnow;
162       itnow = ++entry.find(idnow);
163       if ( itnow == entry.end() ) itnow=entry.begin();
164       return idnow = itnow->first;
165     };
166
167     // Simple print utility
168     void print() {
169       bool finished=false;
170       int ibegin=first();
171       i=ibegin;
172       while (!finished) {
173         cout << "  "<< i << " " << entry[i] <<endl;
174         i=next();
175         if (i == ibegin) finished=true;
176       };       
177     };
178
179     // Special for DRbar running blocks.
180     void setq(double qIn) { qDRbar=qIn; }
181     double q() { return qDRbar; }
182  
183   private:
184     map<int,T> entry;    
185     int idnow;
186     double qDRbar;
187     //Auxiliary vars
188     int i; 
189     T val;
190   };
191
192   // class matrixblock: the generic SLHA matrix 
193   // Explicit sizing required, e.g. matrixblock<4> nmix;
194   template <int size> class matrixblock {    
195   public: 
196     //Constructor. Set uninitialized and explicitly zero.
197     matrixblock<size>() { 
198       initialized=false; 
199       for (i=1;i<=size;i++) {
200         for (j=1;j<=size;j++) {
201           entry[i][j]=0.0;
202         };
203       };
204     };    
205
206     // Assignment
207     matrixblock& operator=(const matrixblock& m) { 
208       if (this != &m) { 
209         for (i=0;i<size;i++) for (j=0;j<=size;j++) entry[i][j] = m(i,j);
210         qDRbar = m.qDRbar; 
211         initialized = m.initialized; 
212       } 
213       return *this; };
214
215     // Does this matrix contain any entries?
216     bool exists() { return initialized; };
217     // Clear initialized flag
218     void clear() { initialized=false; };
219
220     // Set matrix entry
221     int set(int iIn,int jIn, double valIn) { 
222       if (iIn>0 && jIn>0 && iIn<=size && jIn<=size) {
223         entry[iIn][jIn]=valIn;
224         initialized=true;
225         return 0;
226       } else {
227         return -1;
228       };
229     };
230
231     // Set entry from linestream (used during file read)
232     int set(istringstream& linestream) {
233       linestream >> i >> j >> val;
234       return linestream ? set(i,j,val) : -1;
235     };
236
237     // () Overloading: Get entry
238     double operator()(int iIn, int jIn) const {
239       return (iIn <= size && jIn <= size && iIn > 0 && jIn > 0) ? 
240         entry[iIn][jIn] : 0.0;
241     };
242
243     // Set and get scale for DRbar running blocks.
244     void setq(double qIn) { qDRbar=qIn; }
245     double q() { return qDRbar; }
246
247     // Simple print utility, to be elaborated on.
248     void print() {
249       for (i=1;i<=size;i++) {
250         cout << "   "<<i << " " ;
251         for (j=1;j<=size;j++) cout << entry[i][j] << " ";
252         cout << endl;
253       };
254     };
255
256   private:
257     bool initialized;
258     double entry[size+1][size+1];
259     double qDRbar;
260     //Auxiliary vars
261     int i,j; 
262     double val;
263   };
264
265   // class tensorblock: the generic SLHA tensor
266   // Explicit sizing required, e.g. tensorblock<3> rvlam;
267   template <int size> class tensor3block {    
268   public: 
269     //Constructor. Set uninitialized and explicitly zero.
270     tensor3block<size>() { 
271       initialized=false; 
272       for (i=1;i<=size;i++) {
273         for (j=1;j<=size;j++) {
274           for (k=1;k<=size;k++) {
275             entry[i][j][k]=0.0;
276           };
277         };
278       };
279     };    
280     
281     // Assignment
282     tensor3block& operator=(const tensor3block& m) { 
283       if (this != &m) { 
284         for (i=0;i<size;i++) for (j=0;j<=size;j++) for (k=0;k<=size;k++) 
285           entry[i][j][k] = m(i,j,k);
286         qDRbar = m.qDRbar; 
287         initialized = m.initialized; 
288       } 
289       return *this; };
290     
291     // Does this matrix contain any entries?
292     bool exists() { return initialized; };
293     // Clear initialized flag
294     void clear() { initialized=false; };
295     
296     // Set matrix entry
297     int set(int iIn,int jIn, int kIn, double valIn) { 
298       if (iIn>0 && jIn>0 && kIn>0 && iIn<=size && jIn<=size && kIn<=size) {
299         entry[iIn][jIn][kIn]=valIn;
300         initialized=true;
301         return 0;
302       } else {
303         return -1;
304       };
305     };
306
307     // Set entry from linestream (used during file read)
308     int set(istringstream& linestream) {
309       linestream >> i >> j >> k >> val;
310       return linestream ? set(i,j,k,val) : -1;
311     };
312
313     // () Overloading: Get entry
314     double operator()(int iIn, int jIn, int kIn) const {
315       return (iIn <= size && jIn <= size && kIn <= size && iIn > 0 
316         && jIn > 0 && kIn > 0) ? entry[iIn][jIn][kIn] : 0.0;
317     };
318
319     // Set and get scale for DRbar running blocks.
320     void setq(double qIn) { qDRbar=qIn; }
321     double q() { return qDRbar; }
322
323     // Simple print utility, to be elaborated on.
324     void print() {
325       for (i=1;i<=size;i++) {   
326         for (j=1;j<=size;j++) {
327           cout << "   "<<i << " "<<j << " " ;
328           for (k=1;k<=size;k++) {
329             cout << entry[i][j][k] << " ";         
330             cout << endl; 
331           };
332         };
333       };
334     };
335
336   private:
337     bool initialized;
338     double entry[size+1][size+1][size+1];
339     double qDRbar;
340     //Auxiliary vars
341     int i,j,k; 
342     double val;
343   };
344
345   //*************************** DECAY TABLES ***************************//
346
347   class decayChannel {
348   public: 
349
350     decayChannel() : brat(0.0) {};
351     decayChannel(double bratIn, int nDaIn, vector<int> idDaIn, string cIn="") {
352       setChannel(bratIn,nDaIn,idDaIn,cIn);
353     }
354
355     // Functions to set decay channel information
356     void setChannel(double bratIn, int nDaIn, vector<int> idDaIn, string cIn="") {
357       brat    = bratIn;
358       for (int i=0; i<=nDaIn; i++) {
359         if (i < int(idDaIn.size())) idDa.push_back(idDaIn[i]);
360         comment = cIn;
361       }
362     }
363     void setBrat(double bratIn) {brat=bratIn;}
364     void setIdDa(vector<int> idDaIn) {idDa = idDaIn;}
365     
366     // Functions to get decay channel information
367     double getBrat() {return brat;}
368     int getNDa() {return int(idDa.size());}
369     vector<int> getIdDa() {return idDa;}
370     string getComment() {return comment;}
371     
372   private:
373     double brat;
374     vector<int> idDa;
375     string comment;
376       
377   };
378
379   class decayTable {        
380   public: 
381     
382   decayTable() : id(0), width(0.0) {};
383   decayTable(int idIn) : id(idIn), width(0.0) {};
384   decayTable(int idIn, double widthIn) : id(idIn), width(widthIn) {};
385     
386     // Functions to get PDG code (id) and width
387     int    getId() {return id;}
388     double getWidth() {return width;} 
389     
390     // Functions to set PDG code (id) and width
391     void setId(int idIn) {id = idIn;}
392     void setWidth(double widthIn) {width=widthIn;}
393
394     // Function to reset size and width (width -> 0 by default)
395     void reset(double widthIn=0.0) {table.resize(0); width=widthIn;}
396     
397     // Function to add another decay channel
398     void addChannel(decayChannel channelIn) {table.push_back(channelIn);}
399     void addChannel(double bratIn, int nDaIn, vector<int> idDaIn, string cIn="") {
400       decayChannel newChannel(bratIn, nDaIn, idDaIn, cIn);
401       table.push_back(newChannel);
402     }
403
404     // Function to return number of decay channels
405     int size() {return int(table.size());}
406
407     // Function to return a branching ratio
408     double getBrat(int iChannel) {
409       if (iChannel >= 0 && iChannel < int(table.size())) {
410         return table[iChannel].getBrat();
411       } else {
412         return 0.0;
413       }
414     }
415     // Function to return daughter PDG codes 
416     vector<int> getIdDa(int iChannel) {
417       if (iChannel >= 0 && iChannel < int(table.size())) {
418         return table[iChannel].getIdDa();
419       } else {
420         vector<int> dum;
421         return dum;
422       }
423     }
424     // Function to return a decay channel
425     decayChannel getChannel(int iChannel) {
426       if (iChannel >= 0 && iChannel < int(table.size())) {
427         return table[iChannel];
428       } else {
429         decayChannel dum;
430         return dum;
431       }
432     }
433     
434   private:
435     int id;
436     double width;
437     vector<decayChannel> table;
438
439   };
440
441   //*************************** THE SLHA1 BLOCKS ***************************//
442   //blocks for model definition:
443   block<int> modsel;
444   block<int> modsel21;
445   block<double> modsel12;
446   block<double> minpar;
447   block<double> extpar;
448   block<double> sminputs;
449   //blocks for RGE program specific output
450   block<string> spinfo;
451   block<string> spinfo3;
452   block<string> spinfo4;
453   //blocks for DCY program specific output
454   block<string> dcinfo;
455   block<string> dcinfo3;
456   block<string> dcinfo4;
457   //blocks for mass and coupling spectrum
458   block<double> mass;
459   matrixblock<4> nmix;
460   matrixblock<2> umix;
461   matrixblock<2> vmix;
462   matrixblock<2> stopmix;
463   matrixblock<2> sbotmix;
464   matrixblock<2> staumix;
465   block<double> alpha;
466   block<double> hmix;
467   block<double> gauge;
468   block<double> msoft;
469   matrixblock<3> au;
470   matrixblock<3> ad;
471   matrixblock<3> ae;
472   matrixblock<3> yu;
473   matrixblock<3> yd;
474   matrixblock<3> ye;
475
476   //************************ THE SLHA1 DECAY TABLES ************************//
477   vector<decayTable> decays;
478   map<int,int> decayIndices;
479
480   //*************************** THE SLHA2 BLOCKS ***************************//
481   //Additions to SLHA1
482   block<double> qextpar;  
483
484   //FLV Input
485   block<double> vckmin;  // The input CKM Wolfenstein parms.
486   block<double> upmnsin; // The input PMNS PDG parms.
487   matrixblock<3> msq2in; // The input upper off-diagonal msq2
488   matrixblock<3> msu2in; // The input upper off-diagonal msu2
489   matrixblock<3> msd2in; // The input upper off-diagonal msd2
490   matrixblock<3> msl2in; // The input upper off-diagonal msl2
491   matrixblock<3> mse2in; // The input upper off-diagonal mse2
492   matrixblock<3> tuin;   // The input upper off-diagonal TU
493   matrixblock<3> tdin;   // The input upper off-diagonal TD
494   matrixblock<3> tein;   // The input upper off-diagonal TE
495   //FLV Output
496   matrixblock<3> vckm;    // The output DRbar running Re{VCKM} at Q
497   matrixblock<3> upmns;   // The output DRbar running Re{UPMNS} at Q
498   matrixblock<3> msq2;    // The output DRbar running msq2 at Q
499   matrixblock<3> msu2;    // The output DRbar running msu2 at Q
500   matrixblock<3> msd2;    // The output DRbar running msd2 at Q
501   matrixblock<3> msl2;    // The output DRbar running msl2 at Q
502   matrixblock<3> mse2;    // The output DRbar running mse2 at Q
503   matrixblock<3> tu;      // The output DRbar running TU at Q
504   matrixblock<3> td;      // The output DRbar running TD at Q
505   matrixblock<3> te;      // The output DRbar running TE at Q
506   matrixblock<6> usqmix;  // The Re{} up squark mixing matrix
507   matrixblock<6> dsqmix;   // The Re{} down squark mixing matrix
508   matrixblock<6> selmix;   // The Re{} selectron mixing matrix
509   matrixblock<3> snumix;   // The Re{} sneutrino mixing matrix
510   matrixblock<3> snsmix;   // The scalar sneutrino mixing matrix
511   matrixblock<3> snamix;   // The pseudoscalar neutrino mixing matrix
512
513   //RPV Input
514   tensor3block<3> rvlamllein; // The input LNV lambda couplings
515   tensor3block<3> rvlamlqdin; // The input LNV lambda' couplings
516   tensor3block<3> rvlamuddin; // The input BNV lambda'' couplings
517   tensor3block<3> rvtllein;   // The input LNV T couplings
518   tensor3block<3> rvtlqdin;   // The input LNV T' couplings
519   tensor3block<3> rvtuddin;   // The input BNV T'' couplings
520   block<double> rvkappain;    // The input LNV kappa couplings
521   block<double> rvdin;        // The input LNV D terms
522   block<double> rvm2lh1in;    // The input LNV m2LH1 couplings
523   block<double> rvsnvevin;    // The input LNV sneutrino vevs
524   //RPV Output
525   tensor3block<3> rvlamlle;   // The output LNV lambda couplings
526   tensor3block<3> rvlamlqd;   // The output LNV lambda' couplings
527   tensor3block<3> rvlamudd;   // The output BNV lambda'' couplings
528   tensor3block<3> rvtlle;     // The output LNV T couplings
529   tensor3block<3> rvtlqd;     // The output LNV T' couplings
530   tensor3block<3> rvtudd;     // The output BNV T'' couplings
531   block<double> rvkappa;      // The output LNV kappa couplings
532   block<double> rvd;          // The output LNV D terms
533   block<double> rvm2lh1;      // The output LNV m2LH1 couplings
534   block<double> rvsnvev;      // The output LNV sneutrino vevs
535   matrixblock<7> rvnmix;      // The RPV neutralino mixing matrix
536   matrixblock<5> rvumix;      // The RPV chargino L mixing matrix
537   matrixblock<5> rvvmix;      // The RPV chargino R mixing matrix
538   matrixblock<5> rvhmix;      // The RPV neutral scalar mixing matrix
539   matrixblock<5> rvamix;      // The RPV neutral pseudoscalar mixing matrix
540   matrixblock<7> rvlmix;      // The RPV charged fermion mixing matrix
541
542   //CPV Input
543   block<double> imminpar;
544   block<double> imextpar;
545   //CPV Output
546   matrixblock<4> cvhmix;   // The CPV Higgs mixing matrix
547   matrixblock<4> imcvhmix; // Optional: imaginary components
548   matrixblock<3> imau,imad,imae; // Im{} of AU, AD, AE
549   block<double> imhmix;
550   block<double> immsoft;
551
552   //CPV + FLV Input
553   matrixblock<3> immsq2in;  // The Im{} input upper off-diagonal msq2
554   matrixblock<3> immsu2in;  // The Im{} input upper off-diagonal msu2
555   matrixblock<3> immsd2in;  // The Im{} input upper off-diagonal msd2
556   matrixblock<3> immsl2in;  // The Im{} input upper off-diagonal msl2
557   matrixblock<3> immse2in;  // The Im{} input upper off-diagonal mse2
558   matrixblock<3> imtuin,imtdin,imtein; //  The Im{} input upper off-diagonal T
559   //CPV + FLV Output
560   matrixblock<3> imvckm;  // The output DRbar running Im{VCKM} at Q
561   matrixblock<3> imupmns; // The output DRbar running Im{UPMNS} at Q
562   matrixblock<3> immsq2;  // The output DRbar running msq2 at Q
563   matrixblock<3> immsu2;  // The output DRbar running msu2 at Q
564   matrixblock<3> immsd2;  // The output DRbar running msd2 at Q
565   matrixblock<3> immsl2;  // The output DRbar running msl2 at Q
566   matrixblock<3> immse2;  // The output DRbar running mse2 at Q
567   matrixblock<3> imtu,imtd,imte; // Im{} of TU, TD, TE
568   matrixblock<6> imusqmix;// The Im{} up squark mixing matrix
569   matrixblock<6> imdsqmix; // The Im{} down squark mixing matrix
570   matrixblock<6> imselmix; // The Im{} selectron mixing matrix
571   matrixblock<3> imsnumix; // The Im{} sneutrino mixing matrix
572   matrixblock<4> imnmix;   // The Im{} neutralino mixing matrix
573   matrixblock<4> imumix;   // The Im{} chargino L mixing matrix
574   matrixblock<4> imvmix;   // The Im{} chargino R mixing matrix
575
576   //NMSSM Input
577   //    All input is in EXTPAR
578   //NMSSM Output
579   block<double> nmssmrun;  // The block of NMSSM running parameters
580   matrixblock<3> nmhmix;   // The NMSSM scalar Higgs mixing
581   matrixblock<3> nmamix;   // The NMSSM pseudoscalar Higgs mixing
582   matrixblock<5> nmnmix;   // The NMSSM neutralino mixing
583   matrixblock<5> imnmnmix; //   Im{} (for future use)
584
585   //*************************** SET BLOCK VALUE ****************************//
586   template <class T> int set(string,T);
587   template <class T> int set(string,int,T);
588   template <class T> int set(string,int,int,T);
589   template <class T> int set(string,int,int,int,T);
590
591   //***************************** SLHA PRIVATE *****************************//
592 private:
593   //SLHA I/O
594   void message(int, string,string ,int line=0);
595   int verbose;
596   bool headerPrinted, footerPrinted;
597   bool slhaRead, lhefRead, lhefSlha;
598
599 };
600
601 #endif
602
603