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}
t
Test Observable.

◆ 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}