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