JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
RooqqZZ_JHU.cc
Go to the documentation of this file.
1 // Included verbatim from UserCode/scasasso/HZZ4lAnalysis/HZZ4lCommon/src tag V00-00-00
2 
3 
4 /*****************************************************************************
5  * Project: RooFit *
6  * *
7  * This code was autogenerated by RooClassFactory *
8  *****************************************************************************/
9 
10 // Your description goes here...
11 
12 #include <math.h>
13 #include "RooqqZZ_JHU.h"
14 #include "Riostream.h"
15 #include "RooAbsReal.h"
16 #include "RooAbsCategory.h"
17 #include "TMath.h"
18 
19 //ClassImp(RooqqZZ_JHU)
20 
21 RooqqZZ_JHU::RooqqZZ_JHU(const char *name, const char *title,
22  RooAbsReal& _m1j,
23  RooAbsReal& _m2j,
24  RooAbsReal& _h1j,
25  RooAbsReal& _h2j,
26  RooAbsReal& _phij,
27  RooAbsReal& _hsj,
28  RooAbsReal& _Phi1j,
29  RooAbsReal& _mZZj) :
30 RooAbsPdf(name,title),
31 m1j("m1j","m1j",this,_m1j),
32 m2j("m2j","m2j",this,_m2j),
33 h1j("h1j","h1j",this,_h1j),
34 h2j("h2j","h2j",this,_h2j),
35 phij("phij","phij",this,_phij),
36 hsj("hsj","hsj",this,_hsj),
37 Phi1j("Phi1j","Phi1j",this,_Phi1j),
38 mZZj("mZZj","mZZj",this,_mZZj)
39 {
40 }
41 
42 
43 RooqqZZ_JHU::RooqqZZ_JHU(const RooqqZZ_JHU& other, const char* name) :
44 RooAbsPdf(other,name),
45 m1j("m1j",this,other.m1j),
46 m2j("m2j",this,other.m2j),
47 h1j("h1j",this,other.h1j),
48 h2j("h2j",this,other.h2j),
49 phij("phij",this,other.phij),
50 hsj("hsj",this,other.hsj),
51 Phi1j("Phi1j",this,other.Phi1j),
52 mZZj("mZZj",this,other.mZZj)
53 {
54 }
55 
56 
57 
58 Double_t RooqqZZ_JHU::evaluate() const
59 {
60  // Do transformation of the angles from JHU to NW conventions...
61 
62 
63  double m1_c = m1j;
64  double m2_c = m2j;
65  double h1_c = h1j;
66  double h2_c = h2j;
67  double phi_c = phij;
68  double hs_c = hsj;
69  double Phi1_c = Phi1j;
70  double mZZ_c = mZZj;
71 
72  // ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE
73  double up_plusZ = partonicXS(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
74  double up_minusZ = partonicXS(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
75  double down_plusZ = partonicXS(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
76  double down_minusZ = partonicXS(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
77 
78  //std::cout << "up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
79 
80  double upfrac = 0.5;
81  double downfrac = 0.5;
82 
83  //up_minusZ = 0; down_minusZ = 0;
84 
85  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
86 
87  double jacobian = 1./(mZZj*mZZj);
88  return jacobian*totalxs ;
89 }
90 
91 
92 Int_t RooqqZZ_JHU::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
93 {
94 
95  // these are the variables you integrated out!
96  if (matchArgs(allVars,analVars,h1j,h2j,phij,Phi1j)) return 5 ;
97  if (matchArgs(allVars,analVars,h2j,phij,Phi1j)) return 1 ;
98  if (matchArgs(allVars,analVars,h1j,phij,Phi1j)) return 2 ;
99  if (matchArgs(allVars,analVars,h1j,h2j,phij)) return 3 ;
100  if (matchArgs(allVars,analVars,h1j,h2j,Phi1j)) return 4 ;
101  return 0 ;
102 
103 }
104 Double_t RooqqZZ_JHU::analyticalIntegral(Int_t code, const char* rangeName) const
105 {
106 
107  double m1_c = m1j;
108  double m2_c = m2j;
109  double h1_c = h1j;
110  double h2_c = h2j;
111  double phi_c = phij;
112  double hs_c = hsj;
113  double Phi1_c = Phi1j;
114  double mZZ_c = mZZj;
115 
116 
117  switch(code)
118  {
119 
120  case 1:
121  {
122  double up_plusZ = integratePartonicXS_Dh2DphDPh1_1(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
123  double up_minusZ = integratePartonicXS_Dh2DphDPh1_1(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
124  double down_plusZ = integratePartonicXS_Dh2DphDPh1_1(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
125  double down_minusZ = integratePartonicXS_Dh2DphDPh1_1(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
126 
127  double upfrac = 0.5;
128  double downfrac = 0.5;
129 
130  //up_minusZ = 0; down_minusZ = 0;
131 
132  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
133  //std::cout << "case 1: up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
134 
135  double jacobian = 1./(mZZj*mZZj);
136  return jacobian*totalxs ;
137  }
138  case 2:
139  {
140  double up_plusZ = integratePartonicXS_Dh1DphDPh1_2(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
141  double up_minusZ = integratePartonicXS_Dh1DphDPh1_2(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
142  double down_plusZ = integratePartonicXS_Dh1DphDPh1_2(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
143  double down_minusZ = integratePartonicXS_Dh1DphDPh1_2(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
144 
145  double upfrac = 0.5;
146  double downfrac = 0.5;
147 
148  //up_minusZ = 0; down_minusZ = 0;
149 
150  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
151  //std::cout << "case 2: up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
152 
153  double jacobian = 1./(mZZj*mZZj);
154  return jacobian*totalxs ;
155  }
156  case 3:
157  {
158  double up_plusZ = integratePartonicXS_Dh1Dh2Dph_3(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
159  double up_minusZ = integratePartonicXS_Dh1Dh2Dph_3(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
160  double down_plusZ = integratePartonicXS_Dh1Dh2Dph_3(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
161  double down_minusZ = integratePartonicXS_Dh1Dh2Dph_3(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
162 
163  double upfrac = 0.5;
164  double downfrac = 0.5;
165 
166  //up_minusZ = 0; down_minusZ = 0;
167 
168  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
169  //std::cout << "case 3: up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
170 
171  double jacobian = 1./(mZZj*mZZj);
172  return jacobian*totalxs ;
173  }
174  case 4:
175  {
176  double up_plusZ = integratePartonicXS_Dh1Dh2DPh1_4(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
177  double up_minusZ = integratePartonicXS_Dh1Dh2DPh1_4(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
178  double down_plusZ = integratePartonicXS_Dh1Dh2DPh1_4(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
179  double down_minusZ = integratePartonicXS_Dh1Dh2DPh1_4(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
180 
181  double upfrac = 0.5;
182  double downfrac = 0.5;
183 
184  //up_minusZ = 0; down_minusZ = 0;
185 
186  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
187  //std::cout << "case 4: up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
188 
189  double jacobian = 1./(mZZj*mZZj);
190  return jacobian*totalxs ;
191  }
192  case 5:
193  {
194  double up_plusZ = integratePartonicXS_Dh1Dh2DphDPh1_5(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, true);
195  double up_minusZ = integratePartonicXS_Dh1Dh2DphDPh1_5(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, true, false);
196  double down_plusZ = integratePartonicXS_Dh1Dh2DphDPh1_5(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, true);
197  double down_minusZ = integratePartonicXS_Dh1Dh2DphDPh1_5(m1_c, m2_c, h1_c, h2_c, phi_c, hs_c, Phi1_c, mZZ_c, false, false);
198 
199  double upfrac = 0.5;
200  double downfrac = 0.5;
201 
202  //up_minusZ = 0; down_minusZ = 0;
203 
204  double totalxs = upfrac*(up_plusZ + up_minusZ) + downfrac*(down_plusZ + down_minusZ);
205  //std::cout << "case 5: up_plusZ = " << up_plusZ << ", up_minusZ = " << up_minusZ << ", down_plusZ = " << down_plusZ << ", down_minusZ = " << down_minusZ << std::endl;
206 
207  double jacobian = 1./(mZZj*mZZj);
208  return jacobian*totalxs ;
209  }
210 
211  }
212  assert(0) ;
213  return 0 ;
214 }
215 //*/
216 // 45 terms
217 Double_t RooqqZZ_JHU::partonicXS(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
218 {
219 
220  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
221  double s = mZZ_c*mZZ_c;
222  double M1 = m1_c;
223  double M2 = m2_c;
224  double hs = hs_c;
225  double h1 = h1_c;
226  double h2 = h2_c;
227  double phi = phi_c;
228  double Phi1 = Phi1_c;
229 
230  if (flipAxis){
231  hs *= -1.;
232  Phi1 += TMath::Pi();
233  //phi += 0;
234  }
235 
236  // extra definitions needed
237  double shs = Sqrt(1-hs*hs); // sin Theta
238  double sh1 = Sqrt(1-h1*h1); // sin theta1
239  double sh2 = Sqrt(1-h2*h2); // sin theta2
240  double hsdb = 1 - 2*shs*shs; // cos 2*Theta
241  double shsdb = 2*hs*shs; // sin 2*Theta
242  double h1db = 1 - 2*sh1*sh1; // cos 2*theta1
243  double h2db = 1 - 2*sh2*sh2; // cos 2*theta2
244 
245  double x = (M1*M1-M2*M2)/s;
246  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
247  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
248  double MZ = 91.18; // mass of the Z
249  double gamma = 2.475; // width of the Z
250  double xW = 0.2312; // sin2thetaW
251  double Pi = TMath::Pi();
252  double alpha = 1./137; // EM constant
253 
254  double Tf3 = -1/2.; double Qf = -1.;
255  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
256 
257  double qTf3, qQf;
258  if (upType){
259  qTf3 = 1/2.; qQf = 2./3.;
260  }
261  else{
262  qTf3 = -1/2.; qQf = -1./3.;
263  }
264  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
265 
266  // ----
267  const int nterms = 45;
268  double term[nterms];
269 
270  term[0] = (Power(alpha,4)*beta1*Power(M1,2)*Power(M2,2)*(1 + x)*(Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) + Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2))/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
271 
272  term[1] = (Power(alpha,4)*beta1*Power(M1,2)*Power(M2,2)*(1 + x)*(Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) + Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*(4*(-Power(gL,2) + Power(gR,2))*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(4*(-Power(gL,2) + Power(gR,2))*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2))/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
273 
274  term[2] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(M1,2)*Power(M2,2)*(1 + x)*(4*(-Power(gL,2) + Power(gR,2))*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2)*(Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
275 
276  term[3] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(M1,2)*Power(M2,2)*(1 + x)*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(4*(-Power(gL,2) + Power(gR,2))*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2)*(Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
277 
278  term[4] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*shs*(Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) - Power(gqL,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*sh1)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
279 
280  term[5] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(4*(gL - gR)*(gL + gR)*h2 - (Power(gL,2) + Power(gR,2))*(3 + h2db))*shs*(-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) + Power(gqR,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*sh1)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
281 
282  term[6] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*shs*(-(Power(gqR,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) + Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,4))*sh1)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - x)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
283 
284  term[7] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(4*(gL - gR)*(gL + gR)*h2 - (Power(gL,2) + Power(gR,2))*(3 + h2db))*shs*(Power(gqL,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) - Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,4))*sh1)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - x)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
285 
286  term[8] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(-Phi1 + Pi))*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2)*(Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs) + Power(gqL,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*Power(shs,2))*Power(sh1,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
287 
288  term[9] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(-Phi1 + Pi))*(4*(gL - gR)*(gL + gR)*h2 - (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2)*(-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)) - Power(gqR,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*Power(shs,2))*Power(sh1,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
289 
290  term[10] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*(Power(gqR,4)*Power(1 + hs,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) + Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*Power(shs,4))*Power(sh1,2))/(3.145728e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + x,2)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
291 
292  term[11] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*(4*(-Power(gL,2) + Power(gR,2))*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*(Power(gqL,4)*Power(1 + hs,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) + Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*Power(shs,4))*Power(sh1,2))/(3.145728e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + x,2)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
293 
294  term[12] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1)*(4*(gL - gR)*(gL + gR)*h1 - (Power(gL,2) + Power(gR,2))*(3 + h1db))*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) - Power(gqL,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
295 
296  term[13] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1)*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) + Power(gqR,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
297 
298  term[14] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1)*(4*(gL - gR)*(gL + gR)*h1 - (Power(gL,2) + Power(gR,2))*(3 + h1db))*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(Power(gqR,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) - Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
299 
300  term[15] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1)*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(-(Power(gqL,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) + Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
301 
302  term[16] = -(Power(alpha,4)*beta1*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi)*hs*((1 + beta1*beta2)*(Power(gqL,4) - Power(gqR,4))*x - (-1 + beta1*beta2)*(Power(gqL,4) + Power(gqR,4))*hs)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*sh1*sh2)/(49152.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
303 
304  term[17] = -(Power(alpha,4)*beta1*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi)*hs*((1 + beta1*beta2)*(-Power(gqL,4) + Power(gqR,4))*x - (-1 + beta1*beta2)*(Power(gqL,4) + Power(gqR,4))*hs)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*sh1*sh2)/(49152.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
305 
306  term[18] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi + 2*Phi1)*(Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs))*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*sh1*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
307 
308  term[19] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi + 2*Phi1)*(Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs))*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*sh1*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
309 
310  term[20] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi + 2*Phi1)*hs*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*(Power(gqR,4)*Power(1 + hs,2) - Power(gqL,4)*Power(shs,2))*sh1*sh2)/(49152.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
311 
312  term[21] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi + 2*Phi1)*hs*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,2)*(Power(gqL,4)*Power(1 + hs,2) - Power(gqR,4)*Power(shs,2))*sh1*sh2)/(49152.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
313 
314  term[22] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*(Power(gqR,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))*sh1*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
315 
316  term[23] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*Sqrt(1 - Power(beta2,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*(Power(gqL,4)*Power(1 + hs,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))*sh1*sh2)/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 - Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
317 
318  term[24] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1 - 2*(-Phi1 + Pi))*((-Power(gqL,4) + Power(gqR,4))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2)) + 4*(Power(gqL,4) + Power(gqR,4))*(1 + x)*hs)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,3)*Power(sh1,2)*sh2)/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
319 
320  term[25] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1 - 2*(-Phi1 + Pi))*((Power(gqL,4) - Power(gqR,4))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2)) + 4*(Power(gqL,4) + Power(gqR,4))*(1 + x)*hs)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*Power(shs,3)*Power(sh1,2)*sh2)/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
321 
322  term[26] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi + Phi1)*hs*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(-(Power(gqR,4)*Power(1 + hs,2)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) - Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*Power(sh1,2)*sh2)/(98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
323 
324  term[27] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi + Phi1)*hs*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(-(Power(gqL,4)*Power(1 + hs,2)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) - Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*Power(sh1,2)*sh2)/(98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
325 
326  term[28] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1 + 2*(-Phi1 + Pi))*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(Power(gqR,4)*Power(1 + hs,2)*(-((1 + beta1*beta2)*x) + (1 - beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqL,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*Power(sh1,2)*sh2)/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
327 
328  term[29] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*Cos(phi + Phi1 + 2*(-Phi1 + Pi))*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h2)*shs*(Power(gqL,4)*Power(1 + hs,2)*(-((1 + beta1*beta2)*x) + (1 - beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) + Power(gqR,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*Power(sh1,2)*sh2)/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
329 
330  term[30] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(phi + Phi1))*(4*(gL - gR)*(gL + gR)*h1 - (Power(gL,2) + Power(gR,2))*(3 + h1db))*Power(shs,2)*(Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs) + Power(gqL,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*Power(shs,2))*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
331 
332  term[31] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(phi + Phi1))*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*Power(shs,2)*(-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs)) - Power(gqR,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*Power(shs,2))*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
333 
334  term[32] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(4*(-Power(gL,2) + Power(gR,2))*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(Power(gqR,4)*Power(1 + hs,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) + Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*Power(shs,4))*Power(sh2,2))/(3.145728e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 + x)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
335 
336  term[33] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*(Power(gqL,4)*Power(1 + hs,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) + Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*Power(shs,4))*Power(sh2,2))/(3.145728e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 + x)*Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
337 
338  term[34] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1 + 2*(phi + Phi1))*((-Power(gqL,4) + Power(gqR,4))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2)) + 4*(Power(gqL,4) + Power(gqR,4))*(-1 + x)*hs)*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*Power(shs,3)*sh1*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
339 
340  term[35] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1 + 2*(phi + Phi1))*((Power(gqL,4) - Power(gqR,4))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2)) + 4*(Power(gqL,4) + Power(gqR,4))*(-1 + x)*hs)*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*Power(shs,3)*sh1*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
341 
342  term[36] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(Phi1)*hs*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*shs*(-(Power(gqR,4)*Power(1 + hs,2)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) - Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*sh1*Power(sh2,2))/(98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
343 
344  term[37] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(Phi1)*hs*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*shs*(-(Power(gqL,4)*Power(1 + hs,2)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) - Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))*sh1*Power(sh2,2))/(98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + Power(x,2))*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
345 
346  term[38] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1 - 2*(phi + Phi1))*(-Power(gL,2) + Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*shs*(Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) - Power(gqL,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*sh1*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
347 
348  term[39] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1 - 2*(phi + Phi1))*(Power(gL,2) - Power(gR,2) + (Power(gL,2) + Power(gR,2))*h1)*shs*(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) - Power(gqR,4)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2))*sh1*Power(sh2,2))/(393216.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(-1 + x)*Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
349 
350  term[40] = -(Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(phi + Pi))*(x + beta1*beta2*x + hs - beta1*beta2*hs)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*Power(shs,2)*Power(sh1,2)*Power(sh2,2))/(196608.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
351 
352  term[41] = (Power(alpha,4)*beta1*(-1 + Power(beta2,2))*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*Cos(2*(-Phi1 + Pi))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2)*Power(sh1,2)*Power(sh2,2))/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*(1 + x)*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
353 
354  term[42] = (Power(alpha,4)*beta1*(-1 + Power(beta1,2))*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(phi + Phi1))*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x + (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*Power(shs,2)*Power(sh1,2)*Power(sh2,2))/(786432.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + x,2)*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
355 
356  term[43] = -(Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(2*(-phi - 2*Phi1 + Pi))*Power(shs,4)*Power(sh1,2)*Power(sh2,2))/(196608.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2));
357 
358  term[44] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*Power(shsdb,2)*Power(sh1,2)*Power(sh2,2))/(98304.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*(Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,3)*s*Power(-1 + Power(x,2),2)*Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2));
359 
360  double returnVal = 0.;
361  for (int i = 0; i < nterms; i++){
362  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
363  if (term[i] != term[i]) {
364  std::cout << "term[" << i << "] isNaN" << std::endl;
365  term[i] = 0.;
366  }
367  returnVal += term[i];
368  }
369 
370  if (returnVal <= 0) returnVal = 1.e-40;
371 
372  //std::cout << "returnVal = " << returnVal << std::endl;
373 
374  return returnVal;
375 
376 }
377 
378 Double_t RooqqZZ_JHU::integratePartonicXS_Dh2DphDPh1_1(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
379 {
380 
381  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
382  double s = mZZ_c*mZZ_c;
383  double M1 = m1_c;
384  double M2 = m2_c;
385  double hs = hs_c;
386  double h1 = h1_c;
387  //double h2 = h2_c;
388  //double phi = phi_c; NOT USED!!
389  double Phi1 = Phi1_c;
390 
391  if (flipAxis){
392  hs *= -1.;
393  Phi1 += TMath::Pi();
394  //phi += 0;
395  }
396 
397  // extra definitions needed
398  double shs = Sqrt(1-hs*hs); // sin Theta
399  double sh1 = Sqrt(1-h1*h1); // sin theta1
400  //double sh2 = Sqrt(1-h2*h2); // sin theta2
401  //double hsdb = 1 - 2*shs*shs; // cos 2*Theta NOT USED!!
402  double shsdb = 2*hs*shs; // sin 2*Theta
403  double h1db = 1 - 2*sh1*sh1; // cos 2*theta1
404  //double h2db = 1 - 2*sh2*sh2; // cos 2*theta2 NOT USED!!
405 
406  double x = (M1*M1-M2*M2)/s;
407  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
408  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
409  double MZ = 91.18; // mass of the Z
410  double gamma = 2.475; // width of the Z
411  double xW = 0.2312; // sin2thetaW
412  double Pi = TMath::Pi();
413  double alpha = 1./137; // EM constant
414 
415  double Tf3 = -1/2.; double Qf = -1.;
416  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
417 
418  double qTf3, qQf;
419  if (upType){
420  qTf3 = 1/2.; qQf = 2./3.;
421  }
422  else{
423  qTf3 = -1/2.; qQf = -1./3.;
424  }
425  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
426 
427  // ----
428  const int nterms = 9;
429  double term[nterms];
430 
431 
432  term[0] = (Power(alpha,4)*beta1*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*
433  (Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
434  Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*
435  (4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*
436  Power(shs,2))/
437  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
438  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
439  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
440  4*beta1*beta2*Power(shs,2),2)) ;
441 
442  term[1] = (Power(alpha,4)*beta1*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*
443  (Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
444  Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*
445  (4*(-Power(gL,2) + Power(gR,2))*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*
446  Power(shs,2))/
447  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
448  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
449  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
450  4*beta1*beta2*Power(shs,2),2)) ;
451 
452  term[2] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
453  (1 + x)*(4*(-Power(gL,2) + Power(gR,2))*h1 +
454  (Power(gL,2) + Power(gR,2))*(3 + h1db))*Power(shs,2)*
455  (Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/
456  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
457  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
458  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
459  4*beta1*beta2*Power(shs,2),2)) ;
460 
461  term[3] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
462  (1 + x)*(4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*
463  Power(shs,2)*(Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/
464  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
465  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
466  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
467  4*beta1*beta2*Power(shs,2),2)) ;
468 
469  term[4] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
470  (4*(-Power(gL,2) + Power(gR,2))*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*
471  (Power(gqR,4)*Power(1 + hs,4)*
472  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
473  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
474  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
475  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
476  Power(shs,4)))/
477  (589824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
478  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
479  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
480  4*beta1*beta2*Power(shs,2),2)) ;
481 
482  term[5] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
483  (4*(gL - gR)*(gL + gR)*h1 + (Power(gL,2) + Power(gR,2))*(3 + h1db))*
484  (Power(gqL,4)*Power(1 + hs,4)*
485  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
486  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
487  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
488  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
489  Power(shs,4)))/
490  (589824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
491  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
492  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
493  4*beta1*beta2*Power(shs,2),2));
494 
495  term[6] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
496  (1 + x)*(Power(gqR,4)*Power(1 + hs,4)*
497  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
498  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
499  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
500  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
501  Power(shs,4))*Power(sh1,2))/
502  (147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
503  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
504  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
505  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
506 
507  term[7] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
508  Power(M2,2)*(1 + x)*(Power(gqL,4)*Power(1 + hs,4)*
509  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
510  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
511  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
512  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
513  Power(shs,4))*Power(sh1,2))/
514  (147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
515  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
516  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
517  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
518 
519  term[8] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*
520  Power(M1,2)*Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*
521  Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*
522  Power(shsdb,2)*Power(sh1,2))/
523  (18432.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*
524  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
525  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + Power(x,2),2)*
526  Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
527  4*beta1*beta2*Power(shs,2),2));
528 
529 
530  double returnVal = 0.;
531  for (int i = 0; i < nterms; i++){
532  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
533  if (term[i] != term[i]) {
534  //std::cout << "isNaN" << std::endl;
535  term[i] = 0.;
536  }
537  returnVal += term[i];
538  }
539 
540  if (returnVal <= 0) returnVal = 1.e-40;
541  return returnVal;
542 
543 }
544 
545 Double_t RooqqZZ_JHU::integratePartonicXS_Dh1DphDPh1_2(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
546 {
547  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
548  double s = mZZ_c*mZZ_c;
549  double M1 = m1_c;
550  double M2 = m2_c;
551  double hs = hs_c;
552  //double h1 = h1_c;
553  double h2 = h2_c;
554  //double phi = phi_c; NOT USED !!
555  double Phi1 = Phi1_c;
556 
557  if (flipAxis){
558  hs *= -1.;
559  Phi1 += TMath::Pi();
560  //phi += 0;
561  }
562 
563  // extra definitions needed
564  double shs = Sqrt(1-hs*hs); // sin Theta
565  //double sh1 = Sqrt(1-h1*h1); // sin theta1
566  double sh2 = Sqrt(1-h2*h2); // sin theta2
567  //double hsdb = 1 - 2*shs*shs; // cos 2*Theta NOT USED
568  double shsdb = 2*hs*shs; // sin 2*Theta
569  //double h1db = 1 - 2*sh1*sh1; // cos 2*theta1 NOT USED
570  double h2db = 1 - 2*sh2*sh2; // cos 2*theta2
571 
572  double x = (M1*M1-M2*M2)/s;
573  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
574  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
575  double MZ = 91.18; // mass of the Z
576  double gamma = 2.475; // width of the Z
577  double xW = 0.2312; // sin2thetaW
578  double Pi = TMath::Pi();
579  double alpha = 1./137; // EM constant
580 
581  double Tf3 = -1/2.; double Qf = -1.;
582  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
583 
584  double qTf3, qQf;
585  if (upType){
586  qTf3 = 1/2.; qQf = 2./3.;
587  }
588  else{
589  qTf3 = -1/2.; qQf = -1./3.;
590  }
591  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
592 
593  const int nterms = 9;
594  double term[nterms];
595 
596  term[0] = (Power(alpha,4)*beta1*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*
597  (Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
598  Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*
599  (4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*
600  Power(shs,2))/
601  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
602  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
603  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
604  4*beta1*beta2*Power(shs,2),2)) ;
605 
606  term[1] = (Power(alpha,4)*beta1*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*
607  (Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
608  Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*
609  (4*(-Power(gL,2) + Power(gR,2))*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*
610  Power(shs,2))/
611  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
612  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
613  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
614  4*beta1*beta2*Power(shs,2),2)) ;
615 
616  term[2] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
617  (1 + x)*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*
618  Power(shs,2)*(Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/
619  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
620  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
621  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
622  4*beta1*beta2*Power(shs,2),2)) ;
623 
624  term[3] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
625  (1 + x)*(4*(-Power(gL,2) + Power(gR,2))*h2 +
626  (Power(gL,2) + Power(gR,2))*(3 + h2db))*Power(shs,2)*
627  (Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/
628  (73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
629  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
630  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
631  4*beta1*beta2*Power(shs,2),2)) ;
632 
633  term[4] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
634  (1 + x)*(4*(gL - gR)*(gL + gR)*h2 + (Power(gL,2) + Power(gR,2))*(3 + h2db))*
635  (Power(gqR,4)*Power(1 + hs,4)*
636  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
637  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
638  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
639  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
640  Power(shs,4)))/
641  (589824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
642  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
643  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
644  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
645 
646  term[5] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
647  (1 + x)*(4*(-Power(gL,2) + Power(gR,2))*h2 +
648  (Power(gL,2) + Power(gR,2))*(3 + h2db))*
649  (Power(gqL,4)*Power(1 + hs,4)*
650  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
651  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
652  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
653  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
654  Power(shs,4)))/
655  (589824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
656  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
657  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
658  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
659 
660  term[6] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
661  Power(M2,2)*(Power(gqR,4)*Power(1 + hs,4)*
662  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
663  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
664  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
665  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
666  Power(shs,4))*Power(sh2,2))/
667  (147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
668  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
669  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
670  4*beta1*beta2*Power(shs,2),2)) ;
671 
672  term[7] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
673  (Power(gqL,4)*Power(1 + hs,4)*
674  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
675  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
676  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
677  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
678  Power(shs,4))*Power(sh2,2))/
679  (147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
680  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
681  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
682  4*beta1*beta2*Power(shs,2),2)) ;
683 
684  term[8] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
685  Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) +
686  (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*Power(shsdb,2)*
687  Power(sh2,2))/
688  (18432.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*
689  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
690  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + Power(x,2),2)*
691  Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
692  4*beta1*beta2*Power(shs,2),2));
693 
694  double returnVal = 0.;
695  for (int i = 0; i < nterms; i++){
696  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
697  if (term[i] != term[i]) {
698  //std::cout << "isNaN" << std::endl;
699  term[i] = 0.;
700  }
701  returnVal += term[i];
702  }
703 
704  if (returnVal <= 0) returnVal = 1.e-40;
705  return returnVal;
706 
707 }
708 
709 Double_t RooqqZZ_JHU::integratePartonicXS_Dh1Dh2Dph_3(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
710 {
711  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
712  double s = mZZ_c*mZZ_c;
713  double M1 = m1_c;
714  double M2 = m2_c;
715  double hs = hs_c;
716  // double h1 = h1_c;
717  // double h2 = h2_c;
718  //double phi = phi_c; NOT USED !!
719  double Phi1 = Phi1_c;
720 
721  if (flipAxis){
722  hs *= -1.;
723  Phi1 += TMath::Pi();
724  //phi += 0;
725  }
726 
727  /*
728  std::cout << "-----------" << std::endl;
729  std::cout << "h1 = " << h1 << "; " << std::endl;
730  std::cout << "h2 = " << h2 << "; " << std::endl;
731  std::cout << "hs = " << hs << "; " << std::endl;
732  std::cout << "phi = " << phi << "; " << std::endl;
733  std::cout << "Phi1 = " << Phi1 << "; " << std::endl;
734  //
735  */
736 
737  // extra definitions needed
738  double shs = Sqrt(1-hs*hs); // sin Theta
739  //double sh1 = Sqrt(1-h1*h1); // sin theta1
740  //double sh2 = Sqrt(1-h2*h2); // sin theta2
741  double hsdb = 1 - 2*shs*shs; // cos 2*Theta
742  double shsdb = 2*hs*shs; // sin 2*Theta
743  //double h1db = 1 - 2*sh1*sh1; // cos 2*theta1 NOT USED
744  //double h2db = 1 - 2*sh2*sh2; // cos 2*theta2 NOT USED!!
745 
746  double x = (M1*M1-M2*M2)/s;
747  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
748  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
749  double MZ = 91.18; // mass of the Z
750  double gamma = 2.475; // width of the Z
751  double xW = 0.2312; // sin2thetaW
752  double Pi = TMath::Pi();
753  double alpha = 1./137; // EM constant
754 
755  double Tf3 = -1/2.; double Qf = -1.;
756  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
757 
758  double qTf3, qQf;
759  if (upType){
760  qTf3 = 1/2.; qQf = 2./3.;
761  }
762  else{
763  qTf3 = -1/2.; qQf = -1./3.;
764  }
765  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
766 
767  /*
768  std::cout << "-----------" << std::endl;
769  std::cout << "beta1 = " << beta1 << "; " << std::endl;
770  std::cout << "beta2 = " << beta2 << "; " << std::endl;
771  std::cout << "gL = " << gL << "; " << std::endl;
772  std::cout << "gR = " << gR << "; " << std::endl;
773  std::cout << "gqL = " << gqL << "; " << std::endl;
774  std::cout << "gqR = " << gqR << "; " << std::endl;
775  //
776  */
777 
778  const int nterms = 18;
779  double term[nterms];
780 
782  term[0] = (Power(alpha,4)*beta1*(-1 + Power(beta2,2))*(Power(gqL,4) + Power(gqR,4))*
783  Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*Cos(2*(-Phi1 + Pi))*
784  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
785  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)*
786  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
787  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))/
788  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
789  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*(1 + x)*
790  Power(-1 + xW,4)*Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) +
791  Power(x + beta1*beta2*x,2) + 2*beta1*beta2*hsdb,2)) ;
792 
793 
794  term[1] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
795  (Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
796  Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
797  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
798  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
799  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
800  4*beta1*beta2*Power(shs,2),2)) ;
801 
802  term[2] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
803  (Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
804  Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
805  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
806  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
807  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
808  4*beta1*beta2*Power(shs,2),2)) ;
809 
810  term[3] = (Power(alpha,4)*beta1*(1 + beta1*beta2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
811  (1 + x)*Cos(2*(-Phi1 + Pi))*Power(shs,2)*
812  (Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs) +
813  Power(gqL,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*Power(shs,2)))/
814  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
815  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
816  Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
817  4*beta1*beta2*Power(shs,2),2)) ;
818 
819  term[4] = -(Power(alpha,4)*beta1*(1 + beta1*beta2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
820  (1 + x)*Cos(2*(-Phi1 + Pi))*Power(shs,2)*
821  (-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)) -
822  Power(gqR,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*Power(shs,2)))/
823  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
824  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
825  Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
826  4*beta1*beta2*Power(shs,2),2)) ;
827 
828  term[5] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(gL - gR)*(gL + gR)*(Power(gL,2) + Power(gR,2))*
829  Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*shs*
830  (Power(gqR,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*
831  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
832  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) -
833  Power(gqL,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*
834  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
835  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2)
836  ))/(147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
837  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(-1 + x)*Power(-1 + xW,4)*
838  Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
839  4*beta1*beta2*Power(shs,2),2)) ;
840 
841  term[6] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(-Power(gL,2) + Power(gR,2))*
842  (Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*shs*
843  (-(Power(gqL,4)*Power(1 + hs,2)*(x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*
844  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
845  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) +
846  Power(gqR,4)*(x + beta1*beta2*x + hs - beta1*beta2*hs)*
847  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
848  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,2)
849  ))/(147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
850  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(-1 + x)*Power(-1 + xW,4)*
851  Power(xW,4)*(1 + hs)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
852  4*beta1*beta2*Power(shs,2),2)) ;
853 
854  term[7] = (Power(alpha,4)*beta1*Sqrt((1 - Power(beta2,2))/((-1 + Power(beta1,2))*(-1 + Power(beta2,2))))*
855  (-Power(gL,2) + Power(gR,2))*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
856  ((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) +
857  (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(Phi1)*hs*shs*
858  (-(Power(gqR,4)*Power(1 + hs,2)*
859  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
860  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) -
861  Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
862  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))
863  )/(73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
864  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(-1 + Power(x,2))*
865  Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*
866  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
867 
868  term[8] = (Power(alpha,4)*beta1*Sqrt((1 - Power(beta2,2))/((-1 + Power(beta1,2))*(-1 + Power(beta2,2))))*
869  (gL - gR)*(gL + gR)*(Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*
870  ((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) +
871  (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(Phi1)*hs*shs*
872  (-(Power(gqL,4)*Power(1 + hs,2)*
873  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
874  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs)) -
875  Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
876  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,2))
877  )/(73728.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
878  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(-1 + Power(x,2))*
879  Power(-1 + xW,4)*Power(xW,4)*(1 + hs)*
880  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
881 
882  term[9] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
883  Power(M2,2)*(1 + x)*Power(shs,2)*
884  (Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/
885  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
886  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
887  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
888  4*beta1*beta2*Power(shs,2),2)) ;
889 
890  term[10] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
891  Power(M2,2)*(1 + x)*Power(shs,2)*
892  (Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/
893  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
894  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
895  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
896  4*beta1*beta2*Power(shs,2),2)) ;
897 
898  term[11] = -(Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*(-Power(gL,2) + Power(gR,2))*
899  (Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*shs*
900  (-(Power(gqR,4)*Power(1 + hs,4)*
901  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
902  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)) +
903  Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
904  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,4)
905  ))/(147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
906  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 - x)*Power(-1 + xW,4)*
907  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
908  4*beta1*beta2*Power(shs,2),2)) ;
909 
910  term[12] = (Power(alpha,4)*beta1*Sqrt(1 - Power(beta1,2))*(1 + beta1*beta2)*(gL - gR)*(gL + gR)*
911  (Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(Phi1)*shs*
912  (Power(gqL,4)*Power(1 + hs,4)*
913  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
914  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs) -
915  Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
916  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*Power(shs,4)
917  ))/(147456.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
918  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 - x)*Power(-1 + xW,4)*
919  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
920  4*beta1*beta2*Power(shs,2),2)) ;
921 
922  term[13] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
923  (1 + x)*(Power(gqR,4)*Power(1 + hs,4)*
924  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
925  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
926  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
927  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
928  Power(shs,4)))/
929  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
930  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + x,2)*
931  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
932  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
933  ) ;
934 
935  term[14] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
936  Power(M2,2)*(1 + x)*(Power(gqL,4)*Power(1 + hs,4)*
937  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
938  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
939  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
940  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
941  Power(shs,4)))/
942  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
943  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + x,2)*
944  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
945  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
946  ) ;
947 
948  term[15] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
949  Power(M2,2)*(Power(gqR,4)*Power(1 + hs,4)*
950  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
951  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
952  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
953  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
954  Power(shs,4)))/
955  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
956  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*(1 + x)*
957  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
958  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
959  ) ;
960 
961  term[16] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
962  Power(M2,2)*(Power(gqL,4)*Power(1 + hs,4)*
963  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
964  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
965  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
966  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
967  Power(shs,4)))/
968  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
969  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*(1 + x)*
970  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
971  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
972 
973  term[17] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*
974  Power(M1,2)*Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*
975  Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*
976  Power(shsdb,2))/
977  (27648.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*
978  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
979  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*
980  Power(-1 + Power(x,2),2)*Power(-1 + xW,4)*Power(xW,4)*
981  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)) ;
982 
983 
984  double returnVal = 0.;
985  //std::cout << "------" << std::endl;
986  for (int i = 0; i < nterms; i++){
987  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
988  if (term[i] != term[i]) {
989  //std::cout << "isNaN" << std::endl;
990  term[i] = 0.;
991  }
992  returnVal += term[i];
993  }
994 
995  if (returnVal <= 0) returnVal = 1.e-40;
996 
997  //std::cout << "returnVal = " << returnVal << std::endl;
998 
999  return returnVal;
1000 
1001 }
1002 Double_t RooqqZZ_JHU::integratePartonicXS_Dh1Dh2DPh1_4(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
1003 {
1004  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
1005  double s = mZZ_c*mZZ_c;
1006  double M1 = m1_c;
1007  double M2 = m2_c;
1008  double hs = hs_c;
1009  //double h1 = h1_c;
1010  //double h2 = h2_c;
1011  double phi = phi_c;
1012  double Phi1 = Phi1_c;
1013 
1014  if (flipAxis){
1015  hs *= -1.;
1016  Phi1 += TMath::Pi();
1017  //phi += 0;
1018  }
1019 
1020  // extra definitions needed
1021  double shs = Sqrt(1-hs*hs); // sin Theta
1022  //double sh1 = Sqrt(1-h1*h1); // sin theta1
1023  //double sh2 = Sqrt(1-h2*h2); // sin theta2
1024  double hsdb = 1 - 2*shs*shs; // cos 2*Theta
1025  double shsdb = 2*hs*shs; // sin 2*Theta
1026  //double h1db = 1 - 2*sh1*sh1; // cos 2*theta1 NOT USED
1027  //double h2db = 1 - 2*sh2*sh2; // cos 2*theta2 NOT USED
1028 
1029  double x = (M1*M1-M2*M2)/s;
1030  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
1031  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
1032  double MZ = 91.18; // mass of the Z
1033  double gamma = 2.475; // width of the Z
1034  double xW = 0.2312; // sin2thetaW
1035  double Pi = TMath::Pi();
1036  double alpha = 1./137; // EM constant
1037 
1038  double Tf3 = -1/2.; double Qf = -1.;
1039  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
1040 
1041  double qTf3, qQf;
1042  if (upType){
1043  qTf3 = 1/2.; qQf = 2./3.;
1044  }
1045  else{
1046  qTf3 = -1/2.; qQf = -1./3.;
1047  }
1048  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
1049 
1050  const int nterms = 14;
1051  double term[nterms];
1052 
1054  term[0] = -(Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1055  Power(M2,2)*(1 + x)*Cos(2*(phi + Pi))*(x + beta1*beta2*x + hs - beta1*beta2*hs)*
1056  (x + beta1*beta2*x + (-1 + beta1*beta2)*hs)*Power(shs,2))/
1057  (55296.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1058  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
1059  Power(xW,4)*Power(-1 - Power(beta1,2)*Power(beta2,2) + Power(x + beta1*beta2*x,2) +
1060  2*beta1*beta2*hsdb,2)) ;
1061 
1062  term[1] = - (Power(alpha,4)*beta1*Power(gL - gR,2)*Power(gL + gR,2)*
1063  Power(M1,2)*Power(M2,2)*(1 + x)*((beta1 + beta2 - beta1*x + beta2*x)*
1064  Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*
1065  Cos(phi)*hs*((1 + beta1*beta2)*(Power(gqL,4) - Power(gqR,4))*x -
1066  (-1 + beta1*beta2)*(Power(gqL,4) + Power(gqR,4))*hs)*Power(shs,2))/
1067  (98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*
1068  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1069  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*s*(1 - Power(x,2))*
1070  Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1071  4*beta1*beta2*Power(shs,2),2)) ;
1072 
1073  term[2] = - (Power(alpha,4)*beta1*Power(-Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
1074  ((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) +
1075  (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)))*Cos(phi)*hs*
1076  ((1 + beta1*beta2)*(-Power(gqL,4) + Power(gqR,4))*x -
1077  (-1 + beta1*beta2)*(Power(gqL,4) + Power(gqR,4))*hs)*Power(shs,2))/
1078  (98304.*Sqrt((-1 + Power(beta1,2))*(-1 + Power(beta2,2)))*
1079  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1080  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*s*(1 - Power(x,2))*
1081  Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1082  4*beta1*beta2*Power(shs,2),2)) ;
1083 
1084  term[3] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
1085  (Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
1086  Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
1087  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1088  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
1089  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1090  4*beta1*beta2*Power(shs,2),2)) ;
1091 
1092  term[4] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
1093  (Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
1094  Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
1095  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1096  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
1097  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1098  4*beta1*beta2*Power(shs,2),2)) ;
1099 
1100  term[5] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1101  Power(M2,2)*(1 + x)*Power(shs,2)*
1102  (Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/
1103  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1104  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
1105  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1106  4*beta1*beta2*Power(shs,2),2)) ;
1107 
1108  term[6] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1109  Power(M2,2)*(1 + x)*Power(shs,2)*
1110  (Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/
1111  (27648.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1112  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + xW,4)*
1113  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1114  4*beta1*beta2*Power(shs,2),2)) ;
1115 
1116  term[7] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
1117  (1 + x)*(Power(gqR,4)*Power(1 + hs,4)*
1118  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1119  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
1120  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1121  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
1122  Power(shs,4)))/
1123  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1124  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + x,2)*
1125  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1126  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1127  ) ;
1128 
1129  term[8] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1130  Power(M2,2)*(1 + x)*(Power(gqL,4)*Power(1 + hs,4)*
1131  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1132  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
1133  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1134  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
1135  Power(shs,4)))/
1136  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1137  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*Power(-1 + x,2)*
1138  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1139  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1140  ) ;
1141 
1142  term[9] = - (Power(alpha,4)*beta1*Sqrt((1 - Power(beta1,2))*(1 - Power(beta2,2)))*(gL - gR)*(gL + gR)*
1143  (-Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi)*
1144  (Power(gqR,4)*Power(1 + hs,4)*
1145  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1146  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*
1147  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1148  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) +
1149  Power(gqL,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1150  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*
1151  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1152  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))
1153  )/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1154  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*s*(1 - Power(x,2))*
1155  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1156  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1157  );
1158 
1159  term[10] = - (Power(alpha,4)*beta1*Sqrt((1 - Power(beta1,2))*(1 - Power(beta2,2)))*(gL - gR)*(gL + gR)*
1160  (-Power(gL,2) + Power(gR,2))*Power(M1,2)*Power(M2,2)*(1 + x)*Cos(phi)*
1161  (Power(gqL,4)*Power(1 + hs,4)*
1162  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1163  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs)*
1164  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1165  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs) +
1166  Power(gqR,4)*(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1167  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs)*
1168  (2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1169  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs)*Power(shs,4))
1170  )/(1.572864e6*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1171  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*s*(1 - Power(x,2))*
1172  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1173  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1174  ) ;
1175 
1176  term[11] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1177  Power(M2,2)*(Power(gqR,4)*Power(1 + hs,4)*
1178  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1179  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
1180  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1181  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
1182  Power(shs,4)))/
1183  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1184  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*(1 + x)*
1185  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1186  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1187  ) ;
1188 
1189  term[12] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1190  Power(M2,2)*(Power(gqL,4)*Power(1 + hs,4)*
1191  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1192  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
1193  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1194  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
1195  Power(shs,4)))/
1196  (221184.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1197  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*(1 + x)*
1198  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1199  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1200  ) ;
1201 
1202  term[13] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*
1203  Power(M1,2)*Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*
1204  Sqrt(beta1*beta2*(1 - Power(x,2))) + (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*
1205  Power(shsdb,2))/
1206  (27648.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*
1207  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1208  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Power(Pi,2)*s*
1209  Power(-1 + Power(x,2),2)*Power(-1 + xW,4)*Power(xW,4)*
1210  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1211  );
1212 
1213  double returnVal = 0.;
1214  for (int i = 0; i < nterms; i++){
1215  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
1216  if (term[i] != term[i]) {
1217  //std::cout << "isNaN" << std::endl;
1218  term[i] = 0.;
1219  }
1220  returnVal += term[i];
1221  }
1222 
1223  if (returnVal <= 0) returnVal = 1.e-40;
1224  return returnVal;
1225 
1226 }
1227 Double_t RooqqZZ_JHU::integratePartonicXS_Dh1Dh2DphDPh1_5(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
1228 {
1229  // define hs, h1, h2, phi1, phi2 for the NW angular conventions
1230  double s = mZZ_c*mZZ_c;
1231  double M1 = m1_c;
1232  double M2 = m2_c;
1233  double hs = hs_c;
1234  //double h1 = h1_c;
1235  //double h2 = h2_c;
1236  //double phi = phi_c; NOT USED
1237  double Phi1 = Phi1_c;
1238 
1239  if (flipAxis){
1240  hs *= -1.;
1241  Phi1 += TMath::Pi();
1242  //phi += 0;
1243  }
1244 
1245  // extra definitions needed
1246  double shs = Sqrt(1-hs*hs); // sin Theta
1247  //double sh1 = Sqrt(1-h1*h1); // sin theta1
1248  //double sh2 = Sqrt(1-h2*h2); // sin theta2
1249  //double hsdb = 1 - 2*shs*shs; // cos 2*Theta NOT USED!!
1250  double shsdb = 2*hs*shs; // sin 2*Theta
1251  //double h1db = 1 - 2*sh1*sh1; // cos 2*theta1 NOT USED!!
1252  //double h2db = 1 - 2*sh2*sh2; // cos 2*theta2 NOT USED!!
1253 
1254  double x = (M1*M1-M2*M2)/s;
1255  double beta1 = Sqrt(1. - (4.*M1*M1)/(s*(1+x)*(1+x)));
1256  double beta2 = Sqrt(1. - (4.*M2*M2)/(s*(1+x)*(1+x)));
1257  double MZ = 91.18; // mass of the Z
1258  double gamma = 2.475; // width of the Z
1259  double xW = 0.2312; // sin2thetaW
1260  double Pi = TMath::Pi();
1261  double alpha = 1./137; // EM constant
1262 
1263  double Tf3 = -1/2.; double Qf = -1.;
1264  double gL = 2.*(Tf3 - xW*Qf), gR = -2.*xW*Qf;
1265 
1266  double qTf3, qQf;
1267  if (upType){
1268  qTf3 = 1/2.; qQf = 2./3.;
1269  }
1270  else{
1271  qTf3 = -1/2.; qQf = -1./3.;
1272  }
1273  double gqL = 2.*(qTf3 - xW*qQf), gqR = -2.*xW*qQf;
1274 
1275 
1276  const int nterms = 9;
1277  double term[nterms];
1278 
1280  term[0] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
1281  (Power(gqR,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
1282  Power(gqL,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
1283  (13824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1284  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
1285  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1286  4*beta1*beta2*Power(shs,2),2)) ;
1287 
1288  term[1] = (Power(alpha,4)*beta1*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*(1 + x)*
1289  (Power(gqL,4)*Power(x + beta1*beta2*x + hs - beta1*beta2*hs,2) +
1290  Power(gqR,4)*Power(x + beta1*beta2*x + (-1 + beta1*beta2)*hs,2))*Power(shs,2))/
1291  (13824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1292  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
1293  Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1294  4*beta1*beta2*Power(shs,2),2)) ;
1295 
1296  term[2] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1297  Power(M2,2)*(1 + x)*Power(shs,2)*
1298  (Power(gqR,4)*Power(1 + hs,4) + Power(gqL,4)*Power(shs,4)))/
1299  (13824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1300  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
1301  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1302  4*beta1*beta2*Power(shs,2),2)) ;
1303 
1304  term[3] = (Power(alpha,4)*beta1*Power(1 + beta1*beta2,2)*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1305  Power(M2,2)*(1 + x)*Power(shs,2)*
1306  (Power(gqL,4)*Power(1 + hs,4) + Power(gqR,4)*Power(shs,4)))/
1307  (13824.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1308  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + xW,4)*
1309  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1310  4*beta1*beta2*Power(shs,2),2));
1311 
1312  term[4] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
1313  (1 + x)*(Power(gqR,4)*Power(1 + hs,4)*
1314  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1315  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
1316  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1317  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
1318  Power(shs,4)))/
1319  (110592.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1320  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
1321  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1322  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1323  ) ;
1324 
1325  term[5] = (Power(alpha,4)*beta1*(1 - Power(beta1,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1326  Power(M2,2)*(1 + x)*(Power(gqL,4)*Power(1 + hs,4)*
1327  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1328  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(-1 + x)*hs,2) +
1329  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1330  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(-1 + x)*hs,2)*
1331  Power(shs,4)))/
1332  (110592.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1333  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + x,2)*
1334  Power(-1 + xW,4)*Power(xW,4)*Power(1 + hs,2)*
1335  Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) + 4*beta1*beta2*Power(shs,2),2)
1336  ) ;
1337 
1338  term[6] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1339  Power(M2,2)*(Power(gqR,4)*Power(1 + hs,4)*
1340  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1341  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
1342  Power(gqL,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1343  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
1344  Power(shs,4)))/
1345  (110592.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1346  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
1347  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1348  4*beta1*beta2*Power(shs,2),2)) ;
1349 
1350  term[7] = (Power(alpha,4)*beta1*(1 - Power(beta2,2))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*Power(M2,2)*
1351  (Power(gqL,4)*Power(1 + hs,4)*
1352  Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1353  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) - 4*(1 + x)*hs,2) +
1354  Power(gqR,4)*Power(2 + Power(beta1,2) + Power(beta2,2) + 2*(Power(beta1,2) - Power(beta2,2))*x +
1355  (-2 + Power(beta1,2) + Power(beta2,2))*Power(x,2) + 4*(1 + x)*hs,2)*
1356  Power(shs,4)))/
1357  (110592.*(Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1358  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*(1 + x)*Power(-1 + xW,4)*
1359  Power(xW,4)*Power(1 + hs,2)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1360  4*beta1*beta2*Power(shs,2),2)) ;
1361 
1362  term[8] = (Power(alpha,4)*beta1*(Power(gqL,4) + Power(gqR,4))*Power(Power(gL,2) + Power(gR,2),2)*Power(M1,2)*
1363  Power(M2,2)*(1 + x)*Power((beta1 + beta2 - beta1*x + beta2*x)*Sqrt(beta1*beta2*(1 - Power(x,2))) +
1364  (1 + Power(beta1,2)*Power(beta2,2))*(-1 + Power(x,2)),2)*Power(shsdb,2))/
1365  (13824.*(-1 + Power(beta1,2))*(-1 + Power(beta2,2))*
1366  (Power(gamma,2)*Power(MZ,2) + Power(Power(M1,2) - Power(MZ,2),2))*
1367  (Power(gamma,2)*Power(MZ,2) + Power(Power(M2,2) - Power(MZ,2),2))*Pi*s*Power(-1 + Power(x,2),2)*
1368  Power(-1 + xW,4)*Power(xW,4)*Power(Power(-1 + beta1*beta2,2) - Power(x + beta1*beta2*x,2) +
1369  4*beta1*beta2*Power(shs,2),2));
1370 
1371  double returnVal = 0.;
1372  for (int i = 0; i < nterms; i++){
1373  //std::cout << "term[" << i << "] = " << term[i] << std::endl;
1374  if (term[i] != term[i]) {
1375  //std::cout << "isNaN" << std::endl;
1376  term[i] = 0.;
1377  }
1378  returnVal += term[i];
1379  }
1380 
1381  if (returnVal <= 0) returnVal = 1.e-40;
1382  return returnVal;
1383 
1384 }
RooqqZZ_JHU.h
RooqqZZ_JHU::m2j
RooRealProxy m2j
Definition: RooqqZZ_JHU.h:53
RooqqZZ_JHU::h2j
RooRealProxy h2j
Definition: RooqqZZ_JHU.h:55
RooqqZZ_JHU::integratePartonicXS_Dh1Dh2DphDPh1_5
Double_t integratePartonicXS_Dh1Dh2DphDPh1_5(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:1227
RooqqZZ_JHU
Definition: RooqqZZ_JHU.h:21
hto_betacom::beta2
real *8, dimension(3:6) beta2
Definition: CALLING_cpHTO.f:2080
RooqqZZ_JHU::integratePartonicXS_Dh1DphDPh1_2
Double_t integratePartonicXS_Dh1DphDPh1_2(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:545
RooqqZZ_JHU::Phi1j
RooRealProxy Phi1j
Definition: RooqqZZ_JHU.h:58
RooqqZZ_JHU::h1j
RooRealProxy h1j
Definition: RooqqZZ_JHU.h:54
RooqqZZ_JHU::RooqqZZ_JHU
RooqqZZ_JHU()
Definition: RooqqZZ_JHU.h:23
RooqqZZ_JHU::partonicXS
Double_t partonicXS(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:217
RooqqZZ_JHU::hsj
RooRealProxy hsj
Definition: RooqqZZ_JHU.h:57
RooqqZZ_JHU::analyticalIntegral
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Definition: RooqqZZ_JHU.cc:104
RooqqZZ_JHU::getAnalyticalIntegral
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Definition: RooqqZZ_JHU.cc:92
RooqqZZ_JHU::evaluate
Double_t evaluate() const
Definition: RooqqZZ_JHU.cc:58
RooqqZZ_JHU::integratePartonicXS_Dh2DphDPh1_1
Double_t integratePartonicXS_Dh2DphDPh1_1(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:378
RooqqZZ_JHU::integratePartonicXS_Dh1Dh2Dph_3
Double_t integratePartonicXS_Dh1Dh2Dph_3(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:709
RooqqZZ_JHU::phij
RooRealProxy phij
Definition: RooqqZZ_JHU.h:56
dd_global::cout
integer cout
Definition: DD_global.F90:21
RooqqZZ_JHU::m1j
RooRealProxy m1j
Definition: RooqqZZ_JHU.h:52
RooqqZZ_JHU::mZZj
RooRealProxy mZZj
Definition: RooqqZZ_JHU.h:59
hto_betacom::beta1
real *8, dimension(3:6) beta1
Definition: CALLING_cpHTO.f:2080
RooqqZZ_JHU::integratePartonicXS_Dh1Dh2DPh1_4
Double_t integratePartonicXS_Dh1Dh2DPh1_4(double m1_c, double m2_c, double h1_c, double h2_c, double phi_c, double hs_c, double Phi1_c, double mZZ_c, bool upType, bool flipAxis) const
Definition: RooqqZZ_JHU.cc:1002