a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Functions
RunnerTHDMW.cpp File Reference

Go to the source code of this file.

Functions

int JacobianTHDMW (double t, const double y[], double *dfdy, double dfdt[], void *order)
 
int RGEcheckcustodialMW (const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
 
int RGEcheckMW (const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
 
int RGEcheckTHDMW (const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
 
int RGEscustodialMW (double t, const double y[], double beta[], void *flags)
 
int RGEsMW (double t, const double y[], double beta[], void *flags)
 
int RGEsTHDMW (double t, const double y[], double beta[], void *flags)
 

Function Documentation

◆ JacobianTHDMW()

int JacobianTHDMW ( double  t,
const double  y[],
double *  dfdy,
double  dfdt[],
void *  order 
)

Definition at line 329 of file RunnerTHDMW.cpp.

330 {
331  return 0;
332 }

◆ RGEcheckcustodialMW()

int RGEcheckcustodialMW ( const double  InitialValues[],
const double  t1,
const double  Rpeps,
const double  tNLOuni 
)

Definition at line 715 of file RunnerTHDMW.cpp.

716 {
717  int check=0;
718 
719  //perturbativity of the quartic Higgs couplings
720  for(int i=0;i<7;i++)
721  {
722  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
723  }
724 
725  double la1Q = InitialValues[0];
726  double nu1Q = InitialValues[1];
727  double nu2Q = InitialValues[2];
728  double nu4Q = InitialValues[3];
729  double mu1Q = InitialValues[4];
730  double mu3Q = InitialValues[5];
731  double mu4Q = InitialValues[6];
732 
733  //positivity checks
734  double muAtimes2=4.0*mu1Q+2.0*mu3Q+4.0*mu4Q;
735  if(la1Q<0.0) check=1;
736  if(muAtimes2<0.0) check=1;
737  if(5.0*mu1Q+3.0*mu3Q+3.0*mu4Q-fabs(mu1Q)<0.0) check=1;
738  if(sqrt(la1Q*muAtimes2)+nu1Q+nu2Q-fabs(nu2Q)<0.0) check=1;
739  if(la1Q+0.25*muAtimes2+nu1Q+2.0*nu2Q-2.0*fabs(nu4Q)/sqrt(3.0)<0.0) check=1;
740 
741  //unitarity checks
742 
743  double pi=M_PI;
744  gslpp::matrix<gslpp::complex> Smatrix1(2,2,0.), Smatrix2(2,2,0.);
745  gslpp::matrix<gslpp::complex> Seigenvectors1(2,2,0.), Seigenvectors2(2,2,0.);
746  gslpp::vector<double> Seigenvalues1(2,0.), Seigenvalues2(2,0.);
747  gslpp::vector<gslpp::complex> unitarityeigenvalues(5,0.);
748 
749  if(t1>tNLOuni)
750  {
751 
752  //LO part
753  Smatrix1.assign(0,0, 3.0*la1Q/(16.0*pi));
754  Smatrix1.assign(0,1, (2.0*nu1Q+nu2Q)/(8.0*sqrt(2.0)*pi));
755  Smatrix1.assign(1,0, Smatrix1(0,1));
756  Smatrix1.assign(1,1, (26.0*mu1Q+17.0*mu3Q+13.0*mu4Q)/(32.0*pi));
757 
758  Smatrix2.assign(0,0, la1Q/(16.0*pi));
759  Smatrix2.assign(0,1, nu2Q/(8.0*sqrt(2.0)*pi));
760  Smatrix2.assign(1,0, Smatrix2(0,1));
761  Smatrix2.assign(1,1, (14.0*mu1Q+3.0*mu3Q+27.0*mu4Q)/(96.0*pi));
762 
763  Smatrix1.eigensystem(Seigenvectors1, Seigenvalues1);
764  Smatrix2.eigensystem(Seigenvectors2, Seigenvalues2);
765 
766  for (int i=0; i < 2; i++) {
767  unitarityeigenvalues.assign(i, Seigenvalues1(i));
768  unitarityeigenvalues.assign(2+i, Seigenvalues2(i));
769  }
770  unitarityeigenvalues.assign(4, sqrt(15.0)*nu4Q/(16.0*pi));
771 
772  //beta_la1*16pi^2
773  double betala1 = 12.0*la1Q*la1Q + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 8.0*nu2Q*nu2Q;
774  //beta_mu1*16pi^2
775  double betamu1 = 13.0*mu1Q*mu1Q + 6.0*mu1Q*mu3Q + 6.0*mu1Q*mu4Q + 3.0*nu4Q*nu4Q;
776  //beta_mu3*16pi^2
777  double betamu3 = (134.0*mu1Q*mu1Q + 6.0*mu1Q*(39.0*mu3Q + 22.0*mu4Q)
778  + 3.0*(30.0*mu3Q*mu3Q + 39.0*mu3Q*mu4Q + 9.0*mu4Q*mu4Q
779  + 3.0*nu1Q*nu1Q + 3.0*nu1Q*nu2Q - 5.0*nu4Q*nu4Q))/4.5;
780  //beta_mu4*16pi^2
781  double betamu4 = (4.0*mu1Q*mu1Q + 156.0*mu1Q*mu4Q + 54.0*mu3Q*mu4Q + 144.0*mu4Q*mu4Q
782  + 9.0*nu2Q*nu2Q + 6.0*nu4Q*nu4Q)/9.0;
783  //beta_nu1*16pi^2
784  double betanu1 = (18.0*la1Q*nu1Q
785  + 78.0*mu1Q*nu1Q + 51.0*mu3Q*nu1Q + 39.0*mu4Q*nu1Q + 6.0*nu1Q*nu1Q
786  + 6.0*la1Q*nu2Q + 32.0*mu1Q*nu2Q + 24.0*mu3Q*nu2Q + 6.0*mu4Q*nu2Q
787  + 6.0*nu2Q*nu2Q + 10.0*nu4Q*nu4Q)/3.0;
788  //beta_nu2*16pi^2
789  double betanu2 = 2.0*la1Q*nu2Q + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*nu2Q
790  + 4.0*nu1Q*nu2Q + 6.0*nu2Q*nu2Q + (25.0*nu4Q*nu4Q)/3.0;
791  //beta_nu4*16pi^2
792  double betanu4 = 11.0*mu1Q*nu4Q + 3.0*mu3Q*nu4Q + 9.0*mu4Q*nu4Q + 3.0*nu1Q*nu4Q + 9.0*nu2Q*nu4Q;
793 
794 // diagonalization
795  gslpp::matrix<gslpp::complex> Sbmatrix1(2,2,0.), Sbmatrix2(2,2,0.);
796  gslpp::matrix<gslpp::complex> Seigenvectors1T(2,2,0.), Seigenvectors2T(2,2,0.);
797  gslpp::vector<gslpp::complex> Sbeigenvalues1(2,0.), Sbeigenvalues2(2,0.);
798  gslpp::vector<gslpp::complex> betaeigenvalues(5,0.);
799  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(5,0.);
800 
801  Sbmatrix1.assign(0,0, 3.0*betala1/(16.0*pi));
802  Sbmatrix1.assign(0,1, (2.0*betanu1+betanu2)/(8.0*sqrt(2.0)*pi));
803  Sbmatrix1.assign(1,0, Sbmatrix1(0,1));
804  Sbmatrix1.assign(1,1, (26.0*betamu1+17.0*betamu3+13.0*betamu4)/(32.0*pi));
805 
806  Sbmatrix2.assign(0,0, betala1/(16.0*pi));
807  Sbmatrix2.assign(0,1, betanu2/(8.0*sqrt(2.0)*pi));
808  Sbmatrix2.assign(1,0, Sbmatrix2(0,1));
809  Sbmatrix2.assign(1,1, (14.0*betamu1+3.0*betamu3+27.0*betamu4)/(96.0*pi));
810 
811  Seigenvectors1T=Seigenvectors1.hconjugate();
812  Seigenvectors2T=Seigenvectors2.hconjugate();
813 
814  for (int i=0; i < 2; i++) {
815  for (int k=0; k < 2; k++) {
816  for (int l=0; l < 2; l++) {
817  Sbeigenvalues1.assign(i, Sbeigenvalues1(i) + Seigenvectors1T(i,k) * Sbmatrix1(k,l) * Seigenvectors1(l,i) );
818  Sbeigenvalues2.assign(i, Sbeigenvalues2(i) + Seigenvectors2T(i,k) * Sbmatrix2(k,l) * Seigenvectors2(l,i) );
819  }
820  }
821  betaeigenvalues.assign(i, -1.5 * Sbeigenvalues1(i));
822  betaeigenvalues.assign(i+2, -1.5 * Sbeigenvalues2(i));
823  }
824 
825  betaeigenvalues.assign(4, -1.5 * sqrt(15.0)*betanu4/(16.0*pi));
826 
827  for (int i=0; i < 5; i++) {
828  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
829  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
830  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
831  }
832 
833  } //end of the if(t1>tNLOuni)
834  return check;
835 }

◆ RGEcheckMW()

int RGEcheckMW ( const double  InitialValues[],
const double  t1,
const double  Rpeps,
const double  tNLOuni 
)

Definition at line 550 of file RunnerTHDMW.cpp.

551 {
552  int check=0;
553 
554  //perturbativity of the quartic Higgs couplings
555  for(int i=0;i<12;i++)
556  {
557  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
558  }
559 
560  double la1Q = InitialValues[0];
561  double nu1Q = InitialValues[1];
562  double nu2Q = InitialValues[2];
563  double nu3Q = InitialValues[3];
564  double nu4Q = InitialValues[4];
565  double nu5Q = InitialValues[5];
566  double mu1Q = InitialValues[6];
567  double mu2Q = InitialValues[7];
568  double mu3Q = InitialValues[8];
569  double mu4Q = InitialValues[9];
570  double mu5Q = InitialValues[10];
571  double mu6Q = InitialValues[11];
572 
573  //positivity checks
574  double muAtimes2=mu1Q+mu2Q+mu6Q+2.0*(mu3Q+mu4Q+mu5Q);
575  if(la1Q<0.0) check=1;
576  if(muAtimes2<0.0) check=1;
577  if(mu1Q+mu2Q+mu3Q+mu4Q<0.0) check=1;
578  if(14.0*(mu1Q+mu2Q)+5.0*mu6Q+24.0*(mu3Q+mu4Q)-3.0*fabs(2.0*(mu1Q+mu2Q)-mu6Q)<0.0) check=1;
579  if(5.0*(mu1Q+mu2Q+mu6Q)+6.0*(2.0*mu3Q+mu4Q+mu5Q)-fabs(mu1Q+mu2Q+mu6Q)<0.0) check=1;
580  if(sqrt(la1Q*muAtimes2)+nu1Q<0.0) check=1;
581  if(sqrt(la1Q*muAtimes2)+nu1Q+nu2Q-2.0*fabs(nu3Q)<0.0) check=1;
582  if(la1Q+0.25*muAtimes2+nu1Q+nu2Q+2.0*nu3Q-fabs(nu4Q+nu5Q)/sqrt(3.0)<0.0) check=1;
583 
585  //NLO unitarity//
587 
588  double pi=M_PI;
589  gslpp::vector<gslpp::complex> unitarityeigenvalues(8,0.);
590 
591  if(t1>tNLOuni)
592  {
593  //LO part
594  double muA = 4.0*mu1Q+4.0*mu2Q+8.5*mu3Q+5.0*mu4Q+1.5*mu5Q+2.5*mu6Q;
595  double muB = (4.0*mu1Q+4.0*mu2Q+1.5*mu3Q+12.0*mu4Q+1.5*mu5Q-0.5*mu6Q)/3.0;
596  double muC = (-0.5*mu1Q-0.5*mu2Q+1.5*mu3Q+1.5*mu4Q+12.0*mu5Q+4.0*mu6Q)/3.0;
597  double MA1 = 3.0*la1Q + muA - sqrt(9.0*la1Q*la1Q-6.0*la1Q*muA+muA*muA+32.0*nu1Q*nu1Q+32.0*nu1Q*nu2Q+8.0*nu2Q*nu2Q);
598  double MA2 = 3.0*la1Q + muA + sqrt(9.0*la1Q*la1Q-6.0*la1Q*muA+muA*muA+32.0*nu1Q*nu1Q+32.0*nu1Q*nu2Q+8.0*nu2Q*nu2Q);
599  double MB1 = la1Q + muB - sqrt(la1Q*la1Q-2.0*la1Q*muB+muB*muB+8.0*nu2Q*nu2Q);
600  double MB2 = la1Q + muB + sqrt(la1Q*la1Q-2.0*la1Q*muB+muB*muB+8.0*nu2Q*nu2Q);
601  double MC1 = la1Q + muC - sqrt(la1Q*la1Q-2.0*la1Q*muC+muC*muC+32.0*nu3Q*nu3Q);
602  double MC2 = la1Q + muC + sqrt(la1Q*la1Q-2.0*la1Q*muC+muC*muC+32.0*nu3Q*nu3Q);
603  unitarityeigenvalues.assign(0, MA1/(32.0*pi));
604  unitarityeigenvalues.assign(1, MA2/(32.0*pi));
605  unitarityeigenvalues.assign(2, MB1/(32.0*pi));
606  unitarityeigenvalues.assign(3, MB2/(32.0*pi));
607  unitarityeigenvalues.assign(4, MC1/(32.0*pi));
608  unitarityeigenvalues.assign(5, MC2/(32.0*pi));
609  unitarityeigenvalues.assign(6, la1Q/(16.0*pi));
610  unitarityeigenvalues.assign(7, sqrt(15.0)*(nu4Q+nu5Q)/(64.0*pi));
611 
612  //NLO part
613  //beta_la1*16pi^2
614  double betala1 = 12.0*la1Q*la1Q + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 4.0*nu2Q*nu2Q + 16.0*nu3Q*nu3Q;
615  //beta_nu1*16pi^2
616  double betanu1 = 2.0*nu1Q*nu1Q + nu2Q*nu2Q + 4.0*nu3Q*nu3Q + 2.0*la1Q*(3.0*nu1Q+nu2Q)
617  + (7.0*nu4Q*nu4Q - 4.0*nu4Q*nu5Q + 7.0*nu5Q*nu5Q)/3.0
618  + nu1Q*(8.0*mu1Q + 8.0*mu2Q + 17.0*mu3Q + 10.0*mu4Q + 3.0*mu5Q + 5.0*mu6Q)
619  + nu2Q*(8.0*mu1Q + 8.0*mu2Q + 24.0*mu3Q + 3.0*mu4Q
620  + 3.0*mu5Q + 8.0*mu6Q)/3.0;
621  //beta_nu2*16pi^2
622  double betanu2 = 2.0*nu2Q*nu2Q + 4.0*nu1Q*nu2Q + 16.0*nu3Q*nu3Q + 2.0*la1Q*nu2Q
623  + (4.0*nu4Q*nu4Q + 17.0*nu4Q*nu5Q + 4.0*nu5Q*nu5Q)/3.0
624  + nu2Q*(8.0*mu1Q + 8.0*mu2Q + 3.0*mu3Q + 24.0*mu4Q
625  + 3.0*mu5Q - mu6Q)/3.0;
626  //beta_nu3*16pi^2
627  double betanu3 = 2.0*nu3Q*(la1Q + 2.0*nu1Q + 3.0*nu2Q)
628  + (17.0*nu4Q*nu4Q + 16.0*nu4Q*nu5Q + 17.0*nu5Q*nu5Q)/12.0
629  + nu3Q*(-mu1Q - mu2Q + 3.0*mu3Q + 3.0*mu4Q
630  + 24.0*mu5Q + 8.0*mu6Q)/3.0;
631  //beta_nu4*16pi^2
632  double betanu4 = 8.0*nu3Q*nu4Q + 2.0*nu3Q*nu5Q
633  + nu5Q*(2.0*nu2Q - mu2Q + 2.0*mu4Q + 4.0*mu5Q + mu6Q)
634  + nu4Q*(3.0*nu1Q + 2.0*nu2Q + 6.0*mu1Q + 2.0*mu2Q + 3.0*mu3Q
635  + 2.0*mu4Q + mu5Q + mu6Q);
636  //beta_nu5*16pi^2
637  double betanu5 = 2.0*nu3Q*nu4Q + 8.0*nu3Q*nu5Q
638  + nu4Q*(2.0*nu2Q - mu1Q + 2.0*mu4Q + 4.0*mu5Q + mu6Q)
639  + nu5Q*(3.0*nu1Q + 2.0*nu2Q + 6.0*mu1Q + 2.0*mu2Q + 3.0*mu3Q
640  + 2.0*mu4Q + mu5Q + mu6Q);
641  //beta_mu1*16pi^2
642  double betamu1 = 3.0*nu4Q*nu4Q + 7.0*mu1Q*mu1Q
643  + mu1Q*(6.0*mu2Q + 6.0*mu3Q + 4.0*mu4Q - mu5Q - 2.0*mu6Q)
644  + mu2Q*(4.0*mu4Q - mu5Q)
645  - 2.0*mu4Q*mu6Q + 2.0*mu5Q*mu6Q + mu6Q*mu6Q;
646  //beta_mu2*16pi^2
647  double betamu2 = 3.0*nu5Q*nu5Q + 7.0*mu2Q*mu2Q
648  + mu2Q*(6.0*mu1Q + 6.0*mu3Q + 4.0*mu4Q - mu5Q - 2.0*mu6Q)
649  + mu1Q*(4.0*mu4Q - mu5Q)
650  - 2.0*mu4Q*mu6Q + 2.0*mu5Q*mu6Q + mu6Q*mu6Q;
651  //beta_mu3*16pi^2
652  double betamu3 = 20.0*mu3Q*mu3Q
653  + mu3Q*(288.0*mu1Q + 288.0*mu2Q + 360.0*mu4Q + 108.0*mu5Q + 180.0*mu6Q)/18.0
654  + (36.0*nu1Q*nu1Q + 36.0*nu1Q*nu2Q - 24.0*nu4Q*nu4Q - 12.0*nu4Q*nu5Q
655  - 24.0*nu5Q*nu5Q + 62.0*mu1Q*mu1Q + 64.0*mu1Q*mu2Q + 62.0*mu2Q*mu2Q
656  + (96.0*mu4Q + 18.0*mu5Q + 58.0*mu6Q)*(mu1Q + mu2Q)
657  + 54.0*mu4Q*mu4Q + 36.0*mu4Q*mu5Q + 132.0*mu4Q*mu6Q + 18.0*mu5Q*mu5Q
658  + 18.0*mu5Q*mu6Q + 29.0*mu6Q*mu6Q)/18.0;
659  //beta_mu4*16pi^2
660  double betamu4 = nu2Q*nu2Q - (nu4Q*nu4Q - 4.0*nu4Q*nu5Q + nu5Q*nu5Q)/3.0 + 10.0*mu4Q*mu4Q /*mu4Q??*/
661  + mu5Q*(mu1Q + mu2Q + mu6Q)
662  + mu4Q*(4.0*(4.0*mu1Q + 4.0*mu2Q + mu6Q)/3.0 + 2.0*mu5Q + 6.0*mu4Q)
663  + 4.0*mu5Q*mu5Q
664  + (mu1Q*mu1Q + mu2Q*mu2Q - 4.0*mu6Q*(mu1Q+mu2Q) - 2.0*mu6Q*mu6Q)/9.0
665  + 26.0/9.0*mu1Q*mu2Q;
666  //beta_mu5*16pi^2
667  double betamu5 = 4.0*nu3Q*nu3Q - (nu4Q*nu4Q - 4.0*nu4Q*nu5Q + nu5Q*nu5Q)/3.0
668  + mu5Q*((mu1Q + mu2Q + 19.0*mu6Q)/3.0 + 8.0*mu4Q + 6.0*mu3Q)
669  + 2.0*mu4Q*mu6Q + 8.0*mu5Q*mu5Q
670  + (mu1Q*mu1Q + mu2Q*mu2Q - 4.0*mu6Q*(mu1Q+mu2Q) + 7.0*mu6Q*mu6Q)/9.0
671  - 10.0/9.0*mu1Q*mu2Q;
672  //beta_mu6*16pi^2
673  double betamu6 = 0.5*mu6Q*mu6Q + 3.0*nu4Q*nu4Q + 3.0*nu5Q*nu5Q
674  - 2.0*(mu1Q*mu1Q + mu2Q*mu2Q) + 6.0*mu5Q*(mu1Q + mu2Q)
675  + 7.0*mu6Q*(mu1Q + mu2Q + mu3Q);
676 
677  double betamuA = 4.0*betamu1+4.0*betamu2+8.5*betamu3+5.0*betamu4+1.5*betamu5+2.5*betamu6;
678  double betamuB = (4.0*betamu1+4.0*betamu2+1.5*betamu3+12.0*betamu4+1.5*betamu5-0.5*betamu6)/3.0;
679  double betamuC = (-0.5*betamu1-0.5*betamu2+1.5*betamu3+1.5*betamu4+12.0*betamu5+4.0*betamu6)/3.0;
680  double betaMA1 = 3.0*betala1 + betamuA
681  - sqrt(9.0*betala1*betala1-6.0*betala1*betamuA+betamuA*betamuA
682  +32.0*betanu1*betanu1+32.0*betanu1*betanu2+8.0*betanu2*betanu2);
683  double betaMA2 = 3.0*betala1 + betamuA
684  + sqrt(9.0*betala1*betala1-6.0*betala1*betamuA+betamuA*betamuA
685  +32.0*betanu1*betanu1+32.0*betanu1*betanu2+8.0*betanu2*betanu2);
686  double betaMB1 = betala1 + betamuB - sqrt(betala1*betala1-2.0*betala1*betamuB+betamuB*betamuB+8.0*betanu2*betanu2);
687  double betaMB2 = betala1 + betamuB + sqrt(betala1*betala1-2.0*betala1*betamuB+betamuB*betamuB+8.0*betanu2*betanu2);
688  double betaMC1 = betala1 + betamuC - sqrt(betala1*betala1-2.0*betala1*betamuC+betamuC*betamuC+32.0*betanu3*betanu3);
689  double betaMC2 = betala1 + betamuC + sqrt(betala1*betala1-2.0*betala1*betamuC+betamuC*betamuC+32.0*betanu3*betanu3);
690 
691  gslpp::vector<gslpp::complex> betaeigenvalues(8,0.);
692  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(8,0.);
693 
694  betaeigenvalues.assign(0, -1.5 * betaMA1/(32.0*pi));
695  betaeigenvalues.assign(1, -1.5 * betaMA2/(32.0*pi));
696  betaeigenvalues.assign(2, -1.5 * betaMB1/(32.0*pi));
697  betaeigenvalues.assign(3, -1.5 * betaMB2/(32.0*pi));
698  betaeigenvalues.assign(4, -1.5 * betaMC1/(32.0*pi));
699  betaeigenvalues.assign(5, -1.5 * betaMC2/(32.0*pi));
700  betaeigenvalues.assign(6, -1.5 * betala1/(16.0*pi));
701  betaeigenvalues.assign(7, -1.5 * sqrt(15.0)*(betanu4+betanu5)/(64.0*pi));
702 
703  for (int i=0; i < 8; i++) {
704  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
705  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
706  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
707  }
708 
709  } //end of the if(t1>tNLOuni)
710 
711 
712  return check;
713 }

◆ RGEcheckTHDMW()

int RGEcheckTHDMW ( const double  InitialValues[],
const double  t1,
const double  Rpeps,
const double  tNLOuni 
)

Definition at line 334 of file RunnerTHDMW.cpp.

335 {
336  int check=0;
337 
338  //perturbativity of the Yukawa couplings
339 // for(int i=3;i<6;i++)
340 // {
341 // if(fabs(InitialValues[i])>sqrt(4.0*M_PI)) check=1;
342 // }
343 
344  //perturbativity of the quartic Higgs couplings
345  for(int i=0;i<15;i++)
346  {
347  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
348  }
349 
350  double la1Q = InitialValues[0];
351  double la2Q = InitialValues[1];
352  double la3Q = InitialValues[2];
353  double la4Q = InitialValues[3];
354  double mu1Q = InitialValues[4];
355  double mu3Q = InitialValues[5];
356  double mu4Q = InitialValues[6];
357  double nu1Q = InitialValues[7];
358  double om1Q = InitialValues[8];
359  double ka1Q = InitialValues[9];
360  double nu2Q = InitialValues[10];
361  double om2Q = InitialValues[11];
362  double ka2Q = InitialValues[12];
363  double nu4Q = InitialValues[13];
364  double om4Q = InitialValues[14];
365 
366  //positivity checks
367  double muAtimes2=4.0*mu1Q+2.0*mu3Q+4.0*mu4Q;
368  if(la1Q<0.0) check=1;
369  if(la2Q<0.0) check=1;
370  if(muAtimes2<0.0) check=1;
371  if(5.0*mu1Q+3.0*mu3Q+3.0*mu4Q-fabs(mu1Q)<0.0) check=1;
372  if(sqrt(la1Q*muAtimes2)+nu1Q<0.0) check=1;
373  if(sqrt(la1Q*muAtimes2)+nu1Q+2.0*nu2Q<0.0) check=1;
374  if(la1Q+0.25*muAtimes2+nu1Q+2.0*nu2Q-2.0*fabs(nu4Q)/sqrt(3.0)<0.0) check=1;
375  if(la3Q+sqrt(la1Q*la2Q)<0.0) check=1;
376  if(la3Q+2.0*la4Q+sqrt(la1Q*la2Q)<0.0) check=1;
377  if(sqrt(la2Q*muAtimes2)+om1Q<0.0) check=1;
378  if(sqrt(la2Q*muAtimes2)+om1Q+2.0*om2Q<0.0) check=1;
379  if(la2Q+0.25*muAtimes2+om1Q+2.0*om2Q-2.0*fabs(om4Q)/sqrt(3.0)<0.0) check=1;
380 
382  //NLO unitarity//
384 
385  double pi=M_PI;
386  gslpp::matrix<gslpp::complex> Smatrix1(4,4,0.), Smatrix2(4,4,0.);
387  gslpp::matrix<gslpp::complex> Seigenvectors1(4,4,0.), Seigenvectors2(4,4,0.);
388  gslpp::vector<double> Seigenvalues1(4,0.), Seigenvalues2(4,0.);
389  gslpp::vector<gslpp::complex> unitarityeigenvalues(11,0.);
390 
391  if(t1>tNLOuni)
392  {
393 
394  //LO part
395  Smatrix1.assign(0,0, 3.0*la1Q/(16.0*pi));
396  Smatrix1.assign(0,1, (2.0*la3Q+la4Q)/(16.0*pi));
397  Smatrix1.assign(1,0, Smatrix1(0,1));
398  Smatrix1.assign(0,3, (2.0*nu1Q+nu2Q)/(8.0*sqrt(2.0)*pi));
399  Smatrix1.assign(3,0, Smatrix1(0,3));
400  Smatrix1.assign(1,1, 3.0*la2Q/(16.0*pi));
401  Smatrix1.assign(1,3, (2.0*om1Q+om2Q)/(8.0*sqrt(2.0)*pi));
402  Smatrix1.assign(3,1, Smatrix1(1,3));
403  Smatrix1.assign(2,2, (la3Q+5.0*la4Q)/(16.0*pi));
404  Smatrix1.assign(2,3, (4.0*ka1Q+2.0*ka2Q)/(16.0*pi));
405  Smatrix1.assign(3,2, Smatrix1(2,3));
406  Smatrix1.assign(3,3, (26.0*mu1Q+17.0*mu3Q+13.0*mu4Q)/(32.0*pi));
407 
408  Smatrix2.assign(0,0, la1Q/(16.0*pi));
409  Smatrix2.assign(0,1, la4Q/(16.0*pi));
410  Smatrix2.assign(1,0, Smatrix2(0,1));
411  Smatrix2.assign(0,3, nu2Q/(8.0*sqrt(2.0)*pi));
412  Smatrix2.assign(3,0, Smatrix2(0,3));
413  Smatrix2.assign(1,1, la2Q/(16.0*pi));
414  Smatrix2.assign(1,3, om2Q/(8.0*sqrt(2.0)*pi));
415  Smatrix2.assign(3,1, Smatrix2(1,3));
416  Smatrix2.assign(2,2, (la3Q+la4Q)/(16.0*pi));
417  Smatrix2.assign(2,3, ka2Q/(8.0*pi));
418  Smatrix2.assign(3,2, Smatrix2(2,3));
419  Smatrix2.assign(3,3, (14.0*mu1Q+3.0*mu3Q+27.0*mu4Q)/(96.0*pi));
420 
421  Smatrix1.eigensystem(Seigenvectors1, Seigenvalues1);
422  Smatrix2.eigensystem(Seigenvectors2, Seigenvalues2);
423 
424  for (int i=0; i < 4; i++) {
425  unitarityeigenvalues.assign(i, Seigenvalues1(i));
426  unitarityeigenvalues.assign(4+i, Seigenvalues2(i));
427  }
428  unitarityeigenvalues.assign(8, (la3Q-la4Q)/(16.0*pi));
429  unitarityeigenvalues.assign(9, sqrt(15.0)*nu4Q/(16.0*pi));
430  unitarityeigenvalues.assign(10, sqrt(15.0)*om4Q/(16.0*pi));
431 
432  //beta_la1*16pi^2
433  double betala1 = 12.0*la1Q*la1Q + 4.0*la3Q*la3Q + 4.0*la3Q*la4Q + 4.0*la4Q*la4Q
434  + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 8.0*nu2Q*nu2Q;
435  //beta_la2*16pi^2
436  double betala2 = 12.0*la2Q*la2Q + 4.0*la3Q*la3Q + 4.0*la3Q*la4Q + 4.0*la4Q*la4Q
437  + 8.0*om1Q*om1Q + 8.0*om1Q*om2Q + 8.0*om2Q*om2Q;
438  //beta_la3*16pi^2
439  double betala3 = 4.0*la3Q*la3Q + 4.0*la4Q*la4Q + (la1Q+la2Q)*(6.0*la3Q+2.0*la4Q)
440  + 8.0*ka2Q*ka2Q + 8.0*nu1Q*om1Q + 4.0*nu2Q*om1Q + 4.0*nu1Q*om2Q;
441  //beta_la4*16pi^2
442  double betala4 = (la1Q*la4Q + la2Q*la4Q + 4.0*la3Q*la4Q + 6.0*la4Q*la4Q
443  + 4.0*ka1Q*ka1Q + 4.0*ka1Q*ka2Q + 2.0*ka2Q*ka2Q + 2.0*nu2Q*om2Q)*2.0;
444  //beta_mu1*16pi^2
445  double betamu1 = 11.0*mu1Q*mu1Q + 3.0*mu1Q*mu4Q + mu1Q*(2.0*mu1Q+6.0*mu3Q+3.0*mu4Q)
446  + 3.0*nu4Q*nu4Q + 3.0*om4Q*om4Q;
447  //beta_mu3*16pi^2
448  double betamu3 = (18.0*ka1Q*ka1Q + 18.0*ka1Q*ka2Q + 134.0*mu1Q*mu1Q + 6.0*mu1Q*(39.0*mu3Q + 22.0*mu4Q)
449  + 3.0*(30.0*mu3Q*mu3Q + 39.0*mu3Q*mu4Q + 9.0*mu4Q*mu4Q
450  + 3.0*nu1Q*nu1Q + 3.0*nu1Q*nu2Q - 5.0*nu4Q*nu4Q
451  + 3.0*om1Q*om1Q + 3.0*om1Q*om2Q - 5.0*om4Q*om4Q))/4.5;
452  //beta_mu4*16pi^2
453  double betamu4 = (18.0*ka2Q*ka2Q + 4.0*mu1Q*mu1Q + 156.0*mu1Q*mu4Q + 54.0*mu3Q*mu4Q + 144.0*mu4Q*mu4Q
454  + 9.0*nu2Q*nu2Q + 6.0*nu4Q*nu4Q + 9.0*om2Q*om2Q + 6.0*om4Q*om4Q)/9.0;
455  //beta_nu1*16pi^2
456  double betanu1 = (6.0*ka1Q*ka1Q + 6.0*ka2Q*ka2Q + 18.0*la1Q*nu1Q
457  + 78.0*mu1Q*nu1Q + 51.0*mu3Q*nu1Q + 39.0*mu4Q*nu1Q + 6.0*nu1Q*nu1Q
458  + 6.0*la1Q*nu2Q + 32.0*mu1Q*nu2Q + 24.0*mu3Q*nu2Q + 6.0*mu4Q*nu2Q
459  + 6.0*nu2Q*nu2Q + 10.0*nu4Q*nu4Q
460  + 12.0*la3Q*om1Q + 6.0*la4Q*om1Q + 6.0*la3Q*om2Q)/3.0;
461  //beta_om1*16pi^2
462  double betaom1 = (6.0*ka1Q*ka1Q + 6.0*ka2Q*ka2Q
463  + 12.0*la3Q*nu1Q + 6.0*la4Q*nu1Q + 6.0*la3Q*nu2Q
464  + 18.0*la2Q*om1Q + 78.0*mu1Q*om1Q + 51.0*mu3Q*om1Q + 39.0*mu4Q*om1Q + 6.0*om1Q*om1Q
465  + 6.0*la2Q*om2Q + 32.0*mu1Q*om2Q + 24.0*mu3Q*om2Q + 6.0*mu4Q*om2Q + 6.0*om2Q*om2Q
466  + 10.0*om4Q*om4Q)/3.0;
467  //beta_ka1*16pi^2
468  double betaka1 = (6.0*ka1Q*(2.0*la3Q + 10.0*la4Q + 18.0*mu1Q + 17.0*mu3Q + 13.0*mu4Q + 2.0*nu1Q + 2.0*om1Q)
469  + ka2Q*(24.0*la4Q + 64.0*mu1Q + 48.0*mu3Q + 24.0*mu4Q + 9.0*nu2Q + 9.0*om2Q)
470  + 20.0*nu4Q*om4Q)/6.0;
471  //beta_nu2*16pi^2
472  double betanu2 = 4.0*ka1Q*ka2Q + 6.0*ka2Q*ka2Q + 2.0*la1Q*nu2Q + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*nu2Q
473  + 4.0*nu1Q*nu2Q + 6.0*nu2Q*nu2Q + (25.0*nu4Q*nu4Q)/3.0 + 2.0*la4Q*om2Q;
474  //beta_om2*16pi^2
475  double betaom2 = 4.0*ka1Q*ka2Q + 6.0*ka2Q*ka2Q + 2.0*la4Q*nu2Q + 2.0*la2Q*om2Q
476  + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*om2Q + 4.0*om1Q*om2Q + 6.0*om2Q*om2Q
477  + (25.0*om4Q*om4Q)/3.0;
478  //beta_ka2*16pi^2
479  double betaka2 = (ka2Q*(6.0*la3Q + 6.0*la4Q + 14.0*mu1Q + 3.0*mu3Q + 27.0*mu4Q
480  + 6.0*nu1Q + 12.0*nu2Q + 6.0*om1Q + 12.0*om2Q)
481  + 6.0*ka1Q*(nu2Q + om2Q) + 42.0*nu4Q*om4Q)/3.0;
482  //beta_nu4*16pi^2
483  double betanu4 = 11.0*mu1Q*nu4Q + 3.0*mu3Q*nu4Q + 9.0*mu4Q*nu4Q + 3.0*nu1Q*nu4Q + 9.0*nu2Q*nu4Q
484  + 3.0*ka1Q*om4Q + 9.0*ka2Q*om4Q;
485  //beta_om4*16pi^2
486  double betaom4 = 3.0*ka1Q*nu4Q + 9.0*ka2Q*nu4Q
487  + (11.0*mu1Q + 3.0*(mu3Q + 3.0*mu4Q + om1Q + 3.0*om2Q))*om4Q;
488 
489 // diagonalization
490  gslpp::matrix<gslpp::complex> Sbmatrix1(4,4,0.), Sbmatrix2(4,4,0.);
491  gslpp::matrix<gslpp::complex> Seigenvectors1T(4,4,0.), Seigenvectors2T(4,4,0.);
492  gslpp::vector<gslpp::complex> Sbeigenvalues1(4,0.), Sbeigenvalues2(4,0.);
493  gslpp::vector<gslpp::complex> betaeigenvalues(11,0.);
494  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(11,0.);
495 
496  Sbmatrix1.assign(0,0, 3.0*betala1/(16.0*pi));
497  Sbmatrix1.assign(0,1, (2.0*betala3+betala4)/(16.0*pi));
498  Sbmatrix1.assign(1,0, Sbmatrix1(0,1));
499  Sbmatrix1.assign(0,3, (2.0*betanu1+betanu2)/(8.0*sqrt(2.0)*pi));
500  Sbmatrix1.assign(3,0, Sbmatrix1(0,3));
501  Sbmatrix1.assign(1,1, 3.0*betala2/(16.0*pi));
502  Sbmatrix1.assign(1,3, (2.0*betaom1+betaom2)/(8.0*sqrt(2.0)*pi));
503  Sbmatrix1.assign(3,1, Sbmatrix1(1,3));
504  Sbmatrix1.assign(2,2, (betala3+5.0*betala4)/(16.0*pi));
505  Sbmatrix1.assign(2,3, (4.0*betaka1+2.0*betaka2)/(16.0*pi));
506  Sbmatrix1.assign(3,2, Sbmatrix1(2,3));
507  Sbmatrix1.assign(3,3, (26.0*betamu1+17.0*betamu3+13.0*betamu4)/(32.0*pi));
508 
509  Sbmatrix2.assign(0,0, betala1/(16.0*pi));
510  Sbmatrix2.assign(0,1, betala4/(16.0*pi));
511  Sbmatrix2.assign(1,0, Sbmatrix2(0,1));
512  Sbmatrix2.assign(0,3, betanu2/(8.0*sqrt(2.0)*pi));
513  Sbmatrix2.assign(3,0, Sbmatrix2(0,3));
514  Sbmatrix2.assign(1,1, betala2/(16.0*pi));
515  Sbmatrix2.assign(1,3, betaom2/(8.0*sqrt(2.0)*pi));
516  Sbmatrix2.assign(3,1, Sbmatrix2(1,3));
517  Sbmatrix2.assign(2,2, (betala3+betala4)/(16.0*pi));
518  Sbmatrix2.assign(2,3, betaka2/(8.0*pi));
519  Sbmatrix2.assign(3,2, Sbmatrix2(2,3));
520  Sbmatrix2.assign(3,3, (14.0*betamu1+3.0*betamu3+27.0*betamu4)/(96.0*pi));
521 
522  Seigenvectors1T=Seigenvectors1.hconjugate();
523  Seigenvectors2T=Seigenvectors2.hconjugate();
524 
525  for (int i=0; i < 4; i++) {
526  for (int k=0; k < 4; k++) {
527  for (int l=0; l < 4; l++) {
528  Sbeigenvalues1.assign(i, Sbeigenvalues1(i) + Seigenvectors1T(i,k) * Sbmatrix1(k,l) * Seigenvectors1(l,i) );
529  Sbeigenvalues2.assign(i, Sbeigenvalues2(i) + Seigenvectors2T(i,k) * Sbmatrix2(k,l) * Seigenvectors2(l,i) );
530  }
531  }
532  betaeigenvalues.assign(i, -1.5 * Sbeigenvalues1(i));
533  betaeigenvalues.assign(i+4, -1.5 * Sbeigenvalues2(i));
534  }
535 
536  betaeigenvalues.assign(8, -1.5 * (betala3-betala4)/(16.0*pi));
537  betaeigenvalues.assign(9, -1.5 * sqrt(15.0)*betanu4/(16.0*pi));
538  betaeigenvalues.assign(10, -1.5 * sqrt(15.0)*betaom4/(16.0*pi));
539 
540  for (int i=0; i < 11; i++) {
541  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
542  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
543  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
544  }
545 
546  } //end of the if(t1>tNLOuni)
547  return check;
548 }

◆ RGEscustodialMW()

int RGEscustodialMW ( double  t,
const double  y[],
double  beta[],
void *  flags 
)

Definition at line 260 of file RunnerTHDMW.cpp.

261 {
262  (void)(t); /* avoid unused parameter warning */
263  int ord = *(int *)flags;
264 // int type = flag-ord;
265 // double Yb1=0;
266 // double Yb2=0;
267 // double Ytau1=0;
268 // double Ytau2=0;
269  double la1=y[0];
270  double nu1=y[1];
271  double nu2=y[2];
272  double nu4=y[3];
273  double mu1=y[4];
274  double mu3=y[5];
275  double mu4=y[6];
276 
277  double pi=M_PI;
278 
279  //RGE taken from 1303.4848
280 
281  //beta_la1
282  //This was taken from the custodial THDMW!
283  beta[0] = (12.0*la1*la1 + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 8.0*nu2*nu2)/(16.0*pi*pi);
284  //beta_nu1
285  beta[1] = (2.0*nu1*nu1 + 2.0*nu2*nu2 + 2.0*la1*(3.0*nu1+nu2)
286  + 10.0*nu4*nu4/3.0
287  + nu1*(26.0*mu1 + 17.0*mu3 + 13.0*mu4)
288  + nu2*(32.0*mu1 + 24.0*mu3 + 6.0*mu4)/3.0)/(16.0*pi*pi);
289  //beta_nu2
290  beta[2] = (4.0*nu1*nu2 + 6.0*nu2*nu2 + 2.0*la1*nu2 + 25.0*nu4*nu4/3.0
291  + nu2*(14.0*mu1 + 3.0*mu3 + 27.0*mu4)/3.0)/(16.0*pi*pi);
292  //beta_nu4
293  beta[3] = (3.0*nu1*nu4 + 9.0*nu2*nu4
294  + nu4*(11.0*mu1 + 3.0*mu3 + 9.0*mu4))/(16.0*pi*pi);
295  //beta_mu1
296  beta[4] = (3.0*nu4*nu4 + 13.0*mu1*mu1
297  + 6.0*mu1*(mu3+mu4))/(16.0*pi*pi);
298  //beta_mu3
299  beta[5] = (20.0*mu3*mu3
300  + mu3*(52.0*mu1 + 26.0*mu4)
301  + 2.0*nu1*nu1 + 2.0*nu1*nu2 - 10.0*nu4*nu4/3.0
302  + 268.0*mu1*mu1/9.0 + 88.0*mu1*mu4/3.0 + 6.0*mu4*mu4)/(16.0*pi*pi);
303  //beta_mu4
304  //This was taken from the custodial THDMW, because I think the formula from 1303.4848 is incorrect
305  beta[6] = (nu2*nu2 + 2.0*nu4*nu4/3.0 + 16.0*mu4*mu4
306  + 52.0*mu1*mu4/3.0 + 6.0*mu3*mu4
307  + 4.0/9.0*mu1*mu1)/(16.0*pi*pi);
308 
309  if(ord==1){
310  //beta_la1
311  beta[0] += 0.0;
312  //beta_nu1
313  beta[1] += 0.0;
314  //beta_nu2
315  beta[2] += 0.0;
316  //beta_nu4
317  beta[3] += 0.0;
318  //beta_mu1
319  beta[4] += 0.0;
320  //beta_mu3
321  beta[5] += 0.0;
322  //beta_mu4
323  beta[6] += 0.0;
324  }
325 
326  return 0;
327 }

◆ RGEsMW()

int RGEsMW ( double  t,
const double  y[],
double  beta[],
void *  flags 
)

Definition at line 139 of file RunnerTHDMW.cpp.

140 {
141  (void)(t); /* avoid unused parameter warning */
142  int ord = *(int *)flags;
143 // int type = flag-ord;
144 // double Yb1=0;
145 // double Yb2=0;
146 // double Ytau1=0;
147 // double Ytau2=0;
148  double la1=y[0];
149  double nu1=y[1];
150  double nu2=y[2];
151  double nu3=y[3];
152  double nu4=y[4];
153  double nu5=y[5];
154  double mu1=y[6];
155  double mu2=y[7];
156  double mu3=y[8];
157  double mu4=y[9];
158  double mu5=y[10];
159  double mu6=y[11];
160 
161  double pi=M_PI;
162 
163  //RGE taken from 1303.4848
164  //The lambda1 running is taken from (4.1) in 1303.4848. The rest stems from the appendix.
165 
166  //beta_la1
167  beta[0] = (12.0*la1*la1 + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 4.0*nu2*nu2 + 16.0*nu3*nu3)/(16.0*pi*pi);
168  //beta_nu1
169  beta[1] = (2.0*nu1*nu1 + nu2*nu2 + 4.0*nu3*nu3 + 2.0*la1*(3.0*nu1+nu2)
170  + (7.0*nu4*nu4 - 4.0*nu4*nu5 + 7.0*nu5*nu5)/3.0
171  + nu1*(8.0*mu1 + 8.0*mu2 + 17.0*mu3 + 10.0*mu4 + 3.0*mu5 + 5.0*mu6)
172  + nu2*(8.0*mu1 + 8.0*mu2 + 24.0*mu3 + 3.0*mu4
173  + 3.0*mu5 + 8.0*mu6)/3.0)/(16.0*pi*pi);
174  //beta_nu2
175  beta[2] = (2.0*nu2*nu2 + 4.0*nu1*nu2 + 16.0*nu3*nu3 + 2.0*la1*nu2
176  + (4.0*nu4*nu4 + 17.0*nu4*nu5 + 4.0*nu5*nu5)/3.0
177  + nu2*(8.0*mu1 + 8.0*mu2 + 3.0*mu3 + 24.0*mu4
178  + 3.0*mu5 - mu6)/3.0)/(16.0*pi*pi);
179  //beta_nu3
180  beta[3] = (2.0*nu3*(la1 + 2.0*nu1 + 3.0*nu2)
181  + (17.0*nu4*nu4 + 16.0*nu4*nu5 + 17.0*nu5*nu5)/12.0
182  + nu3*(-mu1 - mu2 + 3.0*mu3 + 3.0*mu4
183  + 24.0*mu5 + 8.0*mu6)/3.0)/(16.0*pi*pi);
184  //beta_nu4
185  beta[4] = (8.0*nu3*nu4 + 2.0*nu3*nu5
186  + nu5*(2.0*nu2 - mu2 + 2.0*mu4 + 4.0*mu5 + mu6)
187  + nu4*(3.0*nu1 + 2.0*nu2 + 6.0*mu1 + 2.0*mu2 + 3.0*mu3
188  + 2.0*mu4 + mu5 + mu6))/(16.0*pi*pi);
189  //beta_nu5
190  beta[5] = (2.0*nu3*nu4 + 8.0*nu3*nu5
191  + nu4*(2.0*nu2 - mu1 + 2.0*mu4 + 4.0*mu5 + mu6)
192  + nu5*(3.0*nu1 + 2.0*nu2 + 6.0*mu1 + 2.0*mu2 + 3.0*mu3
193  + 2.0*mu4 + mu5 + mu6))/(16.0*pi*pi);
194  //beta_mu1
195  beta[6] = (3.0*nu4*nu4 + 7.0*mu1*mu1
196  + mu1*(6.0*mu2 + 6.0*mu3 + 4.0*mu4 - mu5 - 2.0*mu6)
197  + mu2*(4.0*mu4 - mu5)
198  - 2.0*mu4*mu6 + 2.0*mu5*mu6 + mu6*mu6)/(16.0*pi*pi);
199  //beta_mu2
200  beta[7] = (3.0*nu5*nu5 + 7.0*mu2*mu2
201  + mu2*(6.0*mu1 + 6.0*mu3 + 4.0*mu4 - mu5 - 2.0*mu6)
202  + mu1*(4.0*mu4 - mu5)
203  - 2.0*mu4*mu6 + 2.0*mu5*mu6 + mu6*mu6)/(16.0*pi*pi);
204  //beta_mu3
205  beta[8] = (20.0*mu3*mu3
206  + mu3*(288.0*mu1 + 288.0*mu2 + 360.0*mu4 + 108.0*mu5 + 180.0*mu6)/18.0
207  + (36.0*nu1*nu1 + 36.0*nu1*nu2 - 24.0*nu4*nu4 - 12.0*nu4*nu5
208  - 24.0*nu5*nu5 + 62.0*mu1*mu1 + 64.0*mu1*mu2 + 62.0*mu2*mu2
209  + (96.0*mu4 + 18.0*mu5 + 58.0*mu6)*(mu1 + mu2)
210  + 54.0*mu4*mu4 + 36.0*mu4*mu5 + 132.0*mu4*mu6 + 18.0*mu5*mu5
211  + 18.0*mu5*mu6 + 29.0*mu6*mu6)/18.0)/(16.0*pi*pi);
212  //beta_mu4
213  beta[9] = (nu2*nu2 - (nu4*nu4 - 4.0*nu4*nu5 + nu5*nu5)/3.0 + 10.0*mu4*mu4 /*mu4??*/
214  + mu5*(mu1 + mu2 + mu6)
215  + mu4*(4.0*(4.0*mu1 + 4.0*mu2 + mu6)/3.0 + 2.0*mu5 + 6.0*mu4)
216  + 4.0*mu5*mu5
217  + (mu1*mu1 + mu2*mu2 - 4.0*mu6*(mu1+mu2) - 2.0*mu6*mu6)/9.0
218  + 26.0/9.0*mu1*mu2)/(16.0*pi*pi);
219  //beta_mu5
220  beta[10] = (4.0*nu3*nu3 - (nu4*nu4 - 4.0*nu4*nu5 + nu5*nu5)/3.0
221  + mu5*((mu1 + mu2 + 19.0*mu6)/3.0 + 8.0*mu4 + 6.0*mu3)
222  + 2.0*mu4*mu6 + 8.0*mu5*mu5
223  + (mu1*mu1 + mu2*mu2 - 4.0*mu6*(mu1+mu2) + 7.0*mu6*mu6)/9.0
224  - 10.0/9.0*mu1*mu2)/(16.0*pi*pi);
225  //beta_mu6
226  beta[11] = (0.5*mu6*mu6 + 3.0*nu4*nu4 + 3.0*nu5*nu5
227  - 2.0*(mu1*mu1 + mu2*mu2) + 6.0*mu5*(mu1 + mu2)
228  + 7.0*mu6*(mu1 + mu2 + mu3))/(16.0*pi*pi);
229 
230  if(ord==1){
231  //beta_la1
232  beta[0] += 0.0;
233  //beta_nu1
234  beta[1] += 0.0;
235  //beta_nu2
236  beta[2] += 0.0;
237  //beta_nu3
238  beta[3] += 0.0;
239  //beta_nu4
240  beta[4] += 0.0;
241  //beta_nu5
242  beta[5] += 0.0;
243  //beta_mu1
244  beta[6] += 0.0;
245  //beta_mu2
246  beta[7] += 0.0;
247  //beta_mu3
248  beta[8] += 0.0;
249  //beta_mu4
250  beta[9] += 0.0;
251  //beta_mu5
252  beta[10] += 0.0;
253  //beta_mu6
254  beta[11] += 0.0;
255  }
256 
257  return 0;
258 }

◆ RGEsTHDMW()

int RGEsTHDMW ( double  t,
const double  y[],
double  beta[],
void *  flags 
)

Definition at line 19 of file RunnerTHDMW.cpp.

20 {
21  (void)(t); /* avoid unused parameter warning */
22  int ord = *(int *)flags;
23 // int type = flag-ord;
24 // double Yb1=0;
25 // double Yb2=0;
26 // double Ytau1=0;
27 // double Ytau2=0;
28  double la1=y[0];
29  double la2=y[1];
30  double la3=y[2];
31  double la4=y[3];
32  double mu1=y[4];
33  double mu3=y[5];
34  double mu4=y[6];
35  double nu1=y[7];
36  double om1=y[8];
37  double ka1=y[9];
38  double nu2=y[10];
39  double om2=y[11];
40  double ka2=y[12];
41  double nu4=y[13];
42  double om4=y[14];
43 
44  double pi=M_PI;
45 
46  //beta_la1
47  beta[0] = (12.0*la1*la1 + 4.0*la3*la3 + 4.0*la3*la4 + 4.0*la4*la4
48  + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 8.0*nu2*nu2)/(16.0*pi*pi);
49  //beta_la2
50  beta[1] = (12.0*la2*la2 + 4.0*la3*la3 + 4.0*la3*la4 + 4.0*la4*la4
51  + 8.0*om1*om1 + 8.0*om1*om2 + 8.0*om2*om2)/(16.0*pi*pi);
52  //beta_la3
53  beta[2] = (4.0*la3*la3 + 4.0*la4*la4 + (la1+la2)*(6.0*la3+2.0*la4)
54  + 8.0*ka2*ka2 + 8.0*nu1*om1 + 4.0*nu2*om1 + 4.0*nu1*om2)/(16.0*pi*pi);
55  //beta_la4
56  beta[3] = (la1*la4 + la2*la4 + 4.0*la3*la4 + 6.0*la4*la4
57  + 4.0*ka1*ka1 + 4.0*ka1*ka2 + 2.0*ka2*ka2 + 2.0*nu2*om2)/(8.0*pi*pi);
58  //beta_mu1
59  beta[4] = (11.0*mu1*mu1 + 3.0*mu1*mu4 + mu1*(2.0*mu1+6.0*mu3+3.0*mu4)
60  + 3.0*nu4*nu4 + 3.0*om4*om4)/(16.0*pi*pi);
61  //beta_mu3
62  beta[5] = (18.0*ka1*ka1 + 18.0*ka1*ka2 + 134.0*mu1*mu1 + 6.0*mu1*(39.0*mu3 + 22.0*mu4)
63  + 3.0*(30.0*mu3*mu3 + 39.0*mu3*mu4 + 9.0*mu4*mu4
64  + 3.0*nu1*nu1 + 3.0*nu1*nu2 - 5.0*nu4*nu4
65  + 3.0*om1*om1 + 3.0*om1*om2 - 5.0*om4*om4))/(72.0*pi*pi);
66  //beta_mu4
67  beta[6] = (18.0*ka2*ka2 + 4.0*mu1*mu1 + 156.0*mu1*mu4 + 54.0*mu3*mu4 + 144.0*mu4*mu4
68  + 9.0*nu2*nu2 + 6.0*nu4*nu4 + 9.0*om2*om2 + 6.0*om4*om4)/(144.0*pi*pi);
69  //beta_nu1
70  beta[7] = (6.0*ka1*ka1 + 6.0*ka2*ka2 + 18.0*la1*nu1
71  + 78.0*mu1*nu1 + 51.0*mu3*nu1 + 39.0*mu4*nu1 + 6.0*nu1*nu1
72  + 6.0*la1*nu2 + 32.0*mu1*nu2 + 24.0*mu3*nu2 + 6.0*mu4*nu2
73  + 6.0*nu2*nu2 + 10.0*nu4*nu4
74  + 12.0*la3*om1 + 6.0*la4*om1 + 6.0*la3*om2)/(48.0*pi*pi);
75  //beta_om1
76  beta[8] = (6.0*ka1*ka1 + 6.0*ka2*ka2
77  + 12.0*la3*nu1 + 6.0*la4*nu1 + 6.0*la3*nu2
78  + 18.0*la2*om1 + 78.0*mu1*om1 + 51.0*mu3*om1 + 39.0*mu4*om1 + 6.0*om1*om1
79  + 6.0*la2*om2 + 32.0*mu1*om2 + 24.0*mu3*om2 + 6.0*mu4*om2 + 6.0*om2*om2
80  + 10.0*om4*om4)/(48.0*pi*pi);
81  //beta_ka1
82  beta[9] = (6.0*ka1*(2.0*la3 + 10.0*la4 + 18.0*mu1 + 17.0*mu3 + 13.0*mu4 + 2.0*nu1 + 2.0*om1)
83  + ka2*(24.0*la4 + 64.0*mu1 + 48.0*mu3 + 24.0*mu4 + 9.0*nu2 + 9.0*om2)
84  + 20.0*nu4*om4)/(96.0*pi*pi);
85  //beta_nu2
86  beta[10] = (4.0*ka1*ka2 + 6.0*ka2*ka2 + 2.0*la1*nu2 + ((14.0*mu1)/3.0 + mu3 + 9.0*mu4)*nu2
87  + 4.0*nu1*nu2 + 6.0*nu2*nu2 + (25.0*nu4*nu4)/3.0 + 2.0*la4*om2)/(16.0*pi*pi);
88  //beta_om2
89  beta[11] = (4.0*ka1*ka2 + 6.0*ka2*ka2 + 2.0*la4*nu2 + 2.0*la2*om2
90  + ((14.0*mu1)/3.0 + mu3 + 9.0*mu4)*om2 + 4.0*om1*om2 + 6.0*om2*om2
91  + (25.0*om4*om4)/3.0)/(16.0*pi*pi);
92  //beta_ka2
93  beta[12] = (ka2*(6.0*la3 + 6.0*la4 + 14.0*mu1 + 3.0*mu3 + 27.0*mu4
94  + 6.0*nu1 + 12.0*nu2 + 6.0*om1 + 12.0*om2)
95  + 6.0*ka1*(nu2 + om2) + 42.0*nu4*om4)/(48.0*pi*pi);
96  //beta_nu4
97  beta[13] = (11.0*mu1*nu4 + 3.0*mu3*nu4 + 9.0*mu4*nu4 + 3.0*nu1*nu4 + 9.0*nu2*nu4
98  + 3.0*ka1*om4 + 9.0*ka2*om4)/(16.0*pi*pi);
99  //beta_om4
100  beta[14] = (3.0*ka1*nu4 + 9.0*ka2*nu4
101  + (11.0*mu1 + 3.0*(mu3 + 3.0*mu4 + om1 + 3.0*om2))*om4)/(16.0*pi*pi);
102 
103  if(ord==1){
104  //beta_la1
105  beta[0] += 0.0;
106  //beta_la2
107  beta[1] += 0.0;
108  //beta_la3
109  beta[2] += 0.0;
110  //beta_la4
111  beta[3] += 0.0;
112  //beta_mu1
113  beta[4] += 0.0;
114  //beta_mu3
115  beta[5] += 0.0;
116  //beta_mu4
117  beta[6] += 0.0;
118  //beta_nu1
119  beta[7] += 0.0;
120  //beta_om1
121  beta[8] += 0.0;
122  //beta_ka1
123  beta[9] += 0.0;
124  //beta_nu2
125  beta[10] += 0.0;
126  //beta_om2
127  beta[11] += 0.0;
128  //beta_ka2
129  beta[12] += 0.0;
130  //beta_nu4
131  beta[13] += 0.0;
132  //beta_om4
133  beta[14] += 0.0;
134  }
135 
136  return 0;
137 }
gslpp::matrix< gslpp::complex >
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
gslpp::vector< gslpp::complex >