15 #include "TLorentzRotation.h"
42 TLorentzVector
const nullFourVector(0, 0, 0, 0);
43 TLorentzVector p1hat, p2hat;
44 if (p1==nullFourVector || p2==nullFourVector)
return;
47 TLorentzVector p12=p1+p2;
48 TLorentzVector diffp2p1=p2-p1;
49 double p1sq = p1.M2();
50 double p2sq = p2.M2();
51 double p1p2 = p1.Dot(p2);
52 double m1sq = m1*fabs(m1);
53 double m2sq = m2*fabs(m2);
54 double p12sq = p12.M2();
56 TLorentzVector avec=(p1sq*p2 - p2sq*p1 + p1p2*diffp2p1);
58 double b = (p12sq + m2sq - m1sq) * (pow(p1p2, 2) - p1sq*p2sq);
59 double c = pow((p12sq + m2sq - m1sq), 2)*p1sq/4. - pow((p1sq + p1p2), 2)*m2sq;
60 double eta = (-b - sqrt(fabs(b*b -4.*a*c)))/(2.*a);
61 double xi = (p12sq + m2sq - m1sq - 2.*eta*(p2sq + p1p2))/(2.*(p1sq + p1p2));
63 p1hat = (1.-xi)*p1 + (1.-eta)*p2;
64 p2hat = xi*p1 + eta*p2;
69 double energy, p3, newp3, ratio;
70 energy = massiveJet.T();
72 newp3 = sqrt(max(pow(energy, 2)-mass*fabs(mass), 0.));
73 ratio = (p3>0. ? (newp3/p3) : 1.);
74 masslessJet.SetXYZT(massiveJet.X()*ratio, massiveJet.Y()*ratio, massiveJet.Z()*ratio, energy);
77 TLorentzVector
const& jet1,
int const& jet1Id,
78 TLorentzVector
const& jet2,
int const& jet2Id,
81 TLorentzVector
const nullFourVector(0, 0, 0, 0);
82 TLorentzVector jet1massless(0, 0, 0, 0), jet2massless(0, 0, 0, 0);
131 pair<TLorentzVector, TLorentzVector> result(jet1massless, jet2massless);
135 if (daughters.size()!=3)
return;
137 daughters.at(1).second, daughters.at(1).first,
138 daughters.at(2).second, daughters.at(2).first
140 daughters.at(1).second = corrWpair.first;
141 daughters.at(2).second = corrWpair.second;
142 TLorentzVector pW = daughters.at(1).second+daughters.at(2).second;
143 TVector3 pW_boost_old = -pW.BoostVector();
145 daughters.at(0).second, daughters.at(0).first,
146 pW, -daughters.at(0).first,
149 daughters.at(0).second=corrbW.first;
151 TVector3 pW_boost_new = pW.BoostVector();
152 for (
unsigned int idau=1; idau<daughters.size(); idau++){
153 daughters.at(idau).second.Boost(pW_boost_old);
154 daughters.at(idau).second.Boost(pW_boost_new);
159 TLorentzVector masslessRealJet(0, 0, 0, 0);
161 else masslessRealJet = realJet;
162 fakeJet = others + masslessRealJet;
163 fakeJet.SetVect(-fakeJet.Vect());
164 fakeJet.SetE(fakeJet.P());
168 std::pair<TLorentzVector, TLorentzVector>
TUtil::ComplexBoost(TVector3
const& beta, TLorentzVector
const& p4){
172 double b2 = bx*bx + by*by + bz*bz;
173 double gammasqinv = 1.-
b2;
176 double bp = bx*p4.X() + by*p4.Y() + bz*p4.Z();
177 double gammap_real=0;
178 double gammap_imag=0;
179 TLorentzVector p4new_real(0, 0, 0, 0), p4new_imag(0, 0, 0, 0);
181 gamma = 1./sqrt(gammasqinv);
182 if (
b2>0.) gammap_real = (gamma-1.)/
b2;
184 p4new_real.SetX(p4.X() + gammap_real*bp*bx + gamma*bx*p4.T());
185 p4new_real.SetY(p4.Y() + gammap_real*bp*by + gamma*by*p4.T());
186 p4new_real.SetZ(p4.Z() + gammap_real*bp*bz + gamma*bz*p4.T());
187 p4new_real.SetT(gamma*(p4.T() + bp));
189 else if (gammasqinv<0.){
190 gamma = -1./sqrt(-gammasqinv);
192 gammap_real = -1./
b2;
193 gammap_imag = gamma/
b2;
196 p4new_real.SetX(p4.X() + gammap_real*bp*bx);
197 p4new_real.SetY(p4.Y() + gammap_real*bp*by);
198 p4new_real.SetZ(p4.Z() + gammap_real*bp*bz);
200 p4new_imag.SetX(gammap_imag*bp*bx + gamma*bx*p4.T());
201 p4new_imag.SetY(gammap_imag*bp*by + gamma*by*p4.T());
202 p4new_imag.SetZ(gammap_imag*bp*bz + gamma*bz*p4.T());
203 p4new_imag.SetT(gamma*(p4.T() + bp));
206 return (pair<TLorentzVector, TLorentzVector>(p4new_real, p4new_imag));
216 TLorentzVector p4M11,
int Z1_lept1Id,
217 TLorentzVector p4M12,
int Z1_lept2Id,
218 TLorentzVector p4M21,
int Z2_lept1Id,
219 TLorentzVector p4M22,
int Z2_lept2Id
221 TLorentzVector
const nullFourVector(0, 0, 0, 0);
222 if (p4M12==nullFourVector || p4M22==nullFourVector){
224 p4M11 = f13Pair.first;
225 p4M21 = f13Pair.second;
227 else if (p4M11==nullFourVector || p4M21==nullFourVector){
229 p4M12 = f24Pair.first;
230 p4M22 = f24Pair.second;
235 p4M11 = f12Pair.first;
236 p4M12 = f12Pair.second;
237 p4M21 = f34Pair.first;
238 p4M22 = f34Pair.second;
242 TLorentzVector p4Z1 = p4M11 + p4M12;
243 TLorentzVector p4Z2 = p4M21 + p4M22;
250 (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0)
252 ((Z1_lept1Id*Z1_lept2Id>0 || (Z1_lept1Id==0 && Z1_lept2Id==0)) && p4M11.Phi()<=p4M12.Phi())
254 ) swap(p4M11, p4M12);
260 (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0)
262 ((Z2_lept1Id*Z2_lept2Id>0 || (Z2_lept1Id==0 && Z2_lept2Id==0)) && p4M21.Phi()<=p4M22.Phi())
264 ) swap(p4M21, p4M22);
269 TLorentzVector p4H = p4Z1 + p4Z2;
274 TVector3 boostX = -(p4H.BoostVector());
275 TLorentzVector thep4Z1inXFrame(p4Z1);
276 TLorentzVector thep4Z2inXFrame(p4Z2);
277 thep4Z1inXFrame.Boost(boostX);
278 thep4Z2inXFrame.Boost(boostX);
279 TVector3 theZ1X_p3 = thep4Z1inXFrame.Vect();
280 TVector3 theZ2X_p3 = thep4Z2inXFrame.Vect();
281 costhetastar = theZ1X_p3.CosTheta();
283 TVector3 boostV1(0, 0, 0);
284 TVector3 boostV2(0, 0, 0);
286 if (!(fabs(Z1_lept1Id)==21 || fabs(Z1_lept1Id)==22 || fabs(Z1_lept2Id)==21 || fabs(Z1_lept2Id)==22)){
287 boostV1 = -(p4Z1.BoostVector());
288 if (boostV1.Mag()>=1.) {
289 MELAout <<
"Warning: Mela::computeAngles: Z1 boost with beta=1, scaling down" << endl;
290 boostV1*=0.9999/boostV1.Mag();
292 TLorentzVector p4M11_BV1(p4M11);
293 TLorentzVector p4M12_BV1(p4M12);
294 TLorentzVector p4M21_BV1(p4M21);
295 TLorentzVector p4M22_BV1(p4M22);
296 p4M11_BV1.Boost(boostV1);
297 p4M12_BV1.Boost(boostV1);
298 p4M21_BV1.Boost(boostV1);
299 p4M22_BV1.Boost(boostV1);
301 TLorentzVector p4V2_BV1 = p4M21_BV1 + p4M22_BV1;
303 costheta1 = -p4V2_BV1.Vect().Unit().Dot(p4M11_BV1.Vect().Unit());
308 if (!(fabs(Z2_lept1Id)==21 || fabs(Z2_lept1Id)==22 || fabs(Z2_lept2Id)==21 || fabs(Z2_lept2Id)==22)){
309 boostV2 = -(p4Z2.BoostVector());
310 if (boostV2.Mag()>=1.) {
311 MELAout <<
"Warning: Mela::computeAngles: Z2 boost with beta=1, scaling down" << endl;
312 boostV2*=0.9999/boostV2.Mag();
314 TLorentzVector p4M11_BV2(p4M11);
315 TLorentzVector p4M12_BV2(p4M12);
316 TLorentzVector p4M21_BV2(p4M21);
317 TLorentzVector p4M22_BV2(p4M22);
318 p4M11_BV2.Boost(boostV2);
319 p4M12_BV2.Boost(boostV2);
320 p4M21_BV2.Boost(boostV2);
321 p4M22_BV2.Boost(boostV2);
323 TLorentzVector p4V1_BV2 = p4M11_BV2 + p4M12_BV2;
325 costheta2 = -p4V1_BV2.Vect().Unit().Dot(p4M21_BV2.Vect().Unit());
330 TLorentzVector p4M11_BX(p4M11);
331 TLorentzVector p4M12_BX(p4M12);
332 TLorentzVector p4M21_BX(p4M21);
333 TLorentzVector p4M22_BX(p4M22);
335 p4M11_BX.Boost(boostX);
336 p4M12_BX.Boost(boostX);
337 p4M21_BX.Boost(boostX);
338 p4M22_BX.Boost(boostX);
339 TLorentzVector p4V1_BX = p4M11_BX + p4M12_BX;
341 TVector3
const beamAxis(0, 0, 1);
342 TVector3 p3V1_BX = p4V1_BX.Vect().Unit();
343 TVector3 normal1_BX = (p4M11_BX.Vect().Cross(p4M12_BX.Vect())).Unit();
344 TVector3 normal2_BX = (p4M21_BX.Vect().Cross(p4M22_BX.Vect())).Unit();
345 TVector3 normalSC_BX = (beamAxis.Cross(p3V1_BX)).Unit();
349 float tmpSgnPhi = p3V1_BX.Dot(normal1_BX.Cross(normal2_BX));
351 if (fabs(tmpSgnPhi)>0.) sgnPhi = tmpSgnPhi/fabs(tmpSgnPhi);
352 float dot_BX12 = normal1_BX.Dot(normal2_BX);
353 if (fabs(dot_BX12)>=1.) dot_BX12 *= 1./fabs(dot_BX12);
354 Phi = sgnPhi * acos(-1.*dot_BX12);
358 float tmpSgnPhi1 = p3V1_BX.Dot(normal1_BX.Cross(normalSC_BX));
360 if (fabs(tmpSgnPhi1)>0.) sgnPhi1 = tmpSgnPhi1/fabs(tmpSgnPhi1);
361 float dot_BX1SC = normal1_BX.Dot(normalSC_BX);
362 if (fabs(dot_BX1SC)>=1.) dot_BX1SC *= 1./fabs(dot_BX1SC);
363 Phi1 = sgnPhi1 * acos(dot_BX1SC);
366 MELAout <<
"WARNING: NaN in computeAngles: "
367 << costhetastar <<
" "
371 << Phi1 <<
" " << endl;
372 MELAout <<
" boostV1: " <<boostV1.Pt() <<
" " << boostV1.Eta() <<
" " << boostV1.Phi() <<
" " << boostV1.Mag() << endl;
373 MELAout <<
" boostV2: " <<boostV2.Pt() <<
" " << boostV2.Eta() <<
" " << boostV2.Phi() <<
" " << boostV2.Mag() << endl;
383 TLorentzVector p4M11,
int Z1_lept1Id,
384 TLorentzVector p4M12,
int Z1_lept2Id,
385 TLorentzVector p4M21,
int Z2_lept1Id,
386 TLorentzVector p4M22,
int Z2_lept2Id
388 TLorentzVector
const nullFourVector(0, 0, 0, 0);
389 if (p4M12==nullFourVector || p4M22==nullFourVector){
391 p4M11 = f13Pair.first;
392 p4M21 = f13Pair.second;
394 else if (p4M11==nullFourVector || p4M21==nullFourVector){
396 p4M12 = f24Pair.first;
397 p4M22 = f24Pair.second;
402 p4M11 = f12Pair.first;
403 p4M12 = f12Pair.second;
404 p4M21 = f34Pair.first;
405 p4M22 = f34Pair.second;
408 TVector3 LabXaxis(1.0, 0.0, 0.0);
409 TVector3 LabYaxis(0.0, 1.0, 0.0);
410 TVector3 LabZaxis(0.0, 0.0, 1.0);
413 float Ebeam = sqrt(pbeam*pbeam + Mprot*Mprot);
415 TLorentzVector targ(0., 0., -pbeam, Ebeam);
416 TLorentzVector beam(0., 0., pbeam, Ebeam);
419 TLorentzVector p4Z1 = p4M11 + p4M12;
420 TLorentzVector p4Z2 = p4M21 + p4M22;
427 (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0)
429 ((Z1_lept1Id*Z1_lept2Id>0 || (Z1_lept1Id==0 && Z1_lept2Id==0)) && p4M11.Phi()<=p4M12.Phi())
431 ) swap(p4M11, p4M12);
437 (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0)
439 ((Z2_lept1Id*Z2_lept2Id>0 || (Z2_lept1Id==0 && Z2_lept2Id==0)) && p4M21.Phi()<=p4M22.Phi())
441 ) swap(p4M21, p4M22);
445 TLorentzVector p4H = p4Z1 + p4Z2;
446 TVector3 boostX = -(p4H.BoostVector());
455 TRotation rotationCS;
457 TLorentzVector beaminX(beam);
458 TLorentzVector targinX(targ);
459 targinX.Boost(boostX);
460 beaminX.Boost(boostX);
463 TVector3 beam_targ_bisecinX((beaminX.Vect().Unit() - targinX.Vect().Unit()).Unit());
466 TVector3 newZaxisCS(beam_targ_bisecinX.Unit());
467 TVector3 newYaxisCS(beaminX.Vect().Unit().Cross(newZaxisCS).Unit());
468 TVector3 newXaxisCS(newYaxisCS.Unit().Cross(newZaxisCS).Unit());
469 rotationCS.RotateAxes(newXaxisCS, newYaxisCS, newZaxisCS);
473 TLorentzVector thep4Z1inXFrame_rotCS(p4Z1);
474 TLorentzVector thep4Z2inXFrame_rotCS(p4Z2);
475 thep4Z1inXFrame_rotCS.Transform(rotationCS);
476 thep4Z2inXFrame_rotCS.Transform(rotationCS);
477 thep4Z1inXFrame_rotCS.Boost(boostX);
478 thep4Z2inXFrame_rotCS.Boost(boostX);
479 TVector3 theZ1XrotCS_p3 = TVector3(thep4Z1inXFrame_rotCS.X(), thep4Z1inXFrame_rotCS.Y(), thep4Z1inXFrame_rotCS.Z());
480 costhetastar = theZ1XrotCS_p3.CosTheta();
484 TLorentzVector p4M11_BX_rotCS(p4M11);
485 TLorentzVector p4M12_BX_rotCS(p4M12);
486 TLorentzVector p4M21_BX_rotCS(p4M21);
487 TLorentzVector p4M22_BX_rotCS(p4M22);
488 p4M11_BX_rotCS.Transform(rotationCS);
489 p4M12_BX_rotCS.Transform(rotationCS);
490 p4M21_BX_rotCS.Transform(rotationCS);
491 p4M22_BX_rotCS.Transform(rotationCS);
492 p4M11_BX_rotCS.Boost(boostX);
493 p4M12_BX_rotCS.Boost(boostX);
494 p4M21_BX_rotCS.Boost(boostX);
495 p4M22_BX_rotCS.Boost(boostX);
497 TLorentzVector p4Z1_BX_rotCS = p4M11_BX_rotCS + p4M12_BX_rotCS;
498 TVector3 p3V1_BX_rotCS = (p4Z1_BX_rotCS.Vect()).Unit();
499 TVector3
const beamAxis(0, 0, 1);
500 TVector3 normal1_BX_rotCS = (p4M11_BX_rotCS.Vect().Cross(p4M12_BX_rotCS.Vect())).Unit();
501 TVector3 normal2_BX_rotCS = (p4M21_BX_rotCS.Vect().Cross(p4M22_BX_rotCS.Vect())).Unit();
502 TVector3 normalSC_BX_rotCS = (beamAxis.Cross(p3V1_BX_rotCS)).Unit();
505 float tmpSgnPhi = p3V1_BX_rotCS.Dot(normal1_BX_rotCS.Cross(normal2_BX_rotCS));
507 if (fabs(tmpSgnPhi)>0.) sgnPhi = tmpSgnPhi/fabs(tmpSgnPhi);
508 float dot_BX12 = normal1_BX_rotCS.Dot(normal2_BX_rotCS);
509 if (fabs(dot_BX12)>=1.) dot_BX12 *= 1./fabs(dot_BX12);
510 Phi = sgnPhi * acos(-1.*dot_BX12);
513 float tmpSgnPhi1 = p3V1_BX_rotCS.Dot(normal1_BX_rotCS.Cross(normalSC_BX_rotCS));
515 if (fabs(tmpSgnPhi1)>0.) sgnPhi1 = tmpSgnPhi1/fabs(tmpSgnPhi1);
516 float dot_BX1SC = normal1_BX_rotCS.Dot(normalSC_BX_rotCS);
517 if (fabs(dot_BX1SC)>=1.) dot_BX1SC *= 1./fabs(dot_BX1SC);
518 Phi1 = sgnPhi1 * acos(dot_BX1SC);
521 if (!(fabs(Z1_lept1Id)==21 || fabs(Z1_lept1Id)==22 || fabs(Z1_lept2Id)==21 || fabs(Z1_lept2Id)==22)){
523 TRotation rotationZ1;
524 TVector3 newZaxisZ1(thep4Z1inXFrame_rotCS.Vect().Unit());
525 TVector3 newXaxisZ1(newYaxisCS.Cross(newZaxisZ1).Unit());
526 TVector3 newYaxisZ1(newZaxisZ1.Cross(newXaxisZ1).Unit());
527 rotationZ1.RotateAxes(newXaxisZ1, newYaxisZ1, newZaxisZ1);
530 TLorentzVector thep4Z1inXFrame_rotCS_rotZ1(thep4Z1inXFrame_rotCS);
531 thep4Z1inXFrame_rotCS_rotZ1.Transform(rotationZ1);
532 TVector3 boostZ1inX_rotCS_rotZ1= -(thep4Z1inXFrame_rotCS_rotZ1.BoostVector());
534 TLorentzVector p4M11_BX_rotCS_rotZ1(p4M11_BX_rotCS);
535 TLorentzVector p4M12_BX_rotCS_rotZ1(p4M12_BX_rotCS);
536 TLorentzVector p4M21_BX_rotCS_rotZ1(p4M21_BX_rotCS);
537 TLorentzVector p4M22_BX_rotCS_rotZ1(p4M22_BX_rotCS);
538 p4M11_BX_rotCS_rotZ1.Transform(rotationZ1);
539 p4M12_BX_rotCS_rotZ1.Transform(rotationZ1);
540 p4M21_BX_rotCS_rotZ1.Transform(rotationZ1);
541 p4M22_BX_rotCS_rotZ1.Transform(rotationZ1);
542 p4M11_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
543 p4M12_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
544 p4M21_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
545 p4M22_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
547 TLorentzVector p4V2_BX_rotCS_rotZ1 = p4M21_BX_rotCS_rotZ1 + p4M22_BX_rotCS_rotZ1;
549 costheta1 = -p4V2_BX_rotCS_rotZ1.Vect().Dot(p4M11_BX_rotCS_rotZ1.Vect())/p4V2_BX_rotCS_rotZ1.Vect().Mag()/p4M11_BX_rotCS_rotZ1.Vect().Mag();
554 if (!(fabs(Z2_lept1Id)==21 || fabs(Z2_lept1Id)==22 || fabs(Z2_lept2Id)==21 || fabs(Z2_lept2Id)==22)){
556 TRotation rotationZ2;
557 TVector3 newZaxisZ2(thep4Z2inXFrame_rotCS.Vect().Unit());
558 TVector3 newXaxisZ2(newYaxisCS.Cross(newZaxisZ2).Unit());
559 TVector3 newYaxisZ2(newZaxisZ2.Cross(newXaxisZ2).Unit());
560 rotationZ2.RotateAxes(newXaxisZ2, newYaxisZ2, newZaxisZ2);
563 TLorentzVector thep4Z2inXFrame_rotCS_rotZ2(thep4Z2inXFrame_rotCS);
564 thep4Z2inXFrame_rotCS_rotZ2.Transform(rotationZ2);
565 TVector3 boostZ2inX_rotCS_rotZ2= -(thep4Z2inXFrame_rotCS_rotZ2.BoostVector());
567 TLorentzVector p4M11_BX_rotCS_rotZ2(p4M11_BX_rotCS);
568 TLorentzVector p4M12_BX_rotCS_rotZ2(p4M12_BX_rotCS);
569 TLorentzVector p4M21_BX_rotCS_rotZ2(p4M21_BX_rotCS);
570 TLorentzVector p4M22_BX_rotCS_rotZ2(p4M22_BX_rotCS);
571 p4M11_BX_rotCS_rotZ2.Transform(rotationZ2);
572 p4M12_BX_rotCS_rotZ2.Transform(rotationZ2);
573 p4M21_BX_rotCS_rotZ2.Transform(rotationZ2);
574 p4M22_BX_rotCS_rotZ2.Transform(rotationZ2);
575 p4M11_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
576 p4M12_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
577 p4M21_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
578 p4M22_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
581 TLorentzVector p4V1_BX_rotCS_rotZ2= p4M11_BX_rotCS_rotZ2 + p4M12_BX_rotCS_rotZ2;
583 costheta2 = -p4V1_BX_rotCS_rotZ2.Vect().Dot(p4M21_BX_rotCS_rotZ2.Vect())/p4V1_BX_rotCS_rotZ2.Vect().Mag()/p4M21_BX_rotCS_rotZ2.Vect().Mag();
588 MELAout <<
"WARNING: NaN in computeAngles: "
589 << costhetastar <<
" "
593 << Phi1 <<
" " << endl;
606 TLorentzVector p4M11,
int Z1_lept1Id,
607 TLorentzVector p4M12,
int Z1_lept2Id,
608 TLorentzVector p4M21,
int Z2_lept1Id,
609 TLorentzVector p4M22,
int Z2_lept2Id,
610 TLorentzVector jet1,
int jet1Id,
611 TLorentzVector jet2,
int jet2Id,
612 TLorentzVector* injet1,
int injet1Id,
613 TLorentzVector* injet2,
int injet2Id
615 TLorentzVector
const nullFourVector(0, 0, 0, 0);
617 TLorentzVector jet1massless, jet2massless;
619 jet1massless = jetPair.first;
620 jet2massless = jetPair.second;
623 TLorentzVector p4Z1 = p4M11 + p4M12;
624 TLorentzVector p4Z2 = p4M21 + p4M22;
625 TLorentzVector pH = p4Z1+p4Z2;
627 if (jet1massless.Z() < jet2massless.Z()) { swap(jet1massless, jet2massless); swap(jet1Id, jet2Id); }
633 TLorentzRotation movingframe;
634 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
635 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
636 movingframe.Boost(-pHJJ_perp.BoostVector());
637 pHJJ.Boost(-pHJJ_perp.BoostVector());
638 movingframe.Boost(-pHJJ.BoostVector());
639 pHJJ.Boost(-pHJJ.BoostVector());
641 TLorentzVector P1(0, 0, pHJJ.T()/2, pHJJ.T()/2);
642 TLorentzVector P2(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
644 P1.Transform(movingframe.Inverse());
645 P2.Transform(movingframe.Inverse());
646 pHJJ.Transform(movingframe.Inverse());
649 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
652 if (P1.Z() < P2.Z()){
654 swap(injet1Id, injet2Id);
657 int diff1Id = jet1Id-injet1Id;
658 int diff2Id = jet2Id-injet2Id;
661 (diff1Id==0 && diff2Id==0 && !(injet1Id==21 || injet2Id==21))
663 ((fabs(diff1Id)==1 || fabs(diff1Id)==3 || fabs(diff1Id)==5) && (fabs(diff2Id)==1 || fabs(diff2Id)==3 || fabs(diff2Id)==5))
666 int diff12Id = jet1Id-injet2Id;
667 int diff21Id = jet2Id-injet1Id;
669 ((diff12Id==0 || diff21Id==0) && !(injet1Id==21 || injet2Id==21))
671 ((fabs(diff12Id)==1 || fabs(diff12Id)==3 || fabs(diff12Id)==5) || (fabs(diff21Id)==1 || fabs(diff21Id)==3 || fabs(diff21Id)==5))
674 swap(injet1Id, injet2Id);
679 TLorentzRotation ZZframe;
680 ZZframe.Boost(-pH.BoostVector());
681 P1.Transform(ZZframe);
682 P2.Transform(ZZframe);
683 p4Z1.Transform(ZZframe);
684 p4Z2.Transform(ZZframe);
685 jet1massless.Transform(ZZframe);
686 jet2massless.Transform(ZZframe);
688 TLorentzVector fermion1, fermion2, antifermion1, antifermion2;
691 if (injet1 && injet1Id<0){
693 antifermion1 = jet1massless;
696 fermion1 = jet1massless;
699 if (injet2 && injet2Id<0){
701 antifermion2 = jet2massless;
704 fermion2 = jet2massless;
709 TLorentzVector V1 = fermion1 + antifermion1;
710 TLorentzVector V2 = fermion2 + antifermion2;
711 TVector3 normvec1 = fermion1.Vect().Cross(antifermion1.Vect()).Unit();
712 TVector3 normvec2 = fermion2.Vect().Cross(antifermion2.Vect()).Unit();
713 TVector3 normvec3 = p4Z2.Vect().Cross(V1.Vect()).Unit();
714 double cosPhi = normvec1.Dot(normvec2);
715 double sgnPhi = normvec1.Cross(normvec2).Dot(V1.Vect());
716 if (fabs(sgnPhi)>0.) sgnPhi = sgnPhi/fabs(sgnPhi);
717 double cosPhi1 = normvec1.Dot(normvec3);
718 double sgnPhi1 = normvec1.Cross(normvec3).Dot(V1.Vect());
719 if (fabs(sgnPhi1)>0.) sgnPhi1 = sgnPhi1/fabs(sgnPhi1);
720 if (fabs(cosPhi)>1) cosPhi *= 1./fabs(cosPhi);
721 if (fabs(cosPhi1)>1) cosPhi1 *= 1./fabs(cosPhi1);
722 Phi = acos(-cosPhi)*sgnPhi;
723 Phi1 = acos(cosPhi1)*sgnPhi1;
724 costhetastar = V1.Vect().Unit().Dot(p4Z2.Vect().Unit());
729 costheta1 = V1.Vect().Unit().Dot(fermion1.Vect().Unit());
730 costheta2 = V2.Vect().Unit().Dot(fermion2.Vect().Unit());
734 float& costheta1_real,
float& costheta1_imag,
735 float& costheta2_real,
float& costheta2_imag,
740 TLorentzVector p4M11,
int Z1_lept1Id,
741 TLorentzVector p4M12,
int Z1_lept2Id,
742 TLorentzVector p4M21,
int Z2_lept1Id,
743 TLorentzVector p4M22,
int Z2_lept2Id,
744 TLorentzVector jet1,
int jet1Id,
745 TLorentzVector jet2,
int jet2Id,
746 TLorentzVector* injet1,
int injet1Id,
747 TLorentzVector* injet2,
int injet2Id
749 TLorentzVector
const nullFourVector(0, 0, 0, 0);
751 TLorentzVector jet1massless, jet2massless;
753 jet1massless = jetPair.first;
754 jet2massless = jetPair.second;
757 TLorentzVector p4Z1 = p4M11 + p4M12;
758 TLorentzVector p4Z2 = p4M21 + p4M22;
759 TLorentzVector pH = p4Z1+p4Z2;
761 if (jet1massless.Z() < jet2massless.Z()) { swap(jet1massless, jet2massless); swap(jet1Id, jet2Id); }
767 TLorentzRotation movingframe;
768 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
769 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
770 movingframe.Boost(-pHJJ_perp.BoostVector());
771 pHJJ.Boost(-pHJJ_perp.BoostVector());
772 movingframe.Boost(-pHJJ.BoostVector());
773 pHJJ.Boost(-pHJJ.BoostVector());
775 TLorentzVector P1(0, 0, pHJJ.T()/2, pHJJ.T()/2);
776 TLorentzVector P2(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
778 P1.Transform(movingframe.Inverse());
779 P2.Transform(movingframe.Inverse());
780 pHJJ.Transform(movingframe.Inverse());
783 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
786 if (P1.Z() < P2.Z()){
788 swap(injet1Id, injet2Id);
791 int diff1Id = jet1Id-injet1Id;
792 int diff2Id = jet2Id-injet2Id;
795 (diff1Id==0 && diff2Id==0 && !(injet1Id==21 || injet2Id==21))
797 ((fabs(diff1Id)==1 || fabs(diff1Id)==3 || fabs(diff1Id)==5) && (fabs(diff2Id)==1 || fabs(diff2Id)==3 || fabs(diff2Id)==5))
800 int diff12Id = jet1Id-injet2Id;
801 int diff21Id = jet2Id-injet1Id;
803 ((diff12Id==0 || diff21Id==0) && !(injet1Id==21 || injet2Id==21))
805 ((fabs(diff12Id)==1 || fabs(diff12Id)==3 || fabs(diff12Id)==5) || (fabs(diff21Id)==1 || fabs(diff21Id)==3 || fabs(diff21Id)==5))
808 swap(injet1Id, injet2Id);
813 TLorentzRotation ZZframe;
814 ZZframe.Boost(-pH.BoostVector());
815 P1.Transform(ZZframe);
816 P2.Transform(ZZframe);
817 p4Z1.Transform(ZZframe);
818 p4Z2.Transform(ZZframe);
819 jet1massless.Transform(ZZframe);
820 jet2massless.Transform(ZZframe);
822 TLorentzVector fermion1, fermion2, antifermion1, antifermion2;
825 if (injet1 && injet1Id<0){
827 antifermion1 = jet1massless;
830 fermion1 = jet1massless;
833 if (injet2 && injet2Id<0){
835 antifermion2 = jet2massless;
838 fermion2 = jet2massless;
843 TLorentzVector V1 = fermion1 + antifermion1;
844 TLorentzVector V2 = fermion2 + antifermion2;
845 TVector3 normvec1 = fermion1.Vect().Cross(antifermion1.Vect()).Unit();
846 TVector3 normvec2 = fermion2.Vect().Cross(antifermion2.Vect()).Unit();
847 TVector3 normvec3 = p4Z2.Vect().Cross(V1.Vect()).Unit();
848 double cosPhi = normvec1.Dot(normvec2);
849 double sgnPhi = normvec1.Cross(normvec2).Dot(V1.Vect());
850 if (fabs(sgnPhi)>0.) sgnPhi = sgnPhi/fabs(sgnPhi);
851 double cosPhi1 = normvec1.Dot(normvec3);
852 double sgnPhi1 = normvec1.Cross(normvec3).Dot(V1.Vect());
853 if (fabs(sgnPhi1)>0.) sgnPhi1 = sgnPhi1/fabs(sgnPhi1);
854 if (fabs(cosPhi)>1) cosPhi *= 1./fabs(cosPhi);
855 if (fabs(cosPhi1)>1) cosPhi1 *= 1./fabs(cosPhi1);
856 Phi = acos(-cosPhi)*sgnPhi;
857 Phi1 = acos(cosPhi1)*sgnPhi1;
858 costhetastar = V1.Vect().Unit().Dot(p4Z2.Vect().Unit());
866 pair<TLorentzVector, TLorentzVector> fermion1_BV1 =
TUtil::ComplexBoost(V1.BoostVector(), fermion1);
867 costheta1_real = -(V2_BV1.first.Vect().Unit().Dot(fermion1_BV1.first.Vect().Unit()) - V2_BV1.second.Vect().Unit().Dot(fermion1_BV1.second.Vect().Unit()));
868 costheta1_imag = -(V2_BV1.first.Vect().Unit().Dot(fermion1_BV1.second.Vect().Unit()) + V2_BV1.second.Vect().Unit().Dot(fermion1_BV1.first.Vect().Unit()));
871 pair<TLorentzVector, TLorentzVector> fermion2_BV2 =
TUtil::ComplexBoost(V2.BoostVector(), fermion2);
872 costheta2_real = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()) - V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()));
873 costheta2_imag = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()) + V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()));
883 TLorentzVector p4M11,
int Z1_lept1Id,
884 TLorentzVector p4M12,
int Z1_lept2Id,
885 TLorentzVector p4M21,
int Z2_lept1Id,
886 TLorentzVector p4M22,
int Z2_lept2Id,
887 TLorentzVector jet1,
int jet1Id,
888 TLorentzVector jet2,
int jet2Id,
889 TLorentzVector* injet1,
int injet1Id,
890 TLorentzVector* injet2,
int injet2Id
892 TLorentzVector
const nullFourVector(0, 0, 0, 0);
894 TLorentzVector jet1massless, jet2massless;
896 jet1massless = jetPair.first;
897 jet2massless = jetPair.second;
900 TLorentzVector p4Z1 = p4M11 + p4M12;
901 TLorentzVector p4Z2 = p4M21 + p4M22;
902 TLorentzVector pH = p4Z1 + p4Z2;
906 (jet1Id*jet2Id<0 && jet1Id<0)
908 ((jet1Id*jet2Id>0 || (jet1Id==0 && jet2Id==0)) && jet1massless.Phi()<=jet2massless.Phi())
910 swap(jet1massless, jet2massless);
911 swap(jet1Id, jet2Id);
919 TLorentzRotation movingframe;
920 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
921 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
922 movingframe.Boost(-pHJJ_perp.BoostVector());
923 pHJJ.Boost(-pHJJ_perp.BoostVector());
924 movingframe.Boost(-pHJJ.BoostVector());
925 pHJJ.Boost(-pHJJ.BoostVector());
927 TLorentzVector P1(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
928 TLorentzVector P2(0, 0, pHJJ.T()/2, pHJJ.T()/2);
930 P1.Transform(movingframe.Inverse());
931 P2.Transform(movingframe.Inverse());
932 pHJJ.Transform(movingframe.Inverse());
935 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
940 (injet1Id*injet2Id<0 && injet1Id>0)
942 (injet1Id*injet2Id>0 && P1.Z()>=P2.Z())
945 swap(injet1Id, injet2Id);
950 TLorentzRotation ZZframe;
951 TVector3
const beamAxis(0, 0, 1);
952 if (p4Z1==nullFourVector || p4Z2==nullFourVector){
953 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
954 if (pNewAxis != nullFourVector.Vect()){
955 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
956 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
958 P1.Transform(ZZframe);
959 P2.Transform(ZZframe);
960 jet1massless = -jet1massless; jet1massless.Transform(ZZframe); jet1massless = -jet1massless;
961 jet2massless = -jet2massless; jet2massless.Transform(ZZframe); jet2massless = -jet2massless;
965 ZZframe.Boost(-pH.BoostVector());
966 p4Z1.Boost(-pH.BoostVector());
967 p4Z2.Boost(-pH.BoostVector());
968 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
969 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
970 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
971 P1.Transform(ZZframe);
972 P2.Transform(ZZframe);
973 jet1massless = -jet1massless; jet1massless.Transform(ZZframe); jet1massless = -jet1massless;
974 jet2massless = -jet2massless; jet2massless.Transform(ZZframe); jet2massless = -jet2massless;
989 m2 = (jet1massless + jet2massless).M();
1013 TLorentzVector p4M11,
int Z1_lept1Id,
1014 TLorentzVector p4M12,
int Z1_lept2Id,
1015 TLorentzVector p4M21,
int Z2_lept1Id,
1016 TLorentzVector p4M22,
int Z2_lept2Id,
1017 TLorentzVector b,
int bId,
1018 TLorentzVector Wplusf,
int WplusfId,
1019 TLorentzVector Wplusfb,
int WplusfbId,
1020 TLorentzVector bbar,
int bbarId,
1021 TLorentzVector Wminusf,
int WminusfId,
1022 TLorentzVector Wminusfb,
int WminusfbId,
1023 TLorentzVector* injet1,
int injet1Id,
1024 TLorentzVector* injet2,
int injet2Id
1026 TLorentzVector
const nullFourVector(0, 0, 0, 0);
1027 TVector3
const beamAxis(0, 0, 1);
1028 TVector3
const xAxis(1, 0, 0);
1031 hs=hincoming=hTT=PhiTT=Phi1
1032 =hbb=hWW=Phibb=Phi1bb
1034 =hWminusf=PhiWminusf=0;
1040 float m1_dummy=0, m2_dummy=0;
1042 hs, hincoming, hTT, PhiTT, Phi1,
1057 WplusfId!=-9000 && WplusfbId!=-9000
1060 (WplusfId*WplusfbId<0 && WplusfId<0)
1062 ((WplusfId*WplusfbId>0 || (WplusfId==0 && WplusfbId==0)) && Wplusf.Phi()<=Wplusfb.Phi())
1065 swap(Wplusf, Wplusfb);
1066 swap(WplusfId, WplusfbId);
1070 WminusfId!=-9000 && WminusfbId!=-9000
1073 (WminusfId*WminusfbId<0 && WminusfId<0)
1075 ((WminusfId*WminusfbId>0 || (WminusfId==0 && WminusfbId==0)) && Wminusf.Phi()<=Wminusfb.Phi())
1078 swap(Wminusf, Wminusfb);
1079 swap(WminusfId, WminusfbId);
1083 if (bId!=-9000 && WplusfId!=-9000 && WplusfbId!=-9000){
1085 tmp_daus.reserve(3);
1086 tmp_daus.emplace_back(bId, b);
1087 tmp_daus.emplace_back(WplusfId, Wplusf);
1088 tmp_daus.emplace_back(WplusfbId, Wplusfb);
1090 b = tmp_daus.at(0).second;
1091 Wplusf = tmp_daus.at(1).second;
1092 Wplusfb = tmp_daus.at(2).second;
1094 else if (WplusfId!=-9000 && WplusfbId!=-9000){
1096 Wplusf = Wffb.first;
1097 Wplusfb = Wffb.second;
1100 if (bbarId!=-9000 && WminusfId!=-9000 && WminusfbId!=-9000){
1102 tmp_daus.reserve(3);
1103 tmp_daus.emplace_back(bbarId, bbar);
1104 tmp_daus.emplace_back(WminusfId, Wminusf);
1105 tmp_daus.emplace_back(WminusfbId, Wminusfb);
1107 bbar = tmp_daus.at(0).second;
1108 Wminusf = tmp_daus.at(1).second;
1109 Wminusfb = tmp_daus.at(2).second;
1111 else if (WminusfId!=-9000 && WminusfbId!=-9000){
1113 Wminusf = Wffb.first;
1114 Wminusfb = Wffb.second;
1119 TLorentzVector p4Z1 = p4M11 + p4M12;
1120 TLorentzVector p4Z2 = p4M21 + p4M22;
1122 TLorentzVector pH = p4Z1 + p4Z2;
1124 TLorentzVector pWplus = Wplusf + Wplusfb;
1125 TLorentzVector pTop = b + pWplus;
1127 TLorentzVector pWminus = Wminusf + Wminusfb;
1128 TLorentzVector pATop = bbar + pWminus;
1130 TLorentzVector pTT = pTop + pATop;
1131 TLorentzVector pTTH = pTT + pH;
1138 TLorentzRotation movingframe;
1139 TLorentzVector pTTH_perp(pTTH.X(), pTTH.Y(), 0, pTTH.T());
1140 movingframe.Boost(-pTTH_perp.BoostVector());
1141 pTTH.Boost(-pTTH_perp.BoostVector());
1142 movingframe.Boost(-pTTH.BoostVector());
1143 pTTH.Boost(-pTTH.BoostVector());
1145 TLorentzVector P1(0, 0, -pTTH.T()/2, pTTH.T()/2);
1146 TLorentzVector P2(0, 0, pTTH.T()/2, pTTH.T()/2);
1148 P1.Transform(movingframe.Inverse());
1149 P2.Transform(movingframe.Inverse());
1150 pTTH.Transform(movingframe.Inverse());
1153 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pTTH.P())<pTTH.P()*1e-4){
1158 (injet1Id*injet2Id<0 && injet1Id>0)
1160 (injet1Id*injet2Id>0 && P1.Z()>=P2.Z())
1163 swap(injet1Id, injet2Id);
1168 TLorentzRotation ZZframe;
1169 bool applyZZframe=
false;
1170 if (p4Z1==nullFourVector || p4Z2==nullFourVector){
1171 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
1172 if (pNewAxis != nullFourVector.Vect()){
1173 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
1174 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
1179 TVector3 pHboost = pH.BoostVector();
1180 ZZframe.Boost(-pHboost);
1181 p4Z1.Boost(-pHboost);
1182 p4Z2.Boost(-pHboost);
1183 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
1184 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
1185 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
1187 p4Z1.Boost(pHboost);
1188 p4Z2.Boost(pHboost);
1192 P1.Transform(ZZframe); P2.Transform(ZZframe);
1193 p4Z1 = -p4Z1; p4Z1.Transform(ZZframe); p4Z1 = -p4Z1;
1194 p4Z2 = -p4Z2; p4Z2.Transform(ZZframe); p4Z2 = -p4Z2;
1195 b = -b; b.Transform(ZZframe); b = -b;
1196 Wplusf = -Wplusf; Wplusf.Transform(ZZframe); Wplusf = -Wplusf;
1197 Wplusfb = -Wplusfb; Wplusfb.Transform(ZZframe); Wplusfb = -Wplusfb;
1198 bbar = -bbar; bbar.Transform(ZZframe); bbar = -bbar;
1199 Wminusf = -Wminusf; Wminusf.Transform(ZZframe); Wminusf = -Wminusf;
1200 Wminusfb = -Wminusfb; Wminusfb.Transform(ZZframe); Wminusfb = -Wminusfb;
1204 pWplus = Wplusf + Wplusfb;
1206 pWminus = Wminusf + Wminusfb;
1207 pATop = bbar + pWminus;
1226 TVector3 pHboost = pH.BoostVector();
1234 Wplusf.Boost(-pHboost);
1235 Wplusfb.Boost(-pHboost);
1236 bbar.Boost(-pHboost);
1237 Wminusf.Boost(-pHboost);
1238 Wminusfb.Boost(-pHboost);
1243 pWplus = Wplusf + Wplusfb;
1245 pWminus = Wminusf + Wminusfb;
1246 pATop = bbar + pWminus;
1253 TLorentzRotation TTframe;
1256 TVector3 nNewZAxis = (-P1-P2-pTop-pATop).Vect().Unit();
1257 if (nNewZAxis != nullFourVector.Vect()){
1258 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1259 TTframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1262 P1.Transform(TTframe);
1263 P2.Transform(TTframe);
1264 b.Transform(TTframe);
1265 Wplusf.Transform(TTframe);
1266 Wplusfb.Transform(TTframe);
1267 bbar.Transform(TTframe);
1268 Wminusf.Transform(TTframe);
1269 Wminusfb.Transform(TTframe);
1272 pWplus = Wplusf + Wplusfb;
1274 pWminus = Wminusf + Wminusfb;
1275 pATop = bbar + pWminus;
1281 int id_dummy_Wplus=24;
1282 int id_dummy_Wminus=-24;
1283 if (WplusfId!=-9000 && WplusfbId!=-9000) id_dummy_Wplus=-9000;
1284 if (WminusfId!=-9000 && WminusfbId!=-9000) id_dummy_Wminus=-9000;
1293 pWplus, id_dummy_Wplus,
1294 pWminus, id_dummy_Wminus
1299 if (WplusfId!=-9000 && WplusfbId!=-9000){
1300 TLorentzVector pATop_tmp = pATop;
1301 TVector3 pWplus_boost = (Wplusf+Wplusfb).BoostVector();
1302 b.Boost(-pWplus_boost);
1303 Wplusf.Boost(-pWplus_boost);
1304 Wplusfb.Boost(-pWplus_boost);
1305 pATop_tmp.Boost(-pWplus_boost);
1307 TLorentzRotation Wplusframe;
1310 TVector3 nNewZAxis = (b+Wplusf+Wplusfb-pATop_tmp).Vect().Unit();
1311 if (nNewZAxis != nullFourVector.Vect()){
1312 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1313 Wplusframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1316 b.Transform(Wplusframe);
1317 Wplusf.Transform(Wplusframe);
1318 Wplusfb.Transform(Wplusframe);
1320 TVector3 n3_Wb = -b.Vect().Unit();
1321 TVector3 nW = ((Wplusf-Wplusfb).Vect().Cross(n3_Wb)).Unit();
1322 TVector3 nS = (beamAxis.Cross(n3_Wb)).Unit();
1325 hWplusf = n3_Wb.Dot(Wplusf.Vect().Unit());
1327 float tmpSgnPhiWplusf = n3_Wb.Dot(nW.Cross(nS));
1328 float sgnPhiWplusf = 0;
1329 if (fabs(tmpSgnPhiWplusf)>0.) sgnPhiWplusf = tmpSgnPhiWplusf/fabs(tmpSgnPhiWplusf);
1330 float dot_nWnS = nW.Dot(nS);
1331 if (fabs(dot_nWnS)>=1.) dot_nWnS *= 1./fabs(dot_nWnS);
1332 PhiWplusf = sgnPhiWplusf * acos(dot_nWnS);
1341 if (WminusfId!=-9000 && WminusfbId!=-9000){
1342 TLorentzVector pTop_tmp = pTop;
1343 TVector3 pWminus_boost = (Wminusf+Wminusfb).BoostVector();
1344 bbar.Boost(-pWminus_boost);
1345 Wminusf.Boost(-pWminus_boost);
1346 Wminusfb.Boost(-pWminus_boost);
1347 pTop_tmp.Boost(-pWminus_boost);
1349 TLorentzRotation Wminusframe;
1352 TVector3 nNewZAxis = (pTop_tmp-bbar-Wminusf-Wminusfb).Vect().Unit();
1353 if (nNewZAxis != nullFourVector.Vect()){
1354 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1355 Wminusframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1358 bbar.Transform(Wminusframe);
1359 Wminusf.Transform(Wminusframe);
1360 Wminusfb.Transform(Wminusframe);
1362 TVector3 n3_Wb = -bbar.Vect().Unit();
1363 TVector3 nW = ((Wminusf-Wminusfb).Vect().Cross(n3_Wb)).Unit();
1364 TVector3 nS = (beamAxis.Cross(n3_Wb)).Unit();
1367 hWminusf = n3_Wb.Dot(Wminusf.Vect().Unit());
1369 float tmpSgnPhiWminusf = n3_Wb.Dot(nW.Cross(nS));
1370 float sgnPhiWminusf = 0;
1371 if (fabs(tmpSgnPhiWminusf)>0.) sgnPhiWminusf = tmpSgnPhiWminusf/fabs(tmpSgnPhiWminusf);
1372 float dot_nWnS = nW.Dot(nS);
1373 if (fabs(dot_nWnS)>=1.) dot_nWnS *= 1./fabs(dot_nWnS);
1374 PhiWminusf = sgnPhiWminusf * acos(dot_nWnS);
1390 const double GeV=1./100.;
1391 double ext_mZ_jhu = ext_mZ*GeV;
1392 double ext_mW_jhu = ext_mW*GeV;
1393 double ext_Gf_jhu = ext_Gf/pow(GeV, 2);
1397 if (ext_ewscheme<-1 || ext_ewscheme>3) ext_ewscheme=3;
1448 const int ipartabs = abs(ipart);
1449 bool runcoupling_mcfm=
false;
1450 bool runcoupling_jhugen=
false;
1473 (ipartabs>=11 && ipartabs<=16)
1475 ipartabs==23 || ipartabs==24 || ipartabs==25
1477 ipartabs==32 || ipartabs==34
1479 runcoupling_jhugen=(
1484 const double GeV=1./100.;
1485 double jinmass = inmass*GeV;
1491 if (runcoupling_mcfm || runcoupling_jhugen){
1498 const int ipartabs = abs(ipart);
1508 const double GeV=1./100.;
1509 double jinwidth = inwidth*GeV;
1513 void TUtil::SetCKMElements(
double* invckm_ud,
double* invckm_us,
double* invckm_cd,
double* invckm_cs,
double* invckm_ts,
double* invckm_tb,
double* invckm_ub,
double* invckm_cb,
double* invckm_td){
1545 }
else if(jquark == 2){
1547 }
else if(jquark == 3){
1550 MELAerr <<
"TUtil::GetMadgraphCKMElement: Invalid second index!" << endl;
1556 }
else if(jquark == 2){
1558 }
else if(jquark == 3){
1561 MELAerr <<
"TUtil::GetMadgraphCKMElement: Invalid second index!" << endl;
1567 }
else if(jquark == 2){
1569 }
else if(jquark == 3){
1572 MELAerr <<
"TUtil::GetMadgraphCKMElement: Invalid second index!" << endl;
1576 MELAerr <<
"TUtil::GetMadgraphCKMElement: Invalid first index!" << endl;
1579 return std::complex<double>(0);
1583 const int ipartabs = abs(ipart);
1604 (ipartabs>=11 && ipartabs<=16)
1606 ipartabs==23 || ipartabs==24 || ipartabs==25
1608 ipartabs==32 || ipartabs==34
1610 const double GeV=1./100.;
1613 double outmass = joutmass/GeV;
1620 const int ipartabs = abs(ipart);
1630 const double GeV=1./100.;
1633 double outwidth = joutwidth/GeV;
1638 if (part==
nullptr)
return GetMass(-9000);
1656 TLorentzVector
const nullFourVector(0, 0, 0, 0);
1665 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5];
1666 Q = fabs(pTotal.M());
1669 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5]+p[6]+p[7];
1670 Q = fabs(pTotal.M());
1673 TLorentzVector qH = p[2]+p[3]+p[4]+p[5];
1674 TLorentzVector qJJ = p[6]+p[7];
1675 Q = fabs(qH.M()+qJJ.M());
1678 TLorentzVector qH = p[2]+p[3]+p[4]+p[5];
1679 Q = fabs(qH.M()+p[6].M()+p[7].M());
1682 for (
int c=2; c<
mxpart; c++)
Q += p[c].Pt();
1686 for (
int c=6; c<
mxpart; c++)
Q = std::max(
Q, p[c].Pt());
1691 for (
int c=7; c<
mxpart; c++){
if (p[c].Pt()>0.)
Q = std::min(
Q, p[c].Pt()); }
1711 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5];
1712 Q = fabs(pTotal.M());
1756 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5];
1757 Q = fabs(pTotal.M());
1763 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5]+p[6]+p[7];
1764 Q = fabs(pTotal.M());
1770 MELAerr <<
"Scaling " << scheme <<
" fails for production " << production <<
", defaulting to dynamic scheme m3456 " << endl;
1771 TLorentzVector pTotal = p[2]+p[3]+p[4]+p[5];
1772 Q = fabs(pTotal.M());
1776 void TUtil::SetAlphaS(
double& Q_ren,
double& Q_fac,
double multiplier_ren,
double multiplier_fac,
int mynloop,
int mynflav,
string mypartons){
1777 bool hasReset=
false;
1778 if (multiplier_ren<=0. || multiplier_fac<=0.){
1779 MELAerr <<
"TUtil::SetAlphaS: Invalid scale multipliers" << endl;
1782 if (Q_ren<=1. || Q_fac<=1. || mynloop<=0 || mypartons.compare(
"Default")==0){
1783 if (Q_ren<0.)
MELAout <<
"TUtil::SetAlphaS: Invalid QCD scale for alpha_s, setting to mH/2..." << endl;
1784 if (Q_fac<0.)
MELAout <<
"TUtil::SetAlphaS: Invalid factorization scale, setting to mH/2..." << endl;
1791 Q_ren *= multiplier_ren;
1792 Q_fac *= multiplier_fac;
1796 bool nflav_is_same = (
nflav_.nflav == mynflav);
1797 if (!nflav_is_same)
MELAout <<
"TUtil::SetAlphaS: nflav=" <<
nflav_.nflav <<
" is the only one supported." << endl;
1799 scale_.musq = Q_ren*Q_ren;
1802 MELAout <<
"TUtil::SetAlphaS: Only nloop=1 is supported!" << endl;
1830 const double GeV=1./100.;
1831 double muren_jhu =
scale_.scale*GeV;
1832 double mufac_jhu =
facscale_.facscale*GeV;
1850 double alphasVal=0, alphasmzVal=0;
1853 if (alphasmz_!=0) *alphasmz_ = alphasmzVal;
1864 unsigned int ndau = mela_event.
pDaughters.size();
1865 int* pId =
new int[ndau];
1866 for (
unsigned int ip=0; ip<ndau; ip++){
1876 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::MCFM_chooser: isWW=" << (int)isWW <<
", isZZ=" << (
int)isZZ <<
", hasZZ4fInterf=" << (int)hasZZ4fInterf << endl;
1938 (ndau>=4 || ndau==2)
2265 MELAerr <<
"TUtil::MCFM_chooser: MCFM does not support QCD JJ+VV->4f interaction yet. Please contact the MELA authors if you need more information." << endl;
2270 MELAerr <<
"TUtil::MCFM_chooser: Can't identify (process, production) = (" <<
process <<
", " << production <<
")" << endl;
2271 MELAerr <<
"TUtil::MCFM_chooser: ndau: " << ndau <<
'\t';
2272 MELAerr <<
"TUtil::MCFM_chooser: isZZ: " << isZZ <<
'\t';
2273 MELAerr <<
"TUtil::MCFM_chooser: isWW: " << isWW <<
'\t';
2274 MELAerr <<
"TUtil::MCFM_chooser: isGG: " << isGG <<
'\t';
2286 std::vector<int>* partOrder, std::vector<int>* apartOrder
2299 TString strplabel[
mxpart];
2300 for (
int ip=0; ip<
mxpart; ip++) strplabel[ip]=
" ";
2303 unsigned int ndau = mela_event.
pDaughters.size();
2304 if (ndau<1)
return false;
2305 unsigned int napart = mela_event.
pAssociated.size();
2320 bool hasZ1 = (isZZ || isZG || isZJJ);
2326 vector<TString> strcfgs;
2327 if (isGG) strcfgs.push_back(TString(
"GG"));
2328 if (isZG) strcfgs.push_back(TString(
"ZG"));
2329 if (isZZ) strcfgs.push_back(TString(
"ZZ"));
2330 if (isZJJ) strcfgs.push_back(TString(
"ZJJ"));
2331 if (isWW) strcfgs.push_back(TString(
"WW"));
2332 if (!strcfgs.empty()){
2334 for (
unsigned int istr=0; istr<strcfgs.size(); istr++){
2335 if (istr==0)
MELAout << strcfgs.at(istr);
2336 else MELAout <<
", " << strcfgs.at(istr);
2340 else MELAout <<
"has no valid decay!";
2345 int* pApartOrder = 0;
2348 pApartOrder =
new int[napart];
2349 pApartId =
new int[napart];
2350 for (
unsigned int ip=0; ip<napart; ip++){
2355 int* pOrder =
new int[ndau];
2356 int* pZOrder =
new int[ndau];
2357 int* pWOrder =
new int[ndau];
2358 int* pId =
new int[ndau];
2359 for (
unsigned int ip=0; ip<ndau; ip++){
2384 bool useQQVVQQstrong =
2390 bool useQQVVQQboth =
2396 bool useQQVVQQany = useQQVVQQ || useQQVVQQstrong || useQQVVQQboth;
2403 !(isWW || isZZ || isZJJ || isZG || isGG)
2409 if (isWW && !hasZ1 && !hasZ2){
2413 swap(pZOrder[0], pZOrder[2]);
2422 if (!hasZ1 && hasZ2){
2423 swap(pZOrder[0], pZOrder[2]);
2424 swap(pZOrder[1], pZOrder[3]);
2428 if (hasZ1)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Found a Z1";
2429 if (hasZ2)
MELAout <<
" and a Z2";
2430 if (hasZ1 || hasZ2)
MELAout <<
" in WW." << endl;
2433 if (isZZ && !hasW1 && !hasW2){
2434 swap(pWOrder[0], pWOrder[2]);
2442 if (V1id<0 || V2id>0){
2443 swap(pWOrder[0], pWOrder[2]);
2444 swap(pWOrder[1], pWOrder[3]);
2450 if (hasW1)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Found a W1(" << V1id <<
")";
2451 if (hasW2)
MELAout <<
" and a W2(" << V2id <<
")";
2452 if (hasW1 || hasW2)
MELAout <<
" in ZZ." << endl;
2464 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Setting up mother labels for MCFM:";
2465 for (
int ip=0; ip<min(2, (
int)mela_event.
pMothers.size()); ip++){
2466 const int& idmot = mela_event.
pMothers.at(ip).first;
2497 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2498 sprintf(
runstring_.runstring, strrun.c_str());
2525 else result =
false;
2537 else result =
false;
2549 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2550 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2551 sprintf(
runstring_.runstring, strrun.c_str());
2556 for (
unsigned int ix=0; ix<napart; ix++){
2558 for (
unsigned int iy=ix+1; iy<napart; iy++){
2562 int i1=ix;
int i2=iy;
2563 if (pApartId[ix]<0){ swap(i1, i2); }
2564 int* tmpOrder =
new int[napart];
2565 tmpOrder[0] = i1; tmpOrder[1] = i2;
2567 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2568 delete[] pApartOrder; pApartOrder = tmpOrder;
2573 int i1=ix;
int i2=iy;
2574 if (pApartId[ix]<0){ swap(i1, i2); }
2575 int* tmpOrder =
new int[napart];
2576 tmpOrder[0] = i1; tmpOrder[1] = i2;
2578 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2579 delete[] pApartOrder; pApartOrder = tmpOrder;
2583 if (hasZll || hasZnn)
break;
2586 unsigned int iout=6;
2587 unsigned int jout=7;
2589 int qqvvqq_apartordering[2]={ -1, -1 };
2592 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2593 qqvvqq_apartordering
2595 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2596 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]){
2598 swap(pApartOrder[0], pApartOrder[1]);
2602 strplabel[iout]=
"el";
2603 strplabel[jout]=
"ea";
2606 strplabel[iout]=
"nl";
2607 strplabel[jout]=
"na";
2609 else result =
false;
2617 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2618 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2619 sprintf(
runstring_.runstring, strrun.c_str());
2622 bool hasWplus=
false;
2623 bool hasWminus=
false;
2624 for (
unsigned int ix=0; ix<napart; ix++){
2626 for (
unsigned int iy=ix+1; iy<napart; iy++){
2638 int i1=ix;
int i2=iy;
2639 if (pApartId[ix]<0){ swap(i1, i2); }
2640 int* tmpOrder =
new int[napart];
2641 tmpOrder[0] = i1; tmpOrder[1] = i2;
2643 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2644 delete[] pApartOrder; pApartOrder = tmpOrder;
2657 int i1=ix;
int i2=iy;
2658 if (pApartId[ix]<0){ swap(i1, i2); }
2659 int* tmpOrder =
new int[napart];
2660 tmpOrder[0] = i1; tmpOrder[1] = i2;
2662 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2663 delete[] pApartOrder; pApartOrder = tmpOrder;
2668 if (hasWplus || hasWminus)
break;
2670 unsigned int iout=6;
2671 unsigned int jout=7;
2673 int qqvvqq_apartordering[2]={ -1, -1 };
2676 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2677 qqvvqq_apartordering
2679 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2680 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]){
2682 swap(pApartOrder[0], pApartOrder[1]);
2686 strplabel[iout]=
"nl";
2687 strplabel[jout]=
"ea";
2689 else if (hasWminus){
2690 strplabel[iout]=
"el";
2691 strplabel[jout]=
"na";
2693 else result =
false;
2701 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2702 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2703 sprintf(
runstring_.runstring, strrun.c_str());
2709 for (
unsigned int ix=0; ix<napart; ix++){
2711 for (
unsigned int iy=ix+1; iy<napart; iy++){
2717 int i1=ix;
int i2=iy;
2718 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2719 int* tmpOrder =
new int[napart];
2720 tmpOrder[0] = i1; tmpOrder[1] = i2;
2722 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2723 delete[] pApartOrder; pApartOrder = tmpOrder;
2728 int i1=ix;
int i2=iy;
2729 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2730 int* tmpOrder =
new int[napart];
2731 tmpOrder[0] = i1; tmpOrder[1] = i2;
2733 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2734 delete[] pApartOrder; pApartOrder = tmpOrder;
2739 int i1=ix;
int i2=iy;
2740 int* tmpOrder =
new int[napart];
2741 tmpOrder[0] = i1; tmpOrder[1] = i2;
2743 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2744 delete[] pApartOrder; pApartOrder = tmpOrder;
2748 if (hasZuu || hasZdd || hasZjj)
break;
2752 int qqvvqq_apartordering[2]={ -1, -1 };
2755 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2756 qqvvqq_apartordering
2758 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2759 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2761 if (hasZuu || hasZdd){
2769 else result =
false;
2777 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2778 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2779 sprintf(
runstring_.runstring, strrun.c_str());
2782 bool hasWplus=
false;
2783 bool hasWminus=
false;
2785 for (
unsigned int ix=0; ix<napart; ix++){
2787 for (
unsigned int iy=ix+1; iy<napart; iy++){
2805 int i1=ix;
int i2=iy;
2806 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2807 int* tmpOrder =
new int[napart];
2808 tmpOrder[0] = i1; tmpOrder[1] = i2;
2810 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2811 delete[] pApartOrder; pApartOrder = tmpOrder;
2824 int i1=ix;
int i2=iy;
2825 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2826 int* tmpOrder =
new int[napart];
2827 tmpOrder[0] = i1; tmpOrder[1] = i2;
2829 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2830 delete[] pApartOrder; pApartOrder = tmpOrder;
2835 int i1=ix;
int i2=iy;
2836 int* tmpOrder =
new int[napart];
2837 tmpOrder[0] = i1; tmpOrder[1] = i2;
2839 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2840 delete[] pApartOrder; pApartOrder = tmpOrder;
2845 if (hasWplus || hasWminus || hasWjj)
break;
2848 int qqvvqq_apartordering[2]={ -1, -1 };
2851 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2852 qqvvqq_apartordering
2854 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2855 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2857 if (hasWplus || hasWminus){
2865 else result =
false;
2873 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2874 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2875 sprintf(
runstring_.runstring, strrun.c_str());
2879 int qqvvqq_apartordering[2]={ -1, -1 };
2882 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2883 qqvvqq_apartordering
2885 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2886 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2897 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2898 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2899 sprintf(
runstring_.runstring, strrun.c_str());
2903 int qqvvqq_apartordering[2]={ -1, -1 };
2906 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2907 qqvvqq_apartordering
2909 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2910 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2921 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2922 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2923 sprintf(
runstring_.runstring, strrun.c_str());
2927 int qqvvqq_apartordering[2]={ -1, -1 };
2930 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2931 qqvvqq_apartordering
2933 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2934 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2955 else if (useQQVVQQany){
2977 else if (useQQVVQQany){
2993 if (jetid==6) jetid = 2;
3003 else if (useQQVVQQany){
3017 double yy_up = pow(gV_up, 2) + pow(gA_up, 2);
3018 double yy_dn = pow(gV_dn, 2) + pow(gA_dn, 2);
3019 double xx_up = gV_up*gA_up;
3020 double xx_dn = gV_dn*gA_dn;
3021 double yy = (2.*yy_up+3.*yy_dn);
3022 double xx = (2.*xx_up+3.*xx_dn);
3023 double discriminant = pow(yy, 2)-4.*pow(xx, 2);
3024 double gVsq = (yy+sqrt(fabs(discriminant)))/2.;
3025 double gAsq = pow(xx, 2)/gVsq;
3026 double gV=-sqrt(gVsq);
3027 double gA=-sqrt(gAsq);
3037 else if (useQQVVQQany){
3053 else if (useQQVVQQany){
3109 if (jetid==6) jetid = 2;
3129 double yy_up = pow(gV_up, 2) + pow(gA_up, 2);
3130 double yy_dn = pow(gV_dn, 2) + pow(gA_dn, 2);
3131 double xx_up = gV_up*gA_up;
3132 double xx_dn = gV_dn*gA_dn;
3133 double yy = (2.*yy_up+3.*yy_dn);
3134 double xx = (2.*xx_up+3.*xx_dn);
3135 double discriminant = pow(yy, 2)-4.*pow(xx, 2);
3136 double gVsq = (yy+sqrt(fabs(discriminant)))/2.;
3137 double gAsq = pow(xx, 2)/gVsq;
3138 double gV=-sqrt(gVsq);
3139 double gA=-sqrt(gAsq);
3161 else if (useQQVVQQany){
3172 int* ordering=
nullptr;
3200 else if (
process ==
TVar::bkgZZ){ ordering = pZOrder;
if (!hasZ1 || !hasZ2) result =
false; }
3201 else if (
process ==
TVar::bkgWW){ ordering = pWOrder;
if (!hasW1 || !hasW2) result =
false; }
3204 else{ ordering = pZOrder;
if (!hasZ1 || !hasZ2) result =
false; }
3221 if (!hasZ1 || !hasZ2 || !hasW1 || !hasW2) result =
false;
3225 if (!hasW1 || !hasW2) result =
false;
3257 if (!hasZ1 || !hasZ2 || !hasW1 || !hasW2) result =
false;
3261 if (!hasW1 || !hasW2) result =
false;
3267 if (partOrder && ordering) {
for (
unsigned int ip=0; ip<ndau; ip++) partOrder->push_back(ordering[ip]); }
3268 if (apartOrder && pApartOrder) {
for (
unsigned int ip=0; ip<napart; ip++) apartOrder->push_back(pApartOrder[ip]); }
3269 for (
int ip=0; ip<
mxpart; ip++) sprintf((
plabel_.plabel)[ip], strplabel[ip].Data());
3272 MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Summary (result=" << (int)result <<
"):\n";
3276 if (hasZ1)
MELAout <<
"\tProcess found a Z1." << endl;
3277 if (hasZ2)
MELAout <<
"\tProcess found a Z2." << endl;
3278 if (hasW1)
MELAout <<
"\tProcess found a W1." << endl;
3279 if (hasW2)
MELAout <<
"\tProcess found a W2." << endl;
3285 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
3287 if (!partOrder->empty())
MELAout <<
"\tpartOrder[" << partOrder->size() <<
"]:\n";
3288 for (
unsigned int ip=0; ip<partOrder->size(); ip++)
MELAout <<
"\t[" << ip <<
"] -> " << partOrder->at(ip) << endl;
3291 if (!apartOrder->empty())
MELAout <<
"\tapartOrder[" << apartOrder->size() <<
"]:\n";
3292 for (
unsigned int ip=0; ip<apartOrder->size(); ip++)
MELAout <<
"\t[" << ip <<
"] -> " << apartOrder->at(ip) << endl;
3300 if (pApartOrder)
delete[] pApartOrder;
3301 if (pApartId)
delete[] pApartId;
3306 if (useQJ)
return TString(
"qj");
3307 else return TString(
"pp");
3310 else if (pid==1)
return TString(
"dq");
3311 else if (pid==2)
return TString(
"uq");
3312 else if (pid==3)
return TString(
"sq");
3313 else if (pid==4)
return TString(
"cq");
3314 else if (pid==5)
return TString(
"bq");
3315 else if (pid==-1)
return TString(
"da");
3316 else if (pid==-2)
return TString(
"ua");
3317 else if (pid==-3)
return TString(
"sa");
3318 else if (pid==-4)
return TString(
"ca");
3319 else if (pid==-5)
return TString(
"ba");
3320 else if (abs(pid)==6)
return TString(
"qj");
3321 else if (pid==11)
return TString(
"el");
3322 else if (pid==13)
return TString(
"ml");
3323 else if (pid==15)
return TString(
"tl");
3324 else if (pid==-11)
return TString(
"ea");
3325 else if (pid==-13)
return TString(
"ma");
3326 else if (pid==-15)
return TString(
"ta");
3327 else if (std::abs(pid)>=12 && std::abs(pid)<=16){
3328 if (!useExtendedConventions){
3329 if (pid>0)
return TString(
"nl");
3330 else return TString(
"na");
3333 if (pid==12)
return TString(
"ne");
3334 else if (pid==14)
return TString(
"nm");
3335 else if (pid==16)
return TString(
"nt");
3336 else if (pid==-12)
return TString(
"ke");
3337 else if (pid==-14)
return TString(
"km");
3338 else return TString(
"kt");
3341 else return TString(
" ");
3345 double minZmassSqr=10*10;
3349 (s[2][3]<minZmassSqr || s[4][5]<minZmassSqr)
3361 (-s[5-1][1-1]<
cutoff_.cutoff)
3362 || (-s[5-1][2-1]<
cutoff_.cutoff)
3363 || (-s[4-1][1-1]<
cutoff_.cutoff)
3364 || (-s[4-1][2-1]<
cutoff_.cutoff)
3365 || (-s[3-1][1-1]<
cutoff_.cutoff)
3366 || (-s[3-1][2-1]<
cutoff_.cutoff)
3367 || (+s[5-1][4-1]<
cutoff_.cutoff)
3368 || (+s[5-1][3-1]<
cutoff_.cutoff)
3369 || (+s[4-1][3-1]<
cutoff_.cutoff)
3377 (-s[5-1][1-1]<
cutoff_.cutoff)
3378 || (-s[5-1][2-1]<
cutoff_.cutoff)
3379 || (-s[6-1][1-1]<
cutoff_.cutoff)
3380 || (-s[6-1][2-1]<
cutoff_.cutoff)
3381 || (+s[6-1][5-1]<
cutoff_.cutoff)
3393 const double GeV = 1./100.;
3394 collider_sqrts *= GeV;
3396 char path_pdf_c[200];
3397 sprintf(path_pdf_c,
"%s", pathtoPDFSet);
3398 int pathpdfLength = strlen(path_pdf_c);
3403 const double GeV = 1./100.;
3409 int iAllow = (doAllow ? 1 : 0);
3545 for (
int im=0; im<2; im++){
3673 for (
int im=0; im<2; im++){
3783 for (
int im=0; im<2; im++){
3859 for (
int im=0; im<2; im++){
3912 for (
int im=0; im<2; im++){
3967 for (
int im=0; im<2; im++){
4043 for (
int im=0; im<2; im++){
4096 for (
int im=0; im<2; im++){
4136 for (
int im=0; im<2; im++){
4168 for (
int im=0; im<2; im++){
4220 const double GeV = 1./100.;
4221 int iWWcoupl = (useWWcoupl ? 1 : 0);
4247 const double& EBEAM,
4253 int nRequested_AssociatedJets=0;
4254 int nRequested_AssociatedLeptons=0;
4255 int AssociationVCompatibility=0;
4259 nRequested_AssociatedJets = 1;
4260 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::SumMatrixElementPDF: Requesting " << nRequested_AssociatedJets <<
" jets" << endl;
4273 nRequested_AssociatedJets = 2;
4274 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::SumMatrixElementPDF: Requesting " << nRequested_AssociatedJets <<
" jets" << endl;
4282 nRequested_AssociatedLeptons = 2;
4288 ) AssociationVCompatibility=24;
4293 ) AssociationVCompatibility=23;
4306 vector<int> partOrder;
4307 vector<int> apartOrder;
4314 if (partOrder.size()!=mela_event.
pDaughters.size()){
4316 MELAerr <<
"TUtil::SumMatrixElementPDF: Ordering size " << partOrder.size() <<
" and number of daughter particles " << mela_event.
pDaughters.size() <<
" are not the same!" << endl;
4321 if (apartOrder.size()!=mela_event.
pAssociated.size()){
4323 MELAerr <<
"TUtil::SumMatrixElementPDF: Ordering size " << apartOrder.size() <<
" and number of associated particles " << mela_event.
pAssociated.size() <<
" are not the same!" << endl;
4330 int NPart=
npart_.npart+2;
4331 double p4[4][
mxpart]={ { 0 } };
4332 double p4_tmp[4][
mxpart]={ { 0 } };
4333 int id[
mxpart];
for (
int ipar=0; ipar<
mxpart; ipar++)
id[ipar]=-9000;
4335 double msq_tmp[
nmsq][
nmsq]={ { 0 } };
4336 int channeltoggle=0;
4338 TLorentzVector MomStore[
mxpart];
4339 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
4343 for (
int ipar=0; ipar<2; ipar++){
4344 if (mela_event.
pMothers.at(ipar).second.T()>0.){
4345 p4[0][ipar] = -mela_event.
pMothers.at(ipar).second.X();
4346 p4[1][ipar] = -mela_event.
pMothers.at(ipar).second.Y();
4347 p4[2][ipar] = -mela_event.
pMothers.at(ipar).second.Z();
4348 p4[3][ipar] = -mela_event.
pMothers.at(ipar).second.T();
4349 MomStore[ipar] = mela_event.
pMothers.at(ipar).second;
4350 id[ipar] = mela_event.
pMothers.at(ipar).first;
4353 p4[0][ipar] = mela_event.
pMothers.at(ipar).second.X();
4354 p4[1][ipar] = mela_event.
pMothers.at(ipar).second.Y();
4355 p4[2][ipar] = mela_event.
pMothers.at(ipar).second.Z();
4356 p4[3][ipar] = mela_event.
pMothers.at(ipar).second.T();
4357 MomStore[ipar] = -mela_event.
pMothers.at(ipar).second;
4358 id[ipar] = mela_event.
pMothers.at(ipar).first;
4368 for (
int ix=0; ix<(int)partOrder.size(); ix++){
4369 int ipar = min((
int)
mxpart, min(NPart, 2))+ix;
4371 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(partOrder.at(ix)).second);
4372 p4[0][ipar] = momTmp->X();
4373 p4[1][ipar] = momTmp->Y();
4374 p4[2][ipar] = momTmp->Z();
4375 p4[3][ipar] = momTmp->T();
4376 MomStore[ipar]=*momTmp;
4377 id[ipar] = mela_event.
pDaughters.at(partOrder.at(ix)).first;
4379 for (
int ix=0; ix<(int)apartOrder.size(); ix++){
4380 int ipar = min((
int)
mxpart, min(NPart, (
int)(partOrder.size()+2)))+ix;
4382 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(apartOrder.at(ix)).second);
4383 p4[0][ipar] = momTmp->X();
4384 p4[1][ipar] = momTmp->Y();
4385 p4[2][ipar] = momTmp->Z();
4386 p4[3][ipar] = momTmp->T();
4387 MomStore[ipar]=*momTmp;
4388 id[ipar] = mela_event.
pAssociated.at(apartOrder.at(ix)).first;
4391 if (verbosity >=
TVar::DEBUG){
for (
int i=0; i<NPart; i++)
MELAout <<
"p["<<i<<
"] (Px, Py, Pz, E):\t" << p4[0][i] <<
'\t' << p4[1][i] <<
'\t' << p4[2][i] <<
'\t' << p4[3][i] << endl; }
4393 double defaultRenScale =
scale_.scale;
4394 double defaultFacScale =
facscale_.facscale;
4396 int defaultNflav =
nflav_.nflav;
4397 string defaultPdflabel =
pdlabel_.pdlabel;
4401 double alphasVal, alphasmzVal;
4411 <<
"TUtil::SumMatrixElementPDF: Set AlphaS:\n"
4412 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
4414 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
4432 else channeltoggle=2;
4439 for (
int iquark=-5; iquark<=5; iquark++){
4440 for (
int jquark=-5; jquark<=5; jquark++){
4450 ) msqjk =
msq[3][7]+
msq[7][3];
4453 ) msqjk =
msq[jquark+5][-jquark+5];
4456 ) msqjk =
msq[-iquark+5][iquark+5];
4457 else msqjk +=
msq[jquark+5][iquark+5];
4461 SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4474 MELAout <<
"\tTUtil::SumMatrixElementPDF: ZZGG && ZZ/WW/ZZ+WW MEs using ZZ (runstring: " <<
runstring_.runstring <<
")" << endl;
4475 for (
int i=0; i<NPart; i++)
MELAout <<
"\tp["<<i<<
"] (Px, Py, Pz, E):\t" << p4[0][i] <<
'\t' << p4[1][i] <<
'\t' << p4[2][i] <<
'\t' << p4[3][i] << endl;
4483 MELAout <<
"\tTUtil::SumMatrixElementPDF: ZZGG && ZZ/WW/ZZ+WW MEs using WW (runstring: " <<
runstring_.runstring <<
")" << endl;
4484 for (
int i=0; i<NPart; i++)
MELAout <<
"\tp["<<i<<
"] (Px, Py, Pz, E):\t" << p4[0][i] <<
'\t' << p4[1][i] <<
'\t' << p4[2][i] <<
'\t' << p4[3][i] << endl;
4490 for (
int iquark=-5; iquark<=5; iquark++){
4491 for (
int jquark=-5; jquark<=5; jquark++){
4492 if (iquark==0 && jquark==0)
continue;
4493 msq[jquark+5][iquark+5] = 0;
4507 SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4552 for (
unsigned int ix=0; ix<4; ix++){
4553 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4554 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4557 for (
int iquark=-5; iquark<=5; iquark++){
for (
int jquark=-5; jquark<=5; jquark++){
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] + msq_tmp[jquark+5][iquark+5]);
if (iquark==jquark)
msq[jquark+5][iquark+5]*=0.5; } }
4560 for (
unsigned int ix=0; ix<4; ix++){
4561 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4562 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4564 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4565 for (
int ip=0; ip<
mxpart; ip++){
4566 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4567 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4571 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4573 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4574 MELAout <<
"\tMEsq initial:" << endl;
4575 for (
int iquark=-5; iquark<=5; iquark++){
4576 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4579 MELAout <<
"\tMEsq added:" << endl;
4580 for (
int iquark=-5; iquark<=5; iquark++){
4581 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4585 for (
int iquark=-5; iquark<=5; iquark++){
4586 for (
int jquark=-5; jquark<=5; jquark++){
4587 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4598 for (
unsigned int ix=0; ix<4; ix++){
4599 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4600 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4603 for (
int iquark=-5; iquark<=5; iquark++){
for (
int jquark=-5; jquark<=5; jquark++){
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] + msq_tmp[jquark+5][iquark+5]);
if (iquark==jquark && iquark!=0)
msq[jquark+5][iquark+5]*=0.5; } }
4606 for (
int ip=0; ip<
mxpart; ip++){
4607 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4608 else sprintf((
plabel_.plabel)[ip], gglabel.Data());
4611 for (
int iquark=-5; iquark<=5; iquark++){
int jquark=-iquark;
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] - msq_tmp[jquark+5][iquark+5]); }
4613 MELAout <<
"TUtil::SumMatrixElementPDF: Subtracting qqb/qbq->gg double-counted contribution:\n";
4615 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4617 for (
int iquark=-5; iquark<=5; iquark++){
4618 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4624 for (
unsigned int ix=0; ix<4; ix++){
4625 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4626 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4628 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4629 for (
int ip=0; ip<
mxpart; ip++){
4630 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4631 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4635 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4637 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4638 MELAout <<
"\tMEsq initial:" << endl;
4639 for (
int iquark=-5; iquark<=5; iquark++){
4640 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4643 MELAout <<
"\tMEsq added:" << endl;
4644 for (
int iquark=-5; iquark<=5; iquark++){
4645 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4649 for (
int iquark=-5; iquark<=5; iquark++){
4650 for (
int jquark=-5; jquark<=5; jquark++){
4651 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4671 for (
unsigned int ix=0; ix<4; ix++){
4672 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4673 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4676 for (
int iquark=-5; iquark<=5; iquark++){
for (
int jquark=-5; jquark<=5; jquark++){
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] + msq_tmp[jquark+5][iquark+5]);
if (iquark==jquark)
msq[jquark+5][iquark+5]*=0.5; } }
4679 for (
unsigned int ix=0; ix<4; ix++){
4680 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4681 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4683 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4684 for (
int ip=0; ip<
mxpart; ip++){
4685 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4686 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4690 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4692 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4693 MELAout <<
"\tMEsq initial:" << endl;
4694 for (
int iquark=-5; iquark<=5; iquark++){
4695 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4698 MELAout <<
"\tMEsq added:" << endl;
4699 for (
int iquark=-5; iquark<=5; iquark++){
4700 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4704 for (
int iquark=-5; iquark<=5; iquark++){
4705 for (
int jquark=-5; jquark<=5; jquark++){
4706 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4716 for (
unsigned int ix=0; ix<4; ix++){
4717 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4718 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4721 for (
int iquark=-5; iquark<=5; iquark++){
for (
int jquark=-5; jquark<=5; jquark++){
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] + msq_tmp[jquark+5][iquark+5]);
if (iquark==jquark && iquark!=0)
msq[jquark+5][iquark+5]*=0.5; } }
4724 for (
int ip=0; ip<
mxpart; ip++){
4725 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4726 else sprintf((
plabel_.plabel)[ip], gglabel.Data());
4729 for (
int iquark=-5; iquark<=5; iquark++){
int jquark=-iquark;
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] - msq_tmp[jquark+5][iquark+5]); }
4731 MELAout <<
"TUtil::SumMatrixElementPDF: Subtracting qqb/qbq->gg double-counted contribution:\n";
4733 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4735 for (
int iquark=-5; iquark<=5; iquark++){
4736 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4742 for (
unsigned int ix=0; ix<4; ix++){
4743 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4744 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4746 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4747 for (
int ip=0; ip<
mxpart; ip++){
4748 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4749 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4753 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4755 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4756 MELAout <<
"\tMEsq initial:" << endl;
4757 for (
int iquark=-5; iquark<=5; iquark++){
4758 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4761 MELAout <<
"\tMEsq added:" << endl;
4762 for (
int iquark=-5; iquark<=5; iquark++){
4763 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4767 for (
int iquark=-5; iquark<=5; iquark++){
4768 for (
int jquark=-5; jquark<=5; jquark++){
4769 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4789 for (
unsigned int ix=0; ix<4; ix++){
4790 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4791 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4794 for (
int iquark=-5; iquark<=5; iquark++){
for (
int jquark=-5; jquark<=5; jquark++){
msq[jquark+5][iquark+5] = (
msq[jquark+5][iquark+5] + msq_tmp[jquark+5][iquark+5]);
if (iquark==jquark)
msq[jquark+5][iquark+5]*=0.5; } }
4797 for (
unsigned int ix=0; ix<4; ix++){
4798 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4799 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4801 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4802 for (
int ip=0; ip<
mxpart; ip++){
4803 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4804 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4808 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4810 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4811 MELAout <<
"\tMEsq initial:" << endl;
4812 for (
int iquark=-5; iquark<=5; iquark++){
4813 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4816 MELAout <<
"\tMEsq added:" << endl;
4817 for (
int iquark=-5; iquark<=5; iquark++){
4818 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4822 for (
int iquark=-5; iquark<=5; iquark++){
4823 for (
int jquark=-5; jquark<=5; jquark++){
4824 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4835 msqjk =
SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4852 if (msqjk != msqjk){
4854 for (
int i=0; i<NPart; i++)
MELAout <<
"p["<<i<<
"] (Px, Py, Pz, E):\t" << p4[0][i] <<
'\t' << p4[1][i] <<
'\t' << p4[2][i] <<
'\t' << p4[3][i] << endl;
4855 MELAout <<
"TUtil::SumMatrixElementPDF: The number of ME instances: " << nInstances <<
". MEsq[ip][jp] = " << endl;
4856 for (
int iquark=-5; iquark<=5; iquark++){
4857 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4863 MELAout <<
"TUtil::SumMatrixElementPDF: The number of ME instances: " << nInstances <<
". MEsq[ip][jp] = " << endl;
4864 for (
int iquark=-5; iquark<=5; iquark++){
4865 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4871 <<
"TUtil::SumMatrixElementPDF: Reset AlphaS:\n"
4872 <<
"\tBefore reset, alphas scale: " <<
scale_.scale
4873 <<
", PDF scale: " <<
facscale_.facscale
4875 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
4877 <<
"TUtil::SumMatrixElementPDF: Reset AlphaS result:\n"
4878 <<
"\tAfter reset, alphas scale: " <<
scale_.scale
4879 <<
", PDF scale: " <<
facscale_.facscale
4882 if (verbosity>=
TVar::DEBUG)
MELAout <<
"End TUtil::SumMatrixElementPDF(" << msqjk <<
")" << endl;
4890 const double& EBEAM,
4893 const double GeV=1./100.;
4896 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::JHUGenMatEl: Non-JHUGen MEs are not supported" << endl;
return MatElSq; }
4931 if (isSpinZero)
MELAout <<
"0";
4932 else if (isSpinOne)
MELAout <<
"1";
4933 else if (isSpinTwo)
MELAout <<
"2";
4939 int MYIDUP_tmp[4]={ 0 };
4940 vector<pair<int, int>> idarray[2];
4941 int MYIDUP[4]={ 0 };
4942 int idfirst[2]={ 0 };
4943 int idsecond[2]={ 0 };
4944 double p4[6][4]={ { 0 } };
4945 TLorentzVector MomStore[
mxpart];
4946 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
4966 for (
int ipar=0; ipar<2; ipar++){
4967 if (mela_event.
pMothers.at(ipar).second.T()>0.){
4968 p4[ipar][0] = -mela_event.
pMothers.at(ipar).second.T()*GeV;
4969 p4[ipar][1] = -mela_event.
pMothers.at(ipar).second.X()*GeV;
4970 p4[ipar][2] = -mela_event.
pMothers.at(ipar).second.Y()*GeV;
4971 p4[ipar][3] = -mela_event.
pMothers.at(ipar).second.Z()*GeV;
4972 MomStore[ipar] = mela_event.
pMothers.at(ipar).second;
4975 p4[ipar][0] = mela_event.
pMothers.at(ipar).second.T()*GeV;
4976 p4[ipar][1] = mela_event.
pMothers.at(ipar).second.X()*GeV;
4977 p4[ipar][2] = mela_event.
pMothers.at(ipar).second.Y()*GeV;
4978 p4[ipar][3] = mela_event.
pMothers.at(ipar).second.Z()*GeV;
4979 MomStore[ipar] = -mela_event.
pMothers.at(ipar).second;
4983 if (ipar==1 && production ==
TVar::ZZINDEPENDENT){
for (
int ix=0; ix<4; ix++){ p4[0][ix] += p4[ipar][ix]; p4[ipar][ix]=0.; } }
4987 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
4988 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
4989 int* idtmp = &(mela_event.
pDaughters.at(ipar).first);
4991 int arrindex = 2*ipar;
4993 else MYIDUP_tmp[arrindex] = 0;
4994 MYIDUP_tmp[arrindex+1] = -9000;
4995 p4[arrindex+2][0] = momTmp->T()*GeV;
4996 p4[arrindex+2][1] = momTmp->X()*GeV;
4997 p4[arrindex+2][2] = momTmp->Y()*GeV;
4998 p4[arrindex+2][3] = momTmp->Z()*GeV;
4999 MomStore[arrindex+2] = *momTmp;
5000 if (verbosity >=
TVar::DEBUG)
MELAout <<
"MYIDUP_tmp[" << arrindex <<
"(" << ipar <<
")" <<
"]=" << MYIDUP_tmp[arrindex] << endl;
5004 for (
unsigned int ipar=0; ipar<4; ipar++){
5006 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
5007 int* idtmp = &(mela_event.
pDaughters.at(ipar).first);
5010 else MYIDUP_tmp[ipar] = 0;
5011 p4[ipar+2][0] = momTmp->T()*GeV;
5012 p4[ipar+2][1] = momTmp->X()*GeV;
5013 p4[ipar+2][2] = momTmp->Y()*GeV;
5014 p4[ipar+2][3] = momTmp->Z()*GeV;
5015 MomStore[ipar+2] = *momTmp;
5017 else MYIDUP_tmp[ipar] = -9000;
5019 if (verbosity >=
TVar::DEBUG)
MELAout <<
"MYIDUP_tmp[" << ipar <<
"]=" << MYIDUP_tmp[ipar] << endl;
5024 double defaultRenScale =
scale_.scale;
5025 double defaultFacScale =
facscale_.facscale;
5027 int defaultNflav =
nflav_.nflav;
5028 string defaultPdflabel =
pdlabel_.pdlabel;
5032 double alphasVal, alphasmzVal;
5042 <<
"TUtil::JHUGenMatEl: Set AlphaS:\n"
5043 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
5045 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5049 for (
int iv=0; iv<2; iv++){
5050 if (MYIDUP_tmp[2*iv+0]!=0 && MYIDUP_tmp[2*iv+1]!=0){
5053 TMath::Sign(1, MYIDUP_tmp[2*iv+0])==TMath::Sign(1, -MYIDUP_tmp[2*iv+1])
5055 (MYIDUP_tmp[2*iv+0]==-9000 || MYIDUP_tmp[2*iv+1]==-9000)
5056 ) idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], MYIDUP_tmp[2*iv+1]));
5058 else if (MYIDUP_tmp[2*iv+0]<0) idarray[iv].push_back(pair<int, int>(-MYIDUP_tmp[2*iv+0], MYIDUP_tmp[2*iv+1]));
5059 else idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], -MYIDUP_tmp[2*iv+1]));
5061 else if (MYIDUP_tmp[2*iv+0]!=0){
5064 int id_dn = TMath::Sign(1, -MYIDUP_tmp[2*iv+0]);
5065 int id_st = TMath::Sign(3, -MYIDUP_tmp[2*iv+0]);
5066 int id_bt = TMath::Sign(5, -MYIDUP_tmp[2*iv+0]);
5067 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_dn));
5068 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_st));
5069 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_bt));
5072 int id_up = TMath::Sign(2, -MYIDUP_tmp[2*iv+0]);
5073 int id_ch = TMath::Sign(4, -MYIDUP_tmp[2*iv+0]);
5074 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_up));
5075 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_ch));
5079 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], -MYIDUP_tmp[2*iv+0]));
5082 else if (MYIDUP_tmp[2*iv+1]!=0){
5085 int id_dn = TMath::Sign(1, -MYIDUP_tmp[2*iv+1]);
5086 int id_st = TMath::Sign(3, -MYIDUP_tmp[2*iv+1]);
5087 int id_bt = TMath::Sign(5, -MYIDUP_tmp[2*iv+1]);
5088 idarray[iv].push_back(pair<int, int>(id_dn, MYIDUP_tmp[2*iv+1]));
5089 idarray[iv].push_back(pair<int, int>(id_st, MYIDUP_tmp[2*iv+1]));
5090 idarray[iv].push_back(pair<int, int>(id_bt, MYIDUP_tmp[2*iv+1]));
5093 int id_up = TMath::Sign(2, -MYIDUP_tmp[2*iv+1]);
5094 int id_ch = TMath::Sign(4, -MYIDUP_tmp[2*iv+1]);
5095 idarray[iv].push_back(pair<int, int>(id_up, MYIDUP_tmp[2*iv+1]));
5096 idarray[iv].push_back(pair<int, int>(id_ch, MYIDUP_tmp[2*iv+1]));
5100 idarray[iv].push_back(pair<int, int>(-MYIDUP_tmp[2*iv+1], MYIDUP_tmp[2*iv+1]));
5105 for (
int iqu=1; iqu<=2; iqu++){
5107 for (
int iqd=1; iqd<=3; iqd++){
5108 int id_dn = iqd*2-1;
5109 if (iv==0){ idarray[iv].push_back(pair<int, int>(id_up, -id_dn)); idarray[iv].push_back(pair<int, int>(-id_dn, id_up)); }
5110 else{ idarray[iv].push_back(pair<int, int>(id_dn, -id_up)); idarray[iv].push_back(pair<int, int>(-id_up, id_dn)); }
5115 for (
int iq=1; iq<=5; iq++){ idarray[iv].push_back(pair<int, int>(iq, -iq)); idarray[iv].push_back(pair<int, int>(-iq, iq)); }
5121 const int InvalidMode=-1, WWMode=0, ZZMode=1, ZgMode=5, ggMode=7;
5122 int VVMode=InvalidMode;
5131 double aL1=0, aL2=0, aR1=0, aR2=0;
5134 for (
unsigned int v1=0; v1<idarray[0].size(); v1++){
5135 for (
unsigned int v2=0; v2<idarray[1].size(); v2++){
5137 if (idarray[0].at(v1).first!=-9000) MYIDUP[0] =
convertLHEreverse(&(idarray[0].at(v1).first));
5139 if (idarray[1].at(v2).first!=-9000) MYIDUP[2] =
convertLHEreverse(&(idarray[1].at(v2).first));
5143 if (verbosity>=
TVar::DEBUG){
for (
unsigned int idau=0; idau<4; idau++)
MELAout <<
"MYIDUP[" << idau <<
"]=" << MYIDUP[idau] << endl; }
5146 for (
int ip=0; ip<2; ip++){ idfirst[ip]=MYIDUP[ip]; idsecond[ip]=MYIDUP[ip+2]; }
5151 MELAout <<
"TUtil::JHUGenMatEl: M_V=" << mv/GeV <<
", Ga_V=" << gv/GeV << endl;
5153 MELAout <<
"TUtil::JHUGenMatEl: M_Vprime=" << mv/GeV <<
", Ga_Vprime=" << gv/GeV << endl;
5157 double aLRtmp[4]={ 0 };
5159 if (idarray[0].size()>1){
5160 aL1 = sqrt(pow(aL1, 2)+pow(aLRtmp[0], 2));
5161 aR1 = sqrt(pow(aR1, 2)+pow(aLRtmp[1], 2));
5167 if (idarray[1].size()>1){
5168 aL2 = sqrt(pow(aL2, 2)+pow(aLRtmp[2], 2));
5169 aR2 = sqrt(pow(aR2, 2)+pow(aLRtmp[3], 2));
5194 if (verbosity >=
TVar::DEBUG)
MELAout <<
"=====\nTUtil::JHUGenMatEl: Instance MatElTmp = " << MatElTmp <<
"\n=====" << endl;
5195 MatElSq += MatElTmp;
5196 if (MatElTmp>0.) nNonZero++;
5199 if (nNonZero>0) MatElSq /= ((double)nNonZero);
5201 MELAout <<
"TUtil::JHUGenMatEl: Number of matrix element instances computed: " << nNonZero << endl;
5202 MELAout <<
"TUtil::JHUGenMatEl: MatElSq after division = " << MatElSq << endl;
5209 <<
"TUtil::JHUGenMatEl: aL1, aR1, aL2, aR2: "
5210 << aL1 <<
", " << aR1 <<
", " << aL2 <<
", " << aR2
5214 int GeVexponent_MEsq = 4-((int)mela_event.
pDaughters.size())*2;
5216 double constant = pow(GeV, -GeVexponent_MEsq);
5217 MatElSq *= constant;
5220 if (mela_event.
pMothers.at(0).first!=0 && mela_event.
pMothers.at(1).first!=0){
5222 if (abs(mela_event.
pMothers.at(0).first)<6) ix=mela_event.
pMothers.at(0).first + 5;
5223 if (abs(mela_event.
pMothers.at(1).first)<6) iy=mela_event.
pMothers.at(1).first + 5;
5224 msq[iy][ix]=MatElSq;
5229 for (
int ix=0; ix<5; ix++){
msq[ix][10-ix]=MatElSq;
msq[10-ix][ix]=MatElSq; }
5234 double fx_dummy[
nmsq]={ 0 }; fx_dummy[5]=1.;
5240 if (verbosity >=
TVar::DEBUG)
MELAout <<
"TUtil::JHUGenMatEl: Final MatElSq = " << MatElSq << endl;
5245 <<
"TUtil::JHUGenMatEl: Reset AlphaS:\n"
5246 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
5248 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
5252 <<
"TUtil::JHUGenMatEl: Reset AlphaS result:\n"
5253 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5261 const double& EBEAM,
5265 if (matrixElement!=
TVar::MADGRAPH){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::MadgraphMatEl: Non-Madgraph MEs are not supported" << endl;
return MatElSq; }
5276 TLorentzVector MomStore[
mxpart];
5277 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
5283 double* p =
new double[4*nPDG];
5287 pdgs[i] = particle.first;
5288 p[i*4+0] = particle.second.E();
5289 p[i*4+1] = particle.second.Px();
5290 p[i*4+2] = particle.second.Py();
5291 p[i*4+3] = particle.second.Pz();
5292 MomStore[i] = particle.second;
5296 bool previously_swapped =
false;
5298 bool swap_spaces =
false;
5299 if((i == 2 || i == 4) && particle.first > 0){
5302 if(verbosity >=
TVar::DEBUG)
MELAout <<
"Swapping daughters at index " << i <<
" and " << i-1 << endl;
5304 pdgs[i] = particle.first;
5305 p[i*4+0] = particle.second.E();
5306 p[i*4+1] = particle.second.Px();
5307 p[i*4+2] = particle.second.Py();
5308 p[i*4+3] = particle.second.Pz();
5309 MomStore[i] = particle.second;
5310 if(previously_swapped){
5312 previously_swapped =
false;
5314 else if(swap_spaces){
5316 previously_swapped =
true;
5321 if(abs(pdgs[2]) > abs(pdgs[4])){
5322 swap(pdgs[2], pdgs[4]);
5323 swap(pdgs[3], pdgs[5]);
5325 swap(p[2*4+i], p[4*4+i]);
5326 swap(p[3*4+i], p[5*4+i]);
5330 double defaultRenScale =
scale_.scale;
5331 double defaultFacScale =
facscale_.facscale;
5333 int defaultNflav =
nflav_.nflav;
5334 string defaultPdflabel =
pdlabel_.pdlabel;
5338 double alphasVal, alphasmzVal;
5348 <<
"TUtil::MadgraphMatEl: Set AlphaS:\n"
5349 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
5351 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5360 MELAout <<
"Raw Input to FORTRAN:" << endl;
5361 for(
int i = 0; i < nPDG; i++){
5362 for(
int j = 0; j < 4; j++){
5372 pdgs, procid, nPDG, p,
5373 alphasVal, scale2, nhel, MatElSq
5376 if(verbosity >=
TVar::DEBUG)
MELAout <<
" smatrixhel returns prob of " << MatElSq << endl;
5386 const double& EBEAM,
5389 const double GeV=1./100.;
5390 double sum_msqjk = 0;
5398 double MatElsq[
nmsq][
nmsq]={ { 0 } };
5399 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
5402 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::HJJMatEl: Non-JHUGen MEs are not supported" << endl;
return sum_msqjk; }
5406 int nRequested_AssociatedJets=2;
5407 if (production ==
TVar::JQCD) nRequested_AssociatedJets=1;
5417 if (mela_event.
pAssociated.size()==0){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::HJJMatEl: Number of associated particles is 0!" << endl;
return sum_msqjk; }
5419 int MYIDUP_tmp[4]={ 0 };
5420 double p4[5][4]={ { 0 } };
5421 double pOneJet[4][4] ={ { 0 } };
5422 TLorentzVector MomStore[
mxpart];
5423 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
5429 for (
int ipar=0; ipar<2; ipar++){
5430 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
5431 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
5433 else MYIDUP_tmp[ipar] = 0;
5434 if (momTmp->T()>0.){
5435 p4[ipar][0] = momTmp->T()*GeV;
5436 p4[ipar][1] = momTmp->X()*GeV;
5437 p4[ipar][2] = momTmp->Y()*GeV;
5438 p4[ipar][3] = momTmp->Z()*GeV;
5439 MomStore[ipar] = (*momTmp);
5442 p4[ipar][0] = -momTmp->T()*GeV;
5443 p4[ipar][1] = -momTmp->X()*GeV;
5444 p4[ipar][2] = -momTmp->Y()*GeV;
5445 p4[ipar][3] = -momTmp->Z()*GeV;
5446 MomStore[ipar] = -(*momTmp);
5447 MYIDUP_tmp[ipar] = -MYIDUP_tmp[ipar];
5450 for (
unsigned int ipar=0; ipar<2; ipar++){
5452 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(ipar).second);
5453 int* idtmp = &(mela_event.
pAssociated.at(ipar).first);
5455 else MYIDUP_tmp[ipar+2] = 0;
5456 p4[ipar+2][0] = momTmp->T()*GeV;
5457 p4[ipar+2][1] = momTmp->X()*GeV;
5458 p4[ipar+2][2] = momTmp->Y()*GeV;
5459 p4[ipar+2][3] = momTmp->Z()*GeV;
5460 MomStore[ipar+6] = (*momTmp);
5462 else MYIDUP_tmp[ipar+2] = -9000;
5464 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
5465 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
5466 p4[4][0] += momTmp->T()*GeV;
5467 p4[4][1] += momTmp->X()*GeV;
5468 p4[4][2] += momTmp->Y()*GeV;
5469 p4[4][3] += momTmp->Z()*GeV;
5470 MomStore[5] = MomStore[5] + (*momTmp);
5473 for (
unsigned int i = 0; i < 4; i++){
5474 if (i<2){
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[i][j]; }
5475 else if (i==2){
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[4][j]; }
5476 else{
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[2][j]; }
5479 for (
unsigned int i=0; i<5; i++)
MELAout <<
"p["<<i<<
"] (Px, Py, Pz, E, M):\t" << p4[i][1]/GeV <<
'\t' << p4[i][2]/GeV <<
'\t' << p4[i][3]/GeV <<
'\t' << p4[i][0]/GeV <<
'\t' << sqrt(fabs(pow(p4[i][0], 2)-pow(p4[i][1], 2)-pow(p4[i][2], 2)-pow(p4[i][3], 2)))/GeV << endl;
5480 MELAout <<
"One-jet momenta: " << endl;
5481 for (
unsigned int i=0; i<4; i++)
MELAout <<
"pOneJet["<<i<<
"] (Px, Py, Pz, E, M):\t" << pOneJet[i][1]/GeV <<
'\t' << pOneJet[i][2]/GeV <<
'\t' << pOneJet[i][3]/GeV <<
'\t' << pOneJet[i][0]/GeV <<
'\t' << sqrt(fabs(pow(pOneJet[i][0], 2)-pow(pOneJet[i][1], 2)-pow(pOneJet[i][2], 2)-pow(pOneJet[i][3], 2)))/GeV << endl;
5484 double defaultRenScale =
scale_.scale;
5485 double defaultFacScale =
facscale_.facscale;
5487 int defaultNflav =
nflav_.nflav;
5488 string defaultPdflabel =
pdlabel_.pdlabel;
5492 double alphasVal, alphasmzVal;
5502 <<
"TUtil::HJJMatEl: Set AlphaS:\n"
5503 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
5505 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5509 bool partonIsUnknown[4];
5510 for (
unsigned int ip=0; ip<4; ip++) partonIsUnknown[ip] = (MYIDUP_tmp[ip]==0);
5516 for (
int isel=-5; isel<=5; isel++){
5517 if (!partonIsUnknown[0] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
continue;
5518 for (
int jsel=-5; jsel<=5; jsel++){
5519 if (!partonIsUnknown[1] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
continue;
5521 if (isel!=0 && jsel!=0) rsel=0;
5522 else if (isel==0) rsel=jsel;
5524 if (!partonIsUnknown[2] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
continue;
5525 MatElsq[jsel+5][isel+5] = MatElsq_tmp[jsel+5][isel+5];
5526 if (verbosity >=
TVar::DEBUG)
MELAout <<
"Channel (isel, jsel)=" << isel <<
", " << jsel << endl;
5532 int nijchannels = ijsel.size();
5533 for (
int ic=0; ic<nijchannels; ic++){
5535 int isel = ijsel[ic][0];
5536 int jsel = ijsel[ic][1];
5537 int code = ijsel[ic][2];
5547 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
5549 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
5552 if (isel==0 && jsel==0){
5563 MatElsq[jsel+5][isel+5] += msq_tmp;
5564 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5572 MatElsq[jsel+5][isel+5] += msq_tmp;
5573 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5584 MatElsq[jsel+5][isel+5] += msq_tmp;
5585 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5589 else if (isel==0 || jsel==0){
5591 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
5593 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && ssel==0) || MYIDUP_tmp[3]==ssel))
5596 MatElsq[jsel+5][isel+5] += msq_tmp;
5597 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5600 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && ssel==0) || MYIDUP_tmp[2]==ssel))
5602 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && rsel==0) || MYIDUP_tmp[3]==rsel))
5605 MatElsq[jsel+5][isel+5] += msq_tmp;
5606 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5609 else if ((isel>0 && jsel<0) || (isel<0 && jsel>0)){
5610 if (code==1 && isel==-jsel){
5618 MatElsq[jsel+5][isel+5] += msq_tmp;
5619 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5622 else if (code==3 && isel==-jsel){
5623 if (abs(isel)!=1){ rsel=1; ssel=-1; }
5624 else{ rsel=2; ssel=-2; }
5627 (partonIsUnknown[2] && partonIsUnknown[3])
5629 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]<0)
5631 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]>0 && partonIsUnknown[3])
5633 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]>0)
5636 MatElsq[jsel+5][isel+5] += msq_tmp;
5637 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5640 (partonIsUnknown[2] && partonIsUnknown[3])
5642 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]>0)
5644 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]<0 && partonIsUnknown[3])
5646 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]<0)
5649 MatElsq[jsel+5][isel+5] += msq_tmp;
5650 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5655 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5657 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5660 MatElsq[jsel+5][isel+5] += msq_tmp;
5661 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5664 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5666 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5669 MatElsq[jsel+5][isel+5] += msq_tmp;
5670 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5676 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5678 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5681 MatElsq[jsel+5][isel+5] += msq_tmp;
5682 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5687 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5689 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5692 MatElsq[jsel+5][isel+5] += msq_tmp;
5693 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5698 if (isel==jsel)
continue;
5699 isel = ijsel[ic][1];
5700 jsel = ijsel[ic][0];
5708 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
5710 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
5713 if (isel==0 || jsel==0){
5715 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
5717 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && ssel==0) || MYIDUP_tmp[3]==ssel))
5720 MatElsq[jsel+5][isel+5] += msq_tmp;
5721 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5724 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && ssel==0) || MYIDUP_tmp[2]==ssel))
5726 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && rsel==0) || MYIDUP_tmp[3]==rsel))
5729 MatElsq[jsel+5][isel+5] += msq_tmp;
5730 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5733 else if ((isel>0 && jsel<0) || (isel<0 && jsel>0)){
5734 if (code==1 && isel==-jsel){
5742 MatElsq[jsel+5][isel+5] += msq_tmp;
5743 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5746 else if(code==3 && isel==-jsel){
5747 if (abs(isel)!=1){ rsel=1; ssel=-1; }
5748 else{ rsel=2; ssel=-2; }
5751 (partonIsUnknown[2] && partonIsUnknown[3])
5753 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]<0)
5755 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]>0 && partonIsUnknown[3])
5757 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]>0)
5760 MatElsq[jsel+5][isel+5] += msq_tmp;
5761 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5764 (partonIsUnknown[2] && partonIsUnknown[3])
5766 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]>0)
5768 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]<0 && partonIsUnknown[3])
5770 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]<0)
5773 MatElsq[jsel+5][isel+5] += msq_tmp;
5774 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5779 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5781 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5784 MatElsq[jsel+5][isel+5] += msq_tmp;
5785 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5788 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5790 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5793 MatElsq[jsel+5][isel+5] += msq_tmp;
5794 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5800 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5802 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5805 MatElsq[jsel+5][isel+5] += msq_tmp;
5806 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5811 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5813 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5816 MatElsq[jsel+5][isel+5] += msq_tmp;
5817 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5824 int isel, jsel, rsel, ssel;
5833 double ckm_takeout=1;
5834 double msq_uu_zz_ijrs1234[2]={ 0 };
5835 double msq_uu_zz_ijrs1243[2]={ 0 };
5836 double msq_dd_zz_ijrs1234[2]={ 0 };
5837 double msq_dd_zz_ijrs1243[2]={ 0 };
5838 double msq_ubarubar_zz_ijrs1234[2]={ 0 };
5839 double msq_ubarubar_zz_ijrs1243[2]={ 0 };
5840 double msq_dbardbar_zz_ijrs1234[2]={ 0 };
5841 double msq_dbardbar_zz_ijrs1243[2]={ 0 };
5842 double msq_uu_zzid_ijrs1234[2]={ 0 };
5843 double msq_uu_zzid_ijrs1243[2]={ 0 };
5844 double msq_dd_zzid_ijrs1234[2]={ 0 };
5845 double msq_dd_zzid_ijrs1243[2]={ 0 };
5846 double msq_ubarubar_zzid_ijrs1234[2]={ 0 };
5847 double msq_ubarubar_zzid_ijrs1243[2]={ 0 };
5848 double msq_dbardbar_zzid_ijrs1234[2]={ 0 };
5849 double msq_dbardbar_zzid_ijrs1243[2]={ 0 };
5850 double msq_udbar_zz_ijrs1234[2]={ 0 };
5851 double msq_udbar_zz_ijrs1243[2]={ 0 };
5852 double msq_dubar_zz_ijrs1234[2]={ 0 };
5853 double msq_dubar_zz_ijrs1243[2]={ 0 };
5855 double msq_uubar_zz_ijrs1234[2]={ 0 };
5856 double msq_uubar_zz_ijrs1243[2]={ 0 };
5857 double msq_ddbar_zz_ijrs1234[2]={ 0 };
5858 double msq_ddbar_zz_ijrs1243[2]={ 0 };
5859 double msq_uubar_ww_ijrs1234[2]={ 0 };
5860 double msq_uubar_ww_ijrs1243[2]={ 0 };
5861 double msq_ddbar_ww_ijrs1234[2]={ 0 };
5862 double msq_ddbar_ww_ijrs1243[2]={ 0 };
5864 double msq_ud_wwonly_ijrs1234[2]={ 0 };
5865 double msq_ud_wwonly_ijrs1243[2]={ 0 };
5866 double msq_ubardbar_wwonly_ijrs1234[2]={ 0 };
5867 double msq_ubardbar_wwonly_ijrs1243[2]={ 0 };
5881 rsel=isel; ssel=jsel;
5885 msq_uu_zz_ijrs1234[1] = msq_uu_zz_ijrs1234[0];
5886 msq_uu_zz_ijrs1243[1] = msq_uu_zz_ijrs1243[0];
5889 rsel=isel; ssel=jsel;
5893 msq_uu_zzid_ijrs1234[1] = msq_uu_zzid_ijrs1234[0];
5894 msq_uu_zzid_ijrs1243[1] = msq_uu_zzid_ijrs1243[0];
5907 rsel=isel; ssel=jsel;
5911 msq_dd_zz_ijrs1234[1] = msq_dd_zz_ijrs1234[0];
5912 msq_dd_zz_ijrs1243[1] = msq_dd_zz_ijrs1243[0];
5915 rsel=isel; ssel=jsel;
5919 msq_dd_zzid_ijrs1234[1] = msq_dd_zzid_ijrs1234[0];
5920 msq_dd_zzid_ijrs1243[1] = msq_dd_zzid_ijrs1243[0];
5933 rsel=isel; ssel=jsel;
5937 msq_ubarubar_zz_ijrs1234[1] = msq_ubarubar_zz_ijrs1234[0];
5938 msq_ubarubar_zz_ijrs1243[1] = msq_ubarubar_zz_ijrs1243[0];
5941 rsel=isel; ssel=jsel;
5945 msq_ubarubar_zzid_ijrs1234[1] = msq_ubarubar_zzid_ijrs1234[0];
5946 msq_ubarubar_zzid_ijrs1243[1] = msq_ubarubar_zzid_ijrs1243[0];
5959 rsel=isel; ssel=jsel;
5963 msq_dbardbar_zz_ijrs1234[1] = msq_dbardbar_zz_ijrs1234[0];
5964 msq_dbardbar_zz_ijrs1243[1] = msq_dbardbar_zz_ijrs1243[0];
5967 rsel=isel; ssel=jsel;
5971 msq_dbardbar_zzid_ijrs1234[1] = msq_dbardbar_zzid_ijrs1234[0];
5972 msq_dbardbar_zzid_ijrs1243[1] = msq_dbardbar_zzid_ijrs1243[0];
5977 (partonIsUnknown[0] && partonIsUnknown[1])
5991 (partonIsUnknown[2] && partonIsUnknown[3])
6005 rsel=isel; ssel=jsel;
6014 (partonIsUnknown[0] && partonIsUnknown[1])
6028 (partonIsUnknown[2] && partonIsUnknown[3])
6042 rsel=isel; ssel=jsel;
6050 (partonIsUnknown[0] && partonIsUnknown[1])
6060 (partonIsUnknown[2] && partonIsUnknown[3])
6068 rsel=isel; ssel=jsel;
6075 (partonIsUnknown[2] && partonIsUnknown[3])
6094 if (abs(isel)>=6)
break;
6103 ckm_takeout = pow(ckm_ir*ckm_js, 2);
6104 for (
unsigned int iswap=0; iswap<2; iswap++){
6105 msq_uubar_ww_ijrs1234[iswap] /= ckm_takeout;
6106 msq_uubar_ww_ijrs1243[iswap] /= ckm_takeout;
6113 (partonIsUnknown[0] && partonIsUnknown[1])
6123 (partonIsUnknown[2] && partonIsUnknown[3])
6131 rsel=isel; ssel=jsel;
6138 (partonIsUnknown[2] && partonIsUnknown[3])
6157 if (abs(isel)>=6)
break;
6166 ckm_takeout = pow(ckm_ir*ckm_js, 2);
6167 for (
unsigned int iswap=0; iswap<2; iswap++){
6168 msq_ddbar_ww_ijrs1234[iswap] /= ckm_takeout;
6169 msq_ddbar_ww_ijrs1243[iswap] /= ckm_takeout;
6177 (partonIsUnknown[0] && partonIsUnknown[1])
6191 (partonIsUnknown[2] && partonIsUnknown[3])
6207 const unsigned int nCombinations = 6;
6208 int ijrs_arr[nCombinations][4] ={
6216 for (
unsigned int ijrs=0; ijrs<nCombinations; ijrs++){
6217 isel = ijrs_arr[ijrs][0];
6218 jsel = ijrs_arr[ijrs][1];
6219 rsel = ijrs_arr[ijrs][2];
6220 ssel = ijrs_arr[ijrs][3];
6223 if (ckm_ir!=0. && ckm_js!=0.)
break;
6225 if (ckm_ir!=0. && ckm_js!=0.){
6230 ckm_takeout = pow(ckm_ir*ckm_js, 2);
6231 for (
unsigned int iswap=0; iswap<2; iswap++){
6232 msq_ud_wwonly_ijrs1234[iswap] /= ckm_takeout;
6233 msq_ud_wwonly_ijrs1243[iswap] /= ckm_takeout;
6240 (partonIsUnknown[0] && partonIsUnknown[1])
6254 (partonIsUnknown[2] && partonIsUnknown[3])
6270 const unsigned int nCombinations = 6;
6271 int ijrs_arr[nCombinations][4] ={
6279 for (
unsigned int ijrs=0; ijrs<nCombinations; ijrs++){
6280 isel = ijrs_arr[ijrs][0];
6281 jsel = ijrs_arr[ijrs][1];
6282 rsel = ijrs_arr[ijrs][2];
6283 ssel = ijrs_arr[ijrs][3];
6286 if (ckm_ir!=0. && ckm_js!=0.)
break;
6288 if (ckm_ir!=0. && ckm_js!=0.){
6293 ckm_takeout = pow(ckm_ir*ckm_js, 2);
6294 for (
unsigned int iswap=0; iswap<2; iswap++){
6295 msq_ubardbar_wwonly_ijrs1234[iswap] /= ckm_takeout;
6296 msq_ubardbar_wwonly_ijrs1243[iswap] /= ckm_takeout;
6302 MELAout <<
"TUtil::HJJMatEl: The pre-computed MEs:" << endl;
6303 for (
unsigned int iswap=0; iswap<2; iswap++){
6305 <<
"\tmsq_uu_zz_ijrs1234[" << iswap <<
"] = " << msq_uu_zz_ijrs1234[iswap] <<
'\n'
6306 <<
"\tmsq_uu_zz_ijrs1243[" << iswap <<
"] = " << msq_uu_zz_ijrs1243[iswap] <<
'\n'
6307 <<
"\tmsq_dd_zz_ijrs1234[" << iswap <<
"] = " << msq_dd_zz_ijrs1234[iswap] <<
'\n'
6308 <<
"\tmsq_dd_zz_ijrs1243[" << iswap <<
"] = " << msq_dd_zz_ijrs1243[iswap] <<
'\n'
6309 <<
"\tmsq_ubarubar_zz_ijrs1234[" << iswap <<
"] = " << msq_ubarubar_zz_ijrs1234[iswap] <<
'\n'
6310 <<
"\tmsq_ubarubar_zz_ijrs1243[" << iswap <<
"] = " << msq_ubarubar_zz_ijrs1243[iswap] <<
'\n'
6311 <<
"\tmsq_dbardbar_zz_ijrs1234[" << iswap <<
"] = " << msq_dbardbar_zz_ijrs1234[iswap] <<
'\n'
6312 <<
"\tmsq_dbardbar_zz_ijrs1243[" << iswap <<
"] = " << msq_dbardbar_zz_ijrs1243[iswap] <<
'\n'
6313 <<
"\tmsq_uu_zzid_ijrs1234[" << iswap <<
"] = " << msq_uu_zzid_ijrs1234[iswap] <<
'\n'
6314 <<
"\tmsq_uu_zzid_ijrs1243[" << iswap <<
"] = " << msq_uu_zzid_ijrs1243[iswap] <<
'\n'
6315 <<
"\tmsq_dd_zzid_ijrs1234[" << iswap <<
"] = " << msq_dd_zzid_ijrs1234[iswap] <<
'\n'
6316 <<
"\tmsq_dd_zzid_ijrs1243[" << iswap <<
"] = " << msq_dd_zzid_ijrs1243[iswap] <<
'\n'
6317 <<
"\tmsq_ubarubar_zzid_ijrs1234[" << iswap <<
"] = " << msq_ubarubar_zzid_ijrs1234[iswap] <<
'\n'
6318 <<
"\tmsq_ubarubar_zzid_ijrs1243[" << iswap <<
"] = " << msq_ubarubar_zzid_ijrs1243[iswap] <<
'\n'
6319 <<
"\tmsq_dbardbar_zzid_ijrs1234[" << iswap <<
"] = " << msq_dbardbar_zzid_ijrs1234[iswap] <<
'\n'
6320 <<
"\tmsq_dbardbar_zzid_ijrs1243[" << iswap <<
"] = " << msq_dbardbar_zzid_ijrs1243[iswap] <<
'\n'
6321 <<
"\tmsq_udbar_zz_ijrs1234[" << iswap <<
"] = " << msq_udbar_zz_ijrs1234[iswap] <<
'\n'
6322 <<
"\tmsq_udbar_zz_ijrs1243[" << iswap <<
"] = " << msq_udbar_zz_ijrs1243[iswap] <<
'\n'
6323 <<
"\tmsq_dubar_zz_ijrs1234[" << iswap <<
"] = " << msq_dubar_zz_ijrs1234[iswap] <<
'\n'
6324 <<
"\tmsq_dubar_zz_ijrs1243[" << iswap <<
"] = " << msq_dubar_zz_ijrs1243[iswap] <<
'\n'
6325 <<
"\tmsq_uubar_zz_ijrs1234[" << iswap <<
"] = " << msq_uubar_zz_ijrs1234[iswap] <<
'\n'
6326 <<
"\tmsq_uubar_zz_ijrs1243[" << iswap <<
"] = " << msq_uubar_zz_ijrs1243[iswap] <<
'\n'
6327 <<
"\tmsq_ddbar_zz_ijrs1234[" << iswap <<
"] = " << msq_ddbar_zz_ijrs1234[iswap] <<
'\n'
6328 <<
"\tmsq_ddbar_zz_ijrs1243[" << iswap <<
"] = " << msq_ddbar_zz_ijrs1243[iswap] <<
'\n'
6329 <<
"\tmsq_uubar_ww_ijrs1234[" << iswap <<
"] = " << msq_uubar_ww_ijrs1234[iswap] <<
'\n'
6330 <<
"\tmsq_uubar_ww_ijrs1243[" << iswap <<
"] = " << msq_uubar_ww_ijrs1243[iswap] <<
'\n'
6331 <<
"\tmsq_ddbar_ww_ijrs1234[" << iswap <<
"] = " << msq_ddbar_ww_ijrs1234[iswap] <<
'\n'
6332 <<
"\tmsq_ddbar_ww_ijrs1243[" << iswap <<
"] = " << msq_ddbar_ww_ijrs1243[iswap] <<
'\n'
6333 <<
"\tmsq_ud_wwonly_ijrs1234[" << iswap <<
"] = " << msq_ud_wwonly_ijrs1234[iswap] <<
'\n'
6334 <<
"\tmsq_ud_wwonly_ijrs1243[" << iswap <<
"] = " << msq_ud_wwonly_ijrs1243[iswap] <<
'\n'
6335 <<
"\tmsq_ubardbar_wwonly_ijrs1234[" << iswap <<
"] = " << msq_ubardbar_wwonly_ijrs1234[iswap] <<
'\n'
6336 <<
"\tmsq_ubardbar_wwonly_ijrs1243[" << iswap <<
"] = " << msq_ubardbar_wwonly_ijrs1243[iswap]
6342 int nijchannels = ijsel.size();
6344 for (
int ic=0; ic<nijchannels; ic++){
6346 isel = ijsel[ic][0];
6347 jsel = ijsel[ic][1];
6348 int code = ijsel[ic][2];
6349 bool ijselIsUpType[2];
6350 bool ijselIsDownType[2];
6351 bool ijselIsParticle[2];
6352 bool ijselIsAntiparticle[2];
6366 ijselIsParticle[0] = (isel>0); ijselIsAntiparticle[0] = (isel<0);
6367 ijselIsParticle[1] = (jsel>0); ijselIsAntiparticle[1] = (jsel<0);
6370 (partonIsUnknown[0] || MYIDUP_tmp[0]==isel)
6372 (partonIsUnknown[1] || MYIDUP_tmp[1]==jsel)
6377 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6379 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6382 if (ijselIsUpType[0] && ijselIsUpType[1]){
6383 if (ijselIsParticle[0] && ijselIsParticle[1]){
6384 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1234[0];
6385 else msq_tmp = msq_uu_zzid_ijrs1234[0];
6387 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6388 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1234[0];
6389 else msq_tmp = msq_ubarubar_zzid_ijrs1234[0];
6391 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_zz_ijrs1234[0];
6393 else if (ijselIsDownType[0] && ijselIsDownType[1]){
6394 if (ijselIsParticle[0] && ijselIsParticle[1]){
6395 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1234[0];
6396 else msq_tmp = msq_dd_zzid_ijrs1234[0];
6398 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6399 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1234[0];
6400 else msq_tmp = msq_dbardbar_zzid_ijrs1234[0];
6402 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_zz_ijrs1234[0];
6404 else if (ijselIsUpType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_udbar_zz_ijrs1234[0];
6405 else if (ijselIsDownType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_dubar_zz_ijrs1234[0];
6406 MatElsq[jsel+5][isel+5] += msq_tmp;
6407 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6412 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6414 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6417 if (ijselIsUpType[0] && ijselIsUpType[1]){
6418 if (ijselIsParticle[0] && ijselIsParticle[1]){
6419 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1243[0];
6420 else msq_tmp = msq_uu_zzid_ijrs1243[0];
6422 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6423 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1243[0];
6424 else msq_tmp = msq_ubarubar_zzid_ijrs1243[0];
6426 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_zz_ijrs1243[0];
6428 else if (ijselIsDownType[0] && ijselIsDownType[1]){
6429 if (ijselIsParticle[0] && ijselIsParticle[1]){
6430 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1243[0];
6431 else msq_tmp = msq_dd_zzid_ijrs1243[0];
6433 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6434 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1243[0];
6435 else msq_tmp = msq_dbardbar_zzid_ijrs1243[0];
6437 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_zz_ijrs1243[0];
6439 else if (ijselIsUpType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_udbar_zz_ijrs1243[0];
6440 else if (ijselIsDownType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_dubar_zz_ijrs1243[0];
6441 MatElsq[jsel+5][isel+5] += msq_tmp;
6442 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6446 vector<int> possible_rsel;
6447 vector<int> possible_ssel;
6448 vector<double> possible_Vsqir;
6449 vector<double> possible_Vsqjs;
6450 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6451 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6452 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6453 rsel=possible_rsel.at(ix);
6455 possible_Vsqir.push_back(ckmval);
6457 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6458 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6459 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6460 ssel=possible_ssel.at(iy);
6462 possible_Vsqjs.push_back(ckmval);
6466 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6467 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6468 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6469 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6470 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6472 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6474 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6477 if (ijselIsUpType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_ww_ijrs1234[0];
6478 else if (ijselIsDownType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_ww_ijrs1234[0];
6480 MatElsq[jsel+5][isel+5] += msq_tmp;
6481 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6486 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6488 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6491 if (ijselIsUpType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_ww_ijrs1243[0];
6492 else if (ijselIsDownType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_ww_ijrs1243[0];
6494 MatElsq[jsel+5][isel+5] += msq_tmp;
6495 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6501 vector<int> possible_rsel;
6502 vector<int> possible_ssel;
6503 vector<double> possible_Vsqir;
6504 vector<double> possible_Vsqjs;
6505 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6506 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6507 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6508 rsel=possible_rsel.at(ix);
6510 possible_Vsqir.push_back(ckmval);
6512 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6513 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6514 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6515 ssel=possible_ssel.at(iy);
6517 possible_Vsqjs.push_back(ckmval);
6521 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6522 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6523 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6524 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6525 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6527 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6529 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6531 if (rsel!=jsel && ssel!=isel){
6533 if (ijselIsUpType[0] && ijselIsDownType[1]){
6534 if (ijselIsParticle[0] && ijselIsParticle[1]) msq_tmp = msq_ud_wwonly_ijrs1234[0];
6535 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ubardbar_wwonly_ijrs1234[0];
6538 MatElsq[jsel+5][isel+5] += msq_tmp;
6539 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6543 MatElsq[jsel+5][isel+5] += msq_tmp;
6544 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6550 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6552 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6554 if (rsel!=jsel && ssel!=isel){
6556 if (ijselIsUpType[0] && ijselIsDownType[1]){
6557 if (ijselIsParticle[0] && ijselIsParticle[1]) msq_tmp = msq_ud_wwonly_ijrs1243[0];
6558 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ubardbar_wwonly_ijrs1243[0];
6561 MatElsq[jsel+5][isel+5] += msq_tmp;
6562 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6566 MatElsq[jsel+5][isel+5] += msq_tmp;
6567 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6575 if (isel==jsel)
continue;
6576 isel = ijsel[ic][1];
6577 jsel = ijsel[ic][0];
6589 ijselIsParticle[0] = (isel>0); ijselIsAntiparticle[0] = (isel<0);
6590 ijselIsParticle[1] = (jsel>0); ijselIsAntiparticle[1] = (jsel<0);
6593 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
6595 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
6601 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6603 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6606 if (ijselIsUpType[1] && ijselIsUpType[0]){
6607 if (ijselIsParticle[1] && ijselIsParticle[0]){
6608 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1234[1];
6609 else msq_tmp = msq_uu_zzid_ijrs1234[1];
6611 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6612 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1234[1];
6613 else msq_tmp = msq_ubarubar_zzid_ijrs1234[1];
6615 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_zz_ijrs1234[1];
6617 else if (ijselIsDownType[1] && ijselIsDownType[0]){
6618 if (ijselIsParticle[1] && ijselIsParticle[0]){
6619 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1234[1];
6620 else msq_tmp = msq_dd_zzid_ijrs1234[1];
6622 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6623 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1234[1];
6624 else msq_tmp = msq_dbardbar_zzid_ijrs1234[1];
6626 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_zz_ijrs1234[1];
6628 else if (ijselIsUpType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_udbar_zz_ijrs1234[1];
6629 else if (ijselIsDownType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_dubar_zz_ijrs1234[1];
6630 MatElsq[jsel+5][isel+5] += msq_tmp;
6631 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6636 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6638 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6641 if (ijselIsUpType[1] && ijselIsUpType[0]){
6642 if (ijselIsParticle[1] && ijselIsParticle[0]){
6643 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1243[1];
6644 else msq_tmp = msq_uu_zzid_ijrs1243[1];
6646 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6647 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1243[1];
6648 else msq_tmp = msq_ubarubar_zzid_ijrs1243[1];
6650 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_zz_ijrs1243[1];
6652 else if (ijselIsDownType[1] && ijselIsDownType[0]){
6653 if (ijselIsParticle[1] && ijselIsParticle[0]){
6654 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1243[1];
6655 else msq_tmp = msq_dd_zzid_ijrs1243[1];
6657 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6658 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1243[1];
6659 else msq_tmp = msq_dbardbar_zzid_ijrs1243[1];
6661 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_zz_ijrs1243[1];
6663 else if (ijselIsUpType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_udbar_zz_ijrs1243[1];
6664 else if (ijselIsDownType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_dubar_zz_ijrs1243[1];
6665 MatElsq[jsel+5][isel+5] += msq_tmp;
6666 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6670 vector<int> possible_rsel;
6671 vector<int> possible_ssel;
6672 vector<double> possible_Vsqir;
6673 vector<double> possible_Vsqjs;
6674 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6675 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6676 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6677 rsel=possible_rsel.at(ix);
6679 possible_Vsqir.push_back(ckmval);
6681 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6682 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6683 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6684 ssel=possible_ssel.at(iy);
6686 possible_Vsqjs.push_back(ckmval);
6690 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6691 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6692 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6693 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6694 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6696 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6698 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6701 if (ijselIsUpType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_ww_ijrs1234[1];
6702 else if (ijselIsDownType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_ww_ijrs1234[1];
6704 MatElsq[jsel+5][isel+5] += msq_tmp;
6705 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6710 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6712 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6715 if (ijselIsUpType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_ww_ijrs1243[1];
6716 else if (ijselIsDownType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_ww_ijrs1243[1];
6718 MatElsq[jsel+5][isel+5] += msq_tmp;
6719 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6725 vector<int> possible_rsel;
6726 vector<int> possible_ssel;
6727 vector<double> possible_Vsqir;
6728 vector<double> possible_Vsqjs;
6729 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6730 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6731 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6732 rsel=possible_rsel.at(ix);
6734 possible_Vsqir.push_back(ckmval);
6736 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6737 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6738 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6739 ssel=possible_ssel.at(iy);
6741 possible_Vsqjs.push_back(ckmval);
6745 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6746 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6747 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6748 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6749 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6751 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6753 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6755 if (rsel!=jsel && ssel!=isel){
6757 if (ijselIsUpType[1] && ijselIsDownType[0]){
6758 if (ijselIsParticle[1] && ijselIsParticle[0]) msq_tmp = msq_ud_wwonly_ijrs1234[1];
6759 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ubardbar_wwonly_ijrs1234[1];
6762 MatElsq[jsel+5][isel+5] += msq_tmp;
6763 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6767 MatElsq[jsel+5][isel+5] += msq_tmp;
6768 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6774 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6776 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6778 if (rsel!=jsel && ssel!=isel){
6780 if (ijselIsUpType[1] && ijselIsDownType[0]){
6781 if (ijselIsParticle[1] && ijselIsParticle[0]) msq_tmp = msq_ud_wwonly_ijrs1243[1];
6782 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ubardbar_wwonly_ijrs1243[1];
6785 MatElsq[jsel+5][isel+5] += msq_tmp;
6786 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6790 MatElsq[jsel+5][isel+5] += msq_tmp;
6791 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6802 int GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets)*2;
6803 double constant = pow(GeV, -GeVexponent_MEsq);
6804 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
6810 for (
int ii = 0; ii <
nmsq; ii++){
for (
int jj = 0; jj <
nmsq; jj++)
MELAout << MatElsq[jj][ii] <<
'\t';
MELAout << endl; }
6812 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
6829 <<
"TUtil::HJJMatEl: Reset AlphaS:\n"
6830 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
6832 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
6836 <<
"TUtil::HJJMatEl: Reset AlphaS result:\n"
6837 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
6846 const double& EBEAM,
6847 bool includeHiggsDecay,
6850 const double GeV=1./100.;
6851 double sum_msqjk = 0;
6859 double MatElsq[
nmsq][
nmsq]={ { 0 } };
6861 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::VHiggsMatEl: Non-JHUGen MEs are not supported" << endl;
return sum_msqjk; }
6864 int nRequested_AssociatedJets=0;
6865 int nRequested_AssociatedLeptons=0;
6866 int nRequested_AssociatedPhotons=0;
6867 int AssociationVCompatibility=0;
6871 nRequested_AssociatedJets=2;
6875 nRequested_AssociatedLeptons=2;
6879 nRequested_AssociatedPhotons=1;
6883 else if (production ==
TVar::GammaH) AssociationVCompatibility=22;
6895 if ((mela_event.
pAssociated.size()<(
unsigned int)(nRequested_AssociatedJets+nRequested_AssociatedLeptons) && production!=
TVar::GammaH) || (mela_event.
pAssociated.size()<(
unsigned int)nRequested_AssociatedPhotons && production ==
TVar::GammaH)){
6897 MELAerr <<
"TUtil::VHiggsMatEl: Number of associated particles (" << mela_event.
pAssociated.size() <<
") is less than ";
6898 if (production!=
TVar::GammaH)
MELAerr << (nRequested_AssociatedJets+nRequested_AssociatedLeptons);
6899 else MELAerr << nRequested_AssociatedPhotons;
6905 int MYIDUP_prod[4]={ 0 };
6906 int MYIDUP_dec[2]={ -9000, -9000 };
6907 double p4[9][4] ={ { 0 } };
6908 double helicities[9] ={ 0 };
6909 int vh_ids[9] ={ 0 };
6910 TLorentzVector MomStore[
mxpart];
6911 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
6919 for (
int ipar=0; ipar<2; ipar++){
6920 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
6921 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
6923 else MYIDUP_prod[ipar] = 0;
6924 if (momTmp->T()>0.){
6925 p4[ipar][0] = momTmp->T()*GeV;
6926 p4[ipar][1] = momTmp->X()*GeV;
6927 p4[ipar][2] = momTmp->Y()*GeV;
6928 p4[ipar][3] = momTmp->Z()*GeV;
6929 MomStore[ipar] = (*momTmp);
6932 p4[ipar][0] = -momTmp->T()*GeV;
6933 p4[ipar][1] = -momTmp->X()*GeV;
6934 p4[ipar][2] = -momTmp->Y()*GeV;
6935 p4[ipar][3] = -momTmp->Z()*GeV;
6936 MomStore[ipar] = -(*momTmp);
6937 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
6941 for (
int ipar=0; ipar<(production!=
TVar::GammaH ? 2 : 1); ipar++){
6942 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(ipar).second);
6943 int* idtmp = &(mela_event.
pAssociated.at(ipar).first);
6945 else MYIDUP_prod[ipar+2] = 0;
6946 p4[ipar+5][0] = momTmp->T()*GeV;
6947 p4[ipar+5][1] = momTmp->X()*GeV;
6948 p4[ipar+5][2] = momTmp->Y()*GeV;
6949 p4[ipar+5][3] = momTmp->Z()*GeV;
6950 MomStore[ipar+6] = (*momTmp);
6956 if (production ==
TVar::GammaH && !
PDGHelpers::isAPhoton(MYIDUP_prod[2])){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::VHiggsMatEl: GammaH associated photon (id=" << MYIDUP_prod[2] <<
") is not a photon! Please fix its id." << endl;
return sum_msqjk; }
6960 for (
int iv=0; iv<2; iv++){
6963 else MYIDUP_dec[iv] = 0;
6966 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
6967 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
6969 p4[ipar+7][0] = momTmp->T()*GeV;
6970 p4[ipar+7][1] = momTmp->X()*GeV;
6971 p4[ipar+7][2] = momTmp->Y()*GeV;
6972 p4[ipar+7][3] = momTmp->Z()*GeV;
6973 MomStore[5] = (*momTmp);
6976 p4[ipar+7][0] = momTmp->T()*GeV;
6977 p4[ipar+7][1] = momTmp->X()*GeV;
6978 p4[ipar+7][2] = momTmp->Y()*GeV;
6979 p4[ipar+7][3] = momTmp->Z()*GeV;
6980 MomStore[2*ipar+2] = (*momTmp);
6985 p4[7][0] += momTmp->T()*GeV;
6986 p4[7][1] += momTmp->X()*GeV;
6987 p4[7][2] += momTmp->Y()*GeV;
6988 p4[7][3] += momTmp->Z()*GeV;
6991 p4[8][0] = momTmp->T()*GeV;
6992 p4[8][1] = momTmp->X()*GeV;
6993 p4[8][2] = momTmp->Y()*GeV;
6994 p4[8][3] = momTmp->Z()*GeV;
6996 MomStore[ipar+2] = (*momTmp);
7000 p4[7][0] += momTmp->T()*GeV;
7001 p4[7][1] += momTmp->X()*GeV;
7002 p4[7][2] += momTmp->Y()*GeV;
7003 p4[7][3] += momTmp->Z()*GeV;
7006 p4[8][0] += momTmp->T()*GeV;
7007 p4[8][1] += momTmp->X()*GeV;
7008 p4[8][2] += momTmp->Y()*GeV;
7009 p4[8][3] += momTmp->Z()*GeV;
7011 MomStore[ipar+2] = (*momTmp);
7014 p4[7][0] += momTmp->T()*GeV;
7015 p4[7][1] += momTmp->X()*GeV;
7016 p4[7][2] += momTmp->Y()*GeV;
7017 p4[7][3] += momTmp->Z()*GeV;
7018 MomStore[5] = MomStore[5] + (*momTmp);
7021 for (
int ix=0; ix<4; ix++){
7022 p4[3][ix] = p4[5][ix] + p4[6][ix];
7023 p4[4][ix] = p4[7][ix] + p4[8][ix];
7024 p4[2][ix] = p4[3][ix] + p4[4][ix];
7031 else vh_ids[3] = 22;
7034 vh_ids[7] = 5; helicities[7] = 1;
7035 vh_ids[8] = -5; helicities[8] = 1;
7037 if (includeHiggsDecay && MYIDUP_dec[0]!=-9000 && MYIDUP_dec[1]!=-9000 && MYIDUP_dec[0]==-MYIDUP_dec[1]){
7042 else if (verbosity>=
TVar::INFO && includeHiggsDecay)
MELAerr <<
"TUtil::VHiggsMatEl: includeHiggsDecay=true is not supported for the present decay mode." << endl;
7045 for (
int i=0; i<9; i++)
MELAout <<
"p4(" << vh_ids[i] <<
") = " << p4[i][0] <<
", " << p4[i][1] <<
", " << p4[i][2] <<
", " << p4[i][3] <<
"\n";
7048 double defaultRenScale =
scale_.scale;
7049 double defaultFacScale =
facscale_.facscale;
7051 int defaultNflav =
nflav_.nflav;
7052 string defaultPdflabel =
pdlabel_.pdlabel;
7056 double alphasVal, alphasmzVal;
7066 <<
"TUtil::VHiggsMatEl: Set AlphaS:\n"
7067 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
7069 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7073 bool partonIsKnown[4];
7074 for (
unsigned int ip=0; ip<4; ip++) partonIsKnown[ip] = (MYIDUP_prod[ip]!=0);
7075 if ((production ==
TVar::Lep_WH || production ==
TVar::Lep_ZH || production ==
TVar::GammaH) && !(partonIsKnown[2] && partonIsKnown[3])){
if (verbosity>=
TVar::INFO)
MELAerr <<
"TUtil::VHiggsMatEl: Final state particles in leptonic/photonic VH have to have a definite id!" << endl;
return sum_msqjk; }
7077 const double allowed_helicities[2] ={ -1, 1 };
7079 vector<pair<int, int>> Hffparticles;
7084 else Hffparticles.push_back(pair<int, int>(-MYIDUP_dec[1], MYIDUP_dec[1]));
7087 double Hffscalesum=0;
7088 for (
int hquark=1; hquark<=5; hquark++){
7090 Hffscalesum += Hffmass;
7092 Hffparticles.push_back(pair<int, int>(hquark, -hquark));
7093 Hffparticles.push_back(pair<int, int>(-hquark, hquark));
7094 Hffscale /= Hffmass;
7097 Hffscale *= Hffscalesum;
7101 MELAout <<
"TUtil::VHiggsMatEl: Outgoing H-> f fbar particles to compute for the ME template:" << endl;
7102 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++)
MELAout <<
"\t - (id8, id9) = (" << Hffparticles.at(ihf).first <<
", " << Hffparticles.at(ihf).second <<
")" << endl;
7103 MELAout <<
"TUtil::VHiggsMatEl: ME scale for the H-> f fbar particles: " << Hffscale << endl;
7108 vector<pair<int, int>> incomingPartons;
7109 if (partonIsKnown[0] && partonIsKnown[1]) incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], MYIDUP_prod[1]));
7111 else if (!partonIsKnown[0] && !partonIsKnown[1]){
7113 incomingPartons.push_back(pair<int, int>(1, -2));
7114 incomingPartons.push_back(pair<int, int>(-2, 1));
7115 incomingPartons.push_back(pair<int, int>(2, -1));
7116 incomingPartons.push_back(pair<int, int>(-1, 2));
7118 else if (!partonIsKnown[1] && partonIsKnown[0]){
7121 for (
int iqf=1; iqf<=5; iqf++){
7122 if (iqf%2==0)
continue;
7123 int jqf = -TMath::Sign(iqf, MYIDUP_prod[0]);
7125 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], jqf));
break; }
7129 for (
int iqf=2; iqf<=4; iqf++){
7130 if (iqf%2==1)
continue;
7131 int jqf = -TMath::Sign(iqf, MYIDUP_prod[0]);
7133 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], jqf));
break; }
7140 for (
int iqf=1; iqf<=5; iqf++){
7141 if (iqf%2==0)
continue;
7142 int jqf = -TMath::Sign(iqf, MYIDUP_prod[1]);
7144 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[1]));
break; }
7148 for (
int iqf=2; iqf<=4; iqf++){
7149 if (iqf%2==1)
continue;
7150 int jqf = -TMath::Sign(iqf, MYIDUP_prod[1]);
7152 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[1]));
break; }
7157 MELAout <<
"TUtil::VHiggsMatEl: Incoming partons to compute for the ME template:" << endl;
7158 for (
auto& inpart:incomingPartons)
MELAout <<
"\t - (id1, id2) = (" << inpart.first <<
", " << inpart.second <<
")" << endl;
7162 vector<pair<int, int>> outgoingPartons;
7163 if (partonIsKnown[2] && partonIsKnown[3]) outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], MYIDUP_prod[3]));
7165 else if (!partonIsKnown[2] && !partonIsKnown[3]){
7167 outgoingPartons.push_back(pair<int, int>(1, -2));
7168 outgoingPartons.push_back(pair<int, int>(-2, 1));
7169 outgoingPartons.push_back(pair<int, int>(2, -1));
7170 outgoingPartons.push_back(pair<int, int>(-1, 2));
7172 else if (!partonIsKnown[3] && partonIsKnown[2]){
7175 for (
int iqf=1; iqf<=5; iqf++){
7176 if (iqf%2==0)
continue;
7177 int jqf = -TMath::Sign(iqf, MYIDUP_prod[2]);
7179 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], jqf));
break; }
7183 for (
int iqf=2; iqf<=4; iqf++){
7184 if (iqf%2==1)
continue;
7185 int jqf = -TMath::Sign(iqf, MYIDUP_prod[2]);
7187 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], jqf));
break; }
7194 for (
int iqf=1; iqf<=5; iqf++){
7195 if (iqf%2==0)
continue;
7196 int jqf = -TMath::Sign(iqf, MYIDUP_prod[3]);
7198 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[3]));
break; }
7202 for (
int iqf=2; iqf<=4; iqf++){
7203 if (iqf%2==1)
continue;
7204 int jqf = -TMath::Sign(iqf, MYIDUP_prod[3]);
7206 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[3]));
break; }
7211 MELAout <<
"TUtil::VHiggsMatEl: Outgoing particles to compute for the ME template:" << endl;
7212 for (
auto& outpart:outgoingPartons)
MELAout <<
"\t - (id6, id7) = (" << outpart.first <<
", " << outpart.second <<
")" << endl;
7215 for (
auto& inpart:incomingPartons){
7216 vh_ids[0] = inpart.first;
7217 vh_ids[1] = inpart.second;
7220 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tIncoming " << vh_ids[0] <<
"," << vh_ids[1] <<
" -> " << vh_ids[2] << endl;
7223 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tNeed to divide the ME by |VCKM_incoming|**2 = " << Vckmsq_in << endl;
7225 for (
auto& outpart:outgoingPartons){
7226 vh_ids[5] = outpart.first;
7227 vh_ids[6] = outpart.second;
7229 if (vh_ids[2]!=vh_ids[3])
continue;
7230 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tOutgoing " << vh_ids[3] <<
" -> " << vh_ids[5] <<
"," << vh_ids[6] << endl;
7234 for (
int h01 = 0; h01 < 2; h01++){
7235 helicities[0] = allowed_helicities[h01];
7236 helicities[1] = -helicities[0];
7237 for (
int h56 = 0; h56 < 2; h56++){
7238 helicities[5] = allowed_helicities[h56];
7239 helicities[6] = -helicities[5];
7244 for (
int h78=0; h78<2; h78++){
7245 helicities[7]=allowed_helicities[h78];
7246 helicities[8]=allowed_helicities[h78];
7247 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++){
7248 vh_ids[7]=Hffparticles.at(ihf).first;
7249 vh_ids[8]=Hffparticles.at(ihf).second;
7250 double msq_inst_LR=0;
7252 msq_inst += msq_inst_LR;
7255 msq_inst *= Hffscale;
7262 double scalesum_out=0;
7263 if (!(partonIsKnown[2] && partonIsKnown[3])){
7265 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tDividing ME by |VCKM_outgoing|**2 = " << Vckmsq_out << endl;
7267 for (
int outgoing1=1; outgoing1<=
nf; outgoing1++){
7268 if (partonIsKnown[2] && outgoing1!=abs(vh_ids[5]))
continue;
7269 if (outgoing1%2!=abs(vh_ids[5])%2 || outgoing1==6)
continue;
7270 int iout = outgoing1 * TMath::Sign(1, vh_ids[5]);
7271 for (
int outgoing2=1; outgoing2<=
nf; outgoing2++){
7272 if (partonIsKnown[3] && outgoing2!=abs(vh_ids[6]))
continue;
7273 if (outgoing2%2!=abs(vh_ids[6])%2 || outgoing2==6)
continue;
7274 int jout = outgoing2 * TMath::Sign(1, vh_ids[6]);
7279 else scalesum_out = 1;
7280 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for outgoing particles: " << scalesum_out << endl;
7283 if (!partonIsKnown[0] || !partonIsKnown[1]){
7284 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tDividing ME by |VCKM_incoming|**2 = " << Vckmsq_in << endl;
7289 for (
int incoming1=1; incoming1<=
nf; incoming1++){
7290 if (partonIsKnown[0] && incoming1!=abs(vh_ids[0]))
continue;
7291 if (incoming1%2!=abs(vh_ids[0])%2 || incoming1==6)
continue;
7292 int iin = incoming1 * TMath::Sign(1, vh_ids[0]);
7293 for (
int incoming2=1; incoming2<=
nf; incoming2++){
7294 if (partonIsKnown[1] && incoming2!=abs(vh_ids[1]))
continue;
7295 if (incoming2%2!=abs(vh_ids[1])%2 || incoming2==6)
continue;
7296 int jin = incoming2 * TMath::Sign(1, vh_ids[1]);
7297 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\t\tiin,jin = " << iin <<
"," << jin << endl;
7299 if (!partonIsKnown[0] || !partonIsKnown[1]){
7301 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for incoming particles: " << scale_in << endl;
7303 MatElsq[jin+5][iin+5] +=
msq * 0.25 * scale_in *scalesum_out;
7312 vector<pair<int, int>> incomingPartons;
7313 if (partonIsKnown[0] && partonIsKnown[1]) incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], MYIDUP_prod[1]));
7314 else if (!partonIsKnown[0] && !partonIsKnown[1]){
7316 incomingPartons.push_back(pair<int, int>(1, -1));
7317 incomingPartons.push_back(pair<int, int>(-1, 1));
7318 incomingPartons.push_back(pair<int, int>(2, -2));
7319 incomingPartons.push_back(pair<int, int>(-2, 2));
7321 else if (!partonIsKnown[1] && partonIsKnown[0])
7322 incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], -MYIDUP_prod[0]));
7324 incomingPartons.push_back(pair<int, int>(-MYIDUP_prod[1], MYIDUP_prod[1]));
7326 MELAout <<
"TUtil::VHiggsMatEl: Incoming partons to compute for the ME template:" << endl;
7327 for (
auto& inpart:incomingPartons)
MELAout <<
"\t - (id1, id2) = (" << inpart.first <<
", " << inpart.second <<
")" << endl;
7331 vector<pair<int, int>> outgoingPartons;
7332 if ((partonIsKnown[2] && partonIsKnown[3]) || production ==
TVar::GammaH) outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], MYIDUP_prod[3]));
7333 else if (!partonIsKnown[2] && !partonIsKnown[3]){
7335 outgoingPartons.push_back(pair<int, int>(1, -1));
7336 outgoingPartons.push_back(pair<int, int>(-1, 1));
7337 outgoingPartons.push_back(pair<int, int>(2, -2));
7338 outgoingPartons.push_back(pair<int, int>(-2, 2));
7340 else if (!partonIsKnown[3] && partonIsKnown[2])
7341 outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], -MYIDUP_prod[2]));
7343 outgoingPartons.push_back(pair<int, int>(-MYIDUP_prod[3], MYIDUP_prod[3]));
7345 MELAout <<
"TUtil::VHiggsMatEl: Outgoing particles to compute for the ME template:" << endl;
7346 for (
unsigned int op=0; op<outgoingPartons.size(); op++)
MELAout <<
"\t - (id6, id7) = (" << outgoingPartons.at(op).first <<
", " << outgoingPartons.at(op).second <<
")" << endl;
7349 for (
auto& inpart:incomingPartons){
7350 vh_ids[0] = inpart.first;
7351 vh_ids[1] = inpart.second;
7354 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tIncoming " << vh_ids[0] <<
"," << vh_ids[1] <<
" -> " << vh_ids[2] << endl;
7358 for (
auto& outpart:outgoingPartons){
7359 vh_ids[5] = outpart.first;
7360 vh_ids[6] = outpart.second;
7364 if (vh_ids[2]!=vh_ids[3])
continue;
7366 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tOutgoing " << vh_ids[3] <<
" -> " << vh_ids[5] <<
"," << vh_ids[6] << endl;
7370 for (
int h01 = 0; h01 < 2; h01++){
7371 helicities[0] = allowed_helicities[h01];
7372 helicities[1] = -helicities[0];
7373 for (
int h56 = 0; h56 < 2; h56++){
7374 helicities[5] = allowed_helicities[h56];
7375 helicities[6] = -helicities[5];
7380 for (
int h78=0; h78<2; h78++){
7381 helicities[7]=allowed_helicities[h78];
7382 helicities[8]=allowed_helicities[h78];
7383 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++){
7384 vh_ids[7]=Hffparticles.at(ihf).first;
7385 vh_ids[8]=Hffparticles.at(ihf).second;
7386 double msq_inst_LR=0;
7388 msq_inst += msq_inst_LR;
7391 msq_inst *= Hffscale;
7399 if (!partonIsKnown[2] && !partonIsKnown[3] && production!=
TVar::GammaH){
7403 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for outgoing particles: " << scale_out << endl;
7406 for (
int incoming1=1; incoming1<=
nf; incoming1++){
7407 if (partonIsKnown[0] && incoming1!=abs(vh_ids[0]))
continue;
7408 if (incoming1%2!=abs(vh_ids[0])%2 || incoming1==6)
continue;
7409 int iin = incoming1 * TMath::Sign(1, vh_ids[0]);
7411 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\t\tiin,jin = " << iin <<
"," << jin << endl;
7412 MatElsq[jin+5][iin+5] +=
msq * 0.25 * scale_out;
7419 int GeVexponent_MEsq;
7420 if (HDKon==0) GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets+nRequested_AssociatedLeptons+nRequested_AssociatedPhotons)*2-4;
7421 else GeVexponent_MEsq = 4-(2+nRequested_AssociatedJets+nRequested_AssociatedLeptons+nRequested_AssociatedPhotons)*2;
7422 double constant = pow(GeV, -GeVexponent_MEsq);
7423 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7426 for (
int ii = 0; ii <
nmsq; ii++){
for (
int jj = 0; jj <
nmsq; jj++)
MELAout << MatElsq[jj][ii] <<
'\t';
MELAout << endl; }
7428 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7438 <<
"TUtil::VHiggsMatEl: Reset AlphaS:\n"
7439 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7441 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7445 <<
"TUtil::VHiggsMatEl: Reset AlphaS result:\n"
7446 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7455 const double& EBEAM,
7456 int topDecay,
int topProcess,
7459 const double GeV=1./100.;
7460 double sum_msqjk = 0;
7461 double MatElsq[
nmsq][
nmsq]={ { 0 } };
7463 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::TTHiggsMatEl: Non-JHUGen MEs are not supported." << endl;
return sum_msqjk; }
7464 if (production!=
TVar::ttH){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::TTHiggsMatEl: Only ttH is supported." << endl;
return sum_msqjk; }
7467 int nRequested_Tops=1;
7468 int nRequested_Antitops=1;
7485 <<
"TUtil::TTHiggsMatEl: Number of stable tops (" << mela_event.
pStableTops.size() <<
")"
7486 <<
" and number of stable antitops (" << mela_event.
pStableAntitops.size() <<
")"
7487 <<
" in ttH process are not 1!" << endl;
7492 <<
"TUtil::TTHiggsMatEl: Number of set of top daughters (" << mela_event.
pTopDaughters.size() <<
")"
7493 <<
" and number of set of antitop daughters (" << mela_event.
pAntitopDaughters.size() <<
")"
7494 <<
" in ttH process are not 1!" << endl;
7499 <<
"TUtil::TTHiggsMatEl: Number of top daughters (" << mela_event.
pTopDaughters.at(0).size() <<
")"
7500 <<
" and number of antitop daughters (" << mela_event.
pAntitopDaughters.at(0).size() <<
")"
7501 <<
" in ttH process are not 3!" << endl;
7507 bool isUnknown[2]={
true,
true };
7511 for (
unsigned int itd=0; itd<mela_event.
pTopDaughters.at(0).size(); itd++) topDaughters.push_back(mela_event.
pTopDaughters.at(0).at(itd));
7515 for (
unsigned int itop=0; itop<mela_event.
pStableTops.size(); itop++) topDaughters.push_back(mela_event.
pStableTops.at(itop));
7519 for (
unsigned int itd=0; itd<topDaughters.size(); itd++){
if (!
PDGHelpers::isAnUnknownJet(topDaughters.at(itd).first)){ isUnknown[0]=
false;
break; } }
7520 for (
unsigned int itd=0; itd<antitopDaughters.size(); itd++){
if (!
PDGHelpers::isAnUnknownJet(antitopDaughters.at(itd).first)){ isUnknown[1]=
false;
break; } }
7527 double p4[13][4]={ { 0 } };
7528 double MYIDUP_prod[2]={ 0 };
7529 TLorentzVector MomStore[
mxpart];
7530 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
7531 for (
int ipar=0; ipar<2; ipar++){
7532 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
7533 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
7535 else MYIDUP_prod[ipar] = 0;
7536 if (momTmp->T()>0.){
7537 p4[ipar][0] = -momTmp->T()*GeV;
7538 p4[ipar][1] = -momTmp->X()*GeV;
7539 p4[ipar][2] = -momTmp->Y()*GeV;
7540 p4[ipar][3] = -momTmp->Z()*GeV;
7541 MomStore[ipar] = (*momTmp);
7544 p4[ipar][0] = momTmp->T()*GeV;
7545 p4[ipar][1] = momTmp->X()*GeV;
7546 p4[ipar][2] = momTmp->Y()*GeV;
7547 p4[ipar][3] = momTmp->Z()*GeV;
7548 MomStore[ipar] = -(*momTmp);
7549 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
7555 const unsigned int t_pos=4;
7556 const unsigned int b_pos=9;
7557 const unsigned int Wp_pos=10;
7558 const unsigned int Wpf_pos=12;
7559 const unsigned int Wpfb_pos=11;
7560 for (
unsigned int ipar=0; ipar<topDaughters.size(); ipar++){
7561 TLorentzVector* momTmp = &(topDaughters.at(ipar).second);
7562 if (topDaughters.size()==1){
7563 p4[t_pos][0] = momTmp->T()*GeV;
7564 p4[t_pos][1] = momTmp->X()*GeV;
7565 p4[t_pos][2] = momTmp->Y()*GeV;
7566 p4[t_pos][3] = momTmp->Z()*GeV;
7570 p4[b_pos][0] = momTmp->T()*GeV;
7571 p4[b_pos][1] = momTmp->X()*GeV;
7572 p4[b_pos][2] = momTmp->Y()*GeV;
7573 p4[b_pos][3] = momTmp->Z()*GeV;
7576 p4[Wpfb_pos][0] = momTmp->T()*GeV;
7577 p4[Wpfb_pos][1] = momTmp->X()*GeV;
7578 p4[Wpfb_pos][2] = momTmp->Y()*GeV;
7579 p4[Wpfb_pos][3] = momTmp->Z()*GeV;
7580 p4[Wp_pos][0] += p4[Wpfb_pos][0];
7581 p4[Wp_pos][1] += p4[Wpfb_pos][1];
7582 p4[Wp_pos][2] += p4[Wpfb_pos][2];
7583 p4[Wp_pos][3] += p4[Wpfb_pos][3];
7586 p4[Wpf_pos][0] = momTmp->T()*GeV;
7587 p4[Wpf_pos][1] = momTmp->X()*GeV;
7588 p4[Wpf_pos][2] = momTmp->Y()*GeV;
7589 p4[Wpf_pos][3] = momTmp->Z()*GeV;
7590 p4[Wp_pos][0] += p4[Wpf_pos][0];
7591 p4[Wp_pos][1] += p4[Wpf_pos][1];
7592 p4[Wp_pos][2] += p4[Wpf_pos][2];
7593 p4[Wp_pos][3] += p4[Wpf_pos][3];
7595 MomStore[6] = MomStore[6] + (*momTmp);
7597 if (topDaughters.size()!=1){
for (
unsigned int ix=0; ix<4; ix++){
for (
unsigned int ip=b_pos; ip<=Wp_pos; ip++) p4[t_pos][ix] = p4[ip][ix]; } }
7601 const unsigned int tb_pos=3;
7602 const unsigned int bb_pos=5;
7603 const unsigned int Wm_pos=6;
7604 const unsigned int Wmf_pos=7;
7605 const unsigned int Wmfb_pos=8;
7606 for (
unsigned int ipar=0; ipar<antitopDaughters.size(); ipar++){
7607 TLorentzVector* momTmp = &(antitopDaughters.at(ipar).second);
7608 if (antitopDaughters.size()==1){
7609 p4[tb_pos][0] = momTmp->T()*GeV;
7610 p4[tb_pos][1] = momTmp->X()*GeV;
7611 p4[tb_pos][2] = momTmp->Y()*GeV;
7612 p4[tb_pos][3] = momTmp->Z()*GeV;
7616 p4[bb_pos][0] = momTmp->T()*GeV;
7617 p4[bb_pos][1] = momTmp->X()*GeV;
7618 p4[bb_pos][2] = momTmp->Y()*GeV;
7619 p4[bb_pos][3] = momTmp->Z()*GeV;
7622 p4[Wmf_pos][0] = momTmp->T()*GeV;
7623 p4[Wmf_pos][1] = momTmp->X()*GeV;
7624 p4[Wmf_pos][2] = momTmp->Y()*GeV;
7625 p4[Wmf_pos][3] = momTmp->Z()*GeV;
7626 p4[Wm_pos][0] += p4[Wmf_pos][0];
7627 p4[Wm_pos][1] += p4[Wmf_pos][1];
7628 p4[Wm_pos][2] += p4[Wmf_pos][2];
7629 p4[Wm_pos][3] += p4[Wmf_pos][3];
7632 p4[Wmfb_pos][0] = momTmp->T()*GeV;
7633 p4[Wmfb_pos][1] = momTmp->X()*GeV;
7634 p4[Wmfb_pos][2] = momTmp->Y()*GeV;
7635 p4[Wmfb_pos][3] = momTmp->Z()*GeV;
7636 p4[Wm_pos][0] += p4[Wmfb_pos][0];
7637 p4[Wm_pos][1] += p4[Wmfb_pos][1];
7638 p4[Wm_pos][2] += p4[Wmfb_pos][2];
7639 p4[Wm_pos][3] += p4[Wmfb_pos][3];
7641 MomStore[7] = MomStore[7] + (*momTmp);
7643 if (antitopDaughters.size()!=1){
for (
unsigned int ix=0; ix<4; ix++){
for (
unsigned int ip=5; ip<=6; ip++) p4[tb_pos][ix] = p4[ip][ix]; } }
7645 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
7646 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
7647 p4[2][0] += momTmp->T()*GeV;
7648 p4[2][1] += momTmp->X()*GeV;
7649 p4[2][2] += momTmp->Y()*GeV;
7650 p4[2][3] += momTmp->Z()*GeV;
7651 MomStore[5] = MomStore[5] + (*momTmp);
7655 for (
int ii=0; ii<13; ii++){
MELAout <<
"p4[" << ii <<
"] = ";
for (
int jj=0; jj<4; jj++)
MELAout << p4[ii][jj]/GeV <<
'\t';
MELAout << endl; }
7658 double defaultRenScale =
scale_.scale;
7659 double defaultFacScale =
facscale_.facscale;
7661 int defaultNflav =
nflav_.nflav;
7662 string defaultPdflabel =
pdlabel_.pdlabel;
7666 double alphasVal, alphasmzVal;
7676 <<
"TUtil::TTHiggsMatEl: Set AlphaS:\n"
7677 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
7679 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7684 for (
unsigned int ib1=0; ib1<3; ib1++){
7685 if (topDaughters.size()==1 && ib1!=0)
continue;
7687 unsigned int b1index;
7688 if (ib1==0) b1index = b_pos;
7689 else if (ib1==1) b1index = Wpf_pos;
7690 else b1index = Wpfb_pos;
7691 for (
unsigned int if1=0; if1<3; if1++){
7692 if (topDaughters.size()==1 && if1!=0)
continue;
7694 unsigned int f1index;
7695 if (if1==0) f1index = Wpf_pos;
7696 else if (if1==2) f1index = b_pos;
7697 else f1index = Wpfb_pos;
7698 for (
unsigned int ifb1=0; ifb1<3; ifb1++){
7699 if (topDaughters.size()==1 && ifb1!=0)
continue;
7701 unsigned int fb1index;
7702 if (ifb1==2) fb1index = Wpf_pos;
7703 else if (ifb1==1) fb1index = b_pos;
7704 else fb1index = Wpfb_pos;
7706 if (b1index==f1index || b1index==fb1index || f1index==fb1index)
continue;
7708 for (
unsigned int ib2=0; ib2<3; ib2++){
7709 if (antitopDaughters.size()==1 && ib2!=0)
continue;
7711 unsigned int b2index;
7712 if (ib2==0) b2index = bb_pos;
7713 else if (ib2==1) b2index = Wmf_pos;
7714 else b2index = Wmfb_pos;
7715 for (
unsigned int if2=0; if2<3; if2++){
7716 if (antitopDaughters.size()==1 && if2!=0)
continue;
7718 unsigned int f2index;
7719 if (if2==0) f2index = Wmf_pos;
7720 else if (if2==2) f2index = bb_pos;
7721 else f2index = Wmfb_pos;
7722 for (
unsigned int ifb2=0; ifb2<3; ifb2++){
7723 if (antitopDaughters.size()==1 && ifb2!=0)
continue;
7725 unsigned int fb2index;
7726 if (ifb2==2) fb2index = Wmf_pos;
7727 else if (ifb2==1) fb2index = bb_pos;
7728 else fb2index = Wmfb_pos;
7730 if (b2index==f2index || b2index==fb2index || f2index==fb2index)
continue;
7732 double p4_current[13][4]={ { 0 } };
7733 for (
unsigned int ix=0; ix<4; ix++){
7734 for (
unsigned int ip=0; ip<=2; ip++) p4_current[ip][ix] = p4[ip][ix];
7735 p4_current[t_pos][ix] = p4[t_pos][ix];
7736 p4_current[b1index][ix] = p4[b_pos][ix];
7737 p4_current[f1index][ix] = p4[Wpf_pos][ix];
7738 p4_current[fb1index][ix] = p4[Wpfb_pos][ix];
7739 p4_current[Wp_pos][ix] = p4_current[Wpf_pos][ix] + p4_current[Wpfb_pos][ix];
7741 p4_current[tb_pos][ix] = p4[tb_pos][ix];
7742 p4_current[b2index][ix] = p4[bb_pos][ix];
7743 p4_current[f2index][ix] = p4[Wmf_pos][ix];
7744 p4_current[fb2index][ix] = p4[Wmfb_pos][ix];
7745 p4_current[Wm_pos][ix] = p4_current[Wmf_pos][ix] + p4_current[Wmfb_pos][ix];
7749 <<
"TUtil::TTHiggsMatEl: Unswapped instance for "
7750 <<
"b(" << b_pos <<
") -> " << b1index <<
", "
7751 <<
"Wpf(" << Wpf_pos <<
") -> " << f1index <<
", "
7752 <<
"Wpfb(" << Wpfb_pos <<
") -> " << fb1index <<
", "
7753 <<
"bb(" << bb_pos <<
") -> " << b2index <<
", "
7754 <<
"Wmf(" << Wmf_pos <<
") -> " << f2index <<
", "
7755 <<
"Wmfb(" << Wmfb_pos <<
") -> " << fb2index << endl;
7756 for (
int ii=0; ii<13; ii++){
MELAout <<
"p4_instance[" << ii <<
"] = ";
for (
int jj=0; jj<4; jj++)
MELAout << p4_current[ii][jj]/GeV <<
'\t';
MELAout << endl; }
7760 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
7761 double MatElsq_tmp_swap[
nmsq][
nmsq]={ { 0 } };
7763 if (isUnknown[0] && isUnknown[1]){
7764 for (
unsigned int ix=0; ix<4; ix++){
7765 swap(p4_current[t_pos][ix], p4_current[tb_pos][ix]);
7766 swap(p4_current[b_pos][ix], p4_current[bb_pos][ix]);
7767 swap(p4_current[Wp_pos][ix], p4_current[Wm_pos][ix]);
7768 swap(p4_current[Wpf_pos][ix], p4_current[Wmf_pos][ix]);
7769 swap(p4_current[Wpfb_pos][ix], p4_current[Wmfb_pos][ix]);
7772 for (
int ix=0; ix<11; ix++){
for (
int iy=0; iy<11; iy++) MatElsq_tmp[iy][ix] = (MatElsq_tmp[iy][ix]+MatElsq_tmp_swap[iy][ix])/2.; }
7774 for (
int ix=0; ix<11; ix++){
for (
int iy=0; iy<11; iy++) MatElsq[iy][ix] += MatElsq_tmp[iy][ix]; }
7777 <<
"TUtil::TTHiggsMatEl: Swapped instance for "
7778 <<
"b(" << b_pos <<
") -> " << b1index <<
", "
7779 <<
"Wpf(" << Wpf_pos <<
") -> " << f1index <<
", "
7780 <<
"Wpfb(" << Wpfb_pos <<
") -> " << fb1index <<
", "
7781 <<
"bb(" << bb_pos <<
") -> " << b2index <<
", "
7782 <<
"Wmf(" << Wmf_pos <<
") -> " << f2index <<
", "
7783 <<
"Wmfb(" << Wmfb_pos <<
") -> " << fb2index << endl;
7784 for (
int ii=0; ii<13; ii++){
MELAout <<
"p4_instance[" << ii <<
"] = ";
for (
int jj=0; jj<4; jj++)
MELAout << p4_current[ii][jj]/GeV <<
'\t';
MELAout << endl; }
7794 int defaultTopDecay=-1;
7797 int GeVexponent_MEsq;
7798 if (topDecay>0) GeVexponent_MEsq = 4-(1+3*(nRequested_Tops+nRequested_Antitops))*2;
7799 else GeVexponent_MEsq = 4-(1+nRequested_Tops+nRequested_Antitops)*2;
7800 double constant = pow(GeV, -GeVexponent_MEsq);
7801 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7803 MELAout <<
"TUtil::TTHiggsMatEl: MEsq[ip][jp] = " << endl;
7804 for (
int iquark=-5; iquark<=5; iquark++){
7805 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << MatElsq[jquark+5][iquark+5] <<
'\t';
7809 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7813 <<
"TUtil::TTHiggsMatEl: Reset AlphaS:\n"
7814 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7816 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7820 <<
"TUtil::TTHiggsMatEl: Reset AlphaS result:\n"
7821 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7830 const double& EBEAM,
7834 const double GeV=1./100.;
7835 double sum_msqjk = 0;
7836 double MatElsq[
nmsq][
nmsq]={ { 0 } };
7837 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
7839 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::BBHiggsMatEl: Non-JHUGen MEs are not supported." << endl;
return sum_msqjk; }
7840 if (production!=
TVar::bbH){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::BBHiggsMatEl: Only bbH is supported." << endl;
return sum_msqjk; }
7843 int nRequested_AssociatedJets=2;
7856 <<
"TUtil::BBHiggsMatEl: Number of stable bs (" << mela_event.
pAssociated.size() <<
")"
7857 <<
" in bbH process is not 2!" << endl;
7863 bool isUnknown[2]; isUnknown[0]=
false; isUnknown[1]=
false;
7866 else antitopDaughters.push_back(mela_event.
pAssociated.at(0));
7868 else topDaughters.push_back(mela_event.
pAssociated.at(1));
7875 double p4[13][4]={ { 0 } };
7876 double MYIDUP_prod[2]={ 0 };
7877 TLorentzVector MomStore[
mxpart];
7878 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
7879 for (
int ipar=0; ipar<2; ipar++){
7880 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
7881 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
7883 else MYIDUP_prod[ipar] = 0;
7884 if (momTmp->T()>0.){
7885 p4[ipar][0] = -momTmp->T()*GeV;
7886 p4[ipar][1] = -momTmp->X()*GeV;
7887 p4[ipar][2] = -momTmp->Y()*GeV;
7888 p4[ipar][3] = -momTmp->Z()*GeV;
7889 MomStore[ipar] = (*momTmp);
7892 p4[ipar][0] = momTmp->T()*GeV;
7893 p4[ipar][1] = momTmp->X()*GeV;
7894 p4[ipar][2] = momTmp->Y()*GeV;
7895 p4[ipar][3] = momTmp->Z()*GeV;
7896 MomStore[ipar] = -(*momTmp);
7897 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
7902 for (
unsigned int ipar=0; ipar<topDaughters.size(); ipar++){
7903 TLorentzVector* momTmp = &(topDaughters.at(ipar).second);
7904 p4[4][0] = momTmp->T()*GeV;
7905 p4[4][1] = momTmp->X()*GeV;
7906 p4[4][2] = momTmp->Y()*GeV;
7907 p4[4][3] = momTmp->Z()*GeV;
7908 MomStore[6] = MomStore[6] + (*momTmp);
7912 for (
unsigned int ipar=0; ipar<antitopDaughters.size(); ipar++){
7913 TLorentzVector* momTmp = &(antitopDaughters.at(ipar).second);
7914 p4[3][0] = momTmp->T()*GeV;
7915 p4[3][1] = momTmp->X()*GeV;
7916 p4[3][2] = momTmp->Y()*GeV;
7917 p4[3][3] = momTmp->Z()*GeV;
7918 MomStore[7] = MomStore[7] + (*momTmp);
7921 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
7922 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
7923 p4[2][0] += momTmp->T()*GeV;
7924 p4[2][1] += momTmp->X()*GeV;
7925 p4[2][2] += momTmp->Y()*GeV;
7926 p4[2][3] += momTmp->Z()*GeV;
7927 MomStore[5] = MomStore[5] + (*momTmp);
7931 for (
int ii=0; ii<13; ii++){
MELAout <<
"p4[" << ii <<
"] = ";
for (
int jj=0; jj<4; jj++)
MELAout << p4[ii][jj]/GeV <<
'\t';
MELAout << endl; }
7934 double defaultRenScale =
scale_.scale;
7935 double defaultFacScale =
facscale_.facscale;
7938 int defaultNflav =
nflav_.nflav;
7939 string defaultPdflabel =
pdlabel_.pdlabel;
7943 double alphasVal, alphasmzVal;
7953 <<
"TUtil::BBHiggsMatEl: Set AlphaS:\n"
7954 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
7956 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7960 if (isUnknown[0] && isUnknown[1]){
7961 for (
unsigned int ix=0; ix<4; ix++) swap(p4[3][ix], p4[4][ix]);
7963 for (
int ix=0; ix<11; ix++){
for (
int iy=0; iy<11; iy++) MatElsq[iy][ix] = (MatElsq[iy][ix]+MatElsq_tmp[iy][ix])/2.; }
7966 int GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets)*2;
7967 double constant = pow(GeV, -GeVexponent_MEsq);
7968 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7970 MELAout <<
"TUtil::BBHiggsMatEl: MEsq[ip][jp] = " << endl;
7971 for (
int iquark=-5; iquark<=5; iquark++){
7972 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << MatElsq[jquark+5][iquark+5] <<
'\t';
7976 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7980 <<
"TUtil::BBHiggsMatEl: Reset AlphaS:\n"
7981 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7983 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7987 <<
"TUtil::BBHiggsMatEl: Reset AlphaS result:\n"
7988 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7996 MELAout <<
"TUtil::WipeMEArray: Initial MEsq:" << endl;
7997 for (
int iquark=-5; iquark<=5; iquark++){
7998 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
8010 for (
int iquark=-5; iquark<=5; iquark++){
8011 for (
int jquark=-5; jquark<=5; jquark++){
8017 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
8019 else msq[jquark+5][iquark+5]=0;
8029 if (
msq[5][5]>0.) nInstances=1;
8047 for (
int iquark=-5; iquark<=5; iquark++){
8048 for (
int jquark=-5; jquark<=5; jquark++){
8054 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
8056 else msq[jquark+5][iquark+5]=0;
8062 for (
int iquark=-5; iquark<=5; iquark++){
8063 for (
int jquark=-5; jquark<=5; jquark++){
8083 (TMath::Sign(1, iquark)==TMath::Sign(1, jquark))
8100 )
msq[jquark+5][iquark+5]=0;
8103 int order[2]={ -1, -1 };
8105 if (order[0]<0 || order[1]<0)
msq[jquark+5][iquark+5]=0;
8107 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
8109 else msq[jquark+5][iquark+5]=0;
8123 xx[0]>1. || xx[0]<=
xmin_.xmin
8125 xx[1]>1. || xx[1]<=
xmin_.xmin
8130 if (xx[0]>1. || xx[1]>1.)
MELAerr <<
"TUtil::CheckPartonMomFraction: At least one of the parton momentum fractions is greater than 1." << endl;
8131 else if (xx[0]<=
xmin_.xmin || xx[1]<=
xmin_.xmin)
MELAerr <<
"TUtil::CheckPartonMomFraction: At least one of the parton momentum fractions is less than or equal to " <<
xmin_.xmin <<
"." << endl;
8132 else MELAerr <<
"TUtil::CheckPartonMomFraction: EBEAM=" << EBEAM <<
"<=0." << endl;
8137 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::CheckPartonMomFraction: xx[0]: " << xx[0] <<
", xx[1] = " << xx[1] <<
", xmin = " <<
xmin_.xmin << endl;
8147 double fx1x2_jhu[2][13]={ { 0 } };
8148 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::ComputePDF: Calling setpdfs with xx[0]: " << xx[0] <<
", xx[1] = " << xx[1] << endl;
8151 for (
int ip=-6; ip<=6; ip++){
8154 if (ip!=0 && (abs(ip)%2==0)) fac=-1;
8155 else if (ip!=0) fac=+1;
8158 fx1[jp+5]=fx1x2_jhu[0][ip+6];
8159 fx2[jp+5]=fx1x2_jhu[1][ip+6];
8171 MELAout <<
"End TUtil::ComputePDF:"<< endl;
8172 for (
int ip=-
nf; ip<=
nf; ip++)
MELAout <<
"(fx1, fx2)[" << ip <<
"] = (" << fx1[ip+5] <<
" , " << fx2[ip+5] <<
")" << endl;
8178 double fx1[
nmsq]={ 0 };
8179 double fx2[
nmsq]={ 0 };
8182 ComputePDF(p0, p1, fx1, fx2, EBEAM, verbosity);
8197 const double GeV=1./100.;
8198 int isch=(int)scheme;
8199 double shat_jhu = pow(
sqrts*GeV, 2);
8221 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: code=" << code << endl;
8223 TLorentzVector
const nullFourVector(0, 0, 0, 0);
8226 vector<int> idVstar; idVstar.reserve(2);
8230 idVstar.push_back(melaCand->
id);
8231 idVstar.push_back(-9000);
8238 for (
unsigned char iv=0; iv<2; iv++){
8241 int idtmp = Vdau->
id;
8243 if (Vdau_i && Vdau_i->passSelection) daughters.push_back(
SimpleParticle_t(Vdau_i->id, Vdau_i->p4));
8247 idVstar.push_back(idtmp);
8250 else idVstar.push_back(-9000);
8253 if (daughters.size()>=2){
8254 size_t nffs = daughters.size()/2;
8255 for (
size_t iv=0; iv<nffs; iv++){
8257 daughters.at(2*iv+0).second, daughters.at(2*iv+0).first,
8258 daughters.at(2*iv+1).second, daughters.at(2*iv+1).first
8260 daughters.at(2*iv+0).second = corrPair.first;
8261 daughters.at(2*iv+1).second = corrPair.second;
8263 if (2*nffs<daughters.size()){
8264 TLorentzVector tmp = nullFourVector;
8267 lastDau.second, lastDau.first,
8270 lastDau.second = corrPair.first;
8280 int nsatisfied_jets=0;
8281 int nsatisfied_lnus=0;
8282 int nsatisfied_gammas=0;
8283 vector<MELAParticle*> candidateVs;
8286 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo!=0 case start" << endl;
8293 if ((abs(idV)==aVhypo || idV==0) && Vdau->getNDaughters()>0 && Vdau->passSelection){
8296 if (!Vdau_i){ doAdd=
false;
break; }
8302 !Vdau_i->passSelection
8308 if (doAdd) candidateVs.push_back(Vdau);
8312 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: candidateVs size = " << candidateVs.size() << endl;
8340 if (associated_tmp.size()>=2 || (associated_tmp.size()==1 &&
PDGHelpers::isAPhoton(associated_tmp.at(0).first))) idVstar.push_back(Vdau->id);
8342 if (associated_tmp.size()>=2){
8343 unsigned int nffs = associated_tmp.size()/2;
8344 for (
unsigned int iv=0; iv<nffs; iv++){
8346 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8347 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8349 associated_tmp.at(2*iv+0).second = corrPair.first;
8350 associated_tmp.at(2*iv+1).second = corrPair.second;
8352 if (2*nffs<associated_tmp.size()){
8353 TLorentzVector tmp = nullFourVector;
8355 associated_tmp.at(associated_tmp.size()-1).second, associated_tmp.at(associated_tmp.size()-1).first,
8358 associated_tmp.at(associated_tmp.size()-1).second = corrPair.first;
8364 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo!=0 case associated.size=" << associated.size() << endl;
8368 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 case begin" << endl;
8379 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 lep case associated_tmp.size=" << associated_tmp.size() << endl;
8382 if (associated_tmp.size()>=1){
8383 size_t nffs = associated_tmp.size()/2;
8384 for (
size_t iv=0; iv<nffs; iv++){
8385 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from lepton pair " << 2*iv+0 <<
'\t' << 2*iv+1 << endl;
8387 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8388 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8390 associated_tmp.at(2*iv+0).second = corrPair.first;
8391 associated_tmp.at(2*iv+1).second = corrPair.second;
8393 if (2*nffs<associated_tmp.size()){
8394 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from last lepton " << associated_tmp.size()-1 << endl;
8395 TLorentzVector tmp = nullFourVector;
8396 SimpleParticle_t& lastAssociated = associated_tmp.at(associated_tmp.size()-1);
8398 lastAssociated.second, lastAssociated.first,
8401 lastAssociated.second = corrPair.first;
8416 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 jet case associated_tmp.size=" << associated_tmp.size() << endl;
8419 if (associated_tmp.size()>=1){
8420 unsigned int nffs = associated_tmp.size()/2;
8421 for (
unsigned int iv=0; iv<nffs; iv++){
8422 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from jet pair " << 2*iv+0 <<
'\t' << 2*iv+1 << endl;
8424 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8425 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8427 associated_tmp.at(2*iv+0).second = corrPair.first;
8428 associated_tmp.at(2*iv+1).second = corrPair.second;
8430 if (2*nffs<associated_tmp.size()){
8431 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from last jet " << associated_tmp.size()-1 << endl;
8432 TLorentzVector tmp = nullFourVector;
8433 SimpleParticle_t& lastAssociated = associated_tmp.at(associated_tmp.size()-1);
8435 lastAssociated.second, lastAssociated.first,
8438 lastAssociated.second = corrPair.first;
8449 nsatisfied_gammas++;
8456 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: associated.size=" << associated.size() << endl;
8459 int nsatisfied_tops=0;
8460 int nsatisfied_antitops=0;
8461 vector<SimpleParticleCollection_t> topDaughters;
8462 vector<SimpleParticleCollection_t> antitopDaughters;
8466 vector<MELATopCandidate_t*> tops;
8467 vector<MELATopCandidate_t*> topbars;
8468 vector<MELATopCandidate_t*> unknowntops;
8473 if (theTop && theTop->passSelection){
8474 vector<MELATopCandidate_t*>* particleArray;
8475 if (theTop->id==6) particleArray = &tops;
8476 else if (theTop->id==-6) particleArray = &topbars;
8477 else particleArray = &unknowntops;
8482 ) particleArray->push_back(theTop);
8485 if (verbosity>=
TVar::DEBUG){
MELAout <<
"TUtil::GetBoostedParticleVectors: tops.size=" << tops.size() <<
", topbars.size=" << topbars.size() <<
", unknowntops.size=" << unknowntops.size() << endl; }
8505 if (vdaughters.size()==3) topDaughters.push_back(vdaughters);
8510 nsatisfied_antitops++;
8514 nsatisfied_antitops++;
8525 if (vdaughters.size()==3) antitopDaughters.push_back(vdaughters);
8538 nsatisfied_antitops++;
8554 if (vdaughters.size()==3) topDaughters.push_back(vdaughters);
8557 nsatisfied_antitops++;
8568 if (vdaughters.size()==3) antitopDaughters.push_back(vdaughters);
8575 MELAout <<
"TUtil::GetBoostedParticleVectors: stableTops.size=" << stableTops.size() << endl;
8576 MELAout <<
"TUtil::GetBoostedParticleVectors: stableAntitops.size=" << stableAntitops.size() << endl;
8577 MELAout <<
"TUtil::GetBoostedParticleVectors: topDaughters.size=" << topDaughters.size() << endl;
8578 for (
unsigned int itop=0; itop<topDaughters.size(); itop++)
MELAout <<
"TUtil::GetBoostedParticleVectors: topDaughters.at(" << itop <<
").size=" << topDaughters.at(itop).size() << endl;
8579 MELAout <<
"TUtil::GetBoostedParticleVectors: antitopDaughters.size=" << antitopDaughters.size() << endl;
8580 for (
unsigned int itop=0; itop<antitopDaughters.size(); itop++)
MELAout <<
"TUtil::GetBoostedParticleVectors: antitopDaughters.at(" << itop <<
").size=" << antitopDaughters.at(itop).size() << endl;
8585 TLorentzVector pTotal(0, 0, 0, 0);
8589 for (
SimpleParticle_t const& sp:stableAntitops) pTotal = pTotal + sp.second;
8594 double qX = pTotal.X();
8595 double qY = pTotal.Y();
8596 double qE = pTotal.T();
8597 if ((qX*qX+qY*qY)>0.){
8598 TVector3 boostV(-qX/qE, -qY/qE, 0.);
8605 pTotal.Boost(boostV);
8613 double sysPz= pTotal.Z();
8614 double sysE = pTotal.T();
8615 double pz0 = (sysE+sysPz)/2.;
8616 double pz1 = -(sysE-sysPz)/2.;
8619 int motherId[2]={ 0, 0 };
8621 for (
int ip=0; ip<2; ip++) motherId[ip]=melaCand->
getMother(ip)->
id;
8624 if (TMath::Sign(1., melaCand->
getMother(0)->
z()-melaCand->
getMother(1)->
z())!=TMath::Sign(1., pz0-pz1)){ swap(pz0, pz1); swap(E0, E1); }
8626 if ((motherId[0]<0 && motherId[1]>=0) || (motherId[1]>0 && motherId[0]<=0)){
8629 swap(motherId[0], motherId[1]);
8632 TLorentzVector pM[2];
8633 pM[0].SetXYZT(0., 0., pz0, E0);
8634 pM[1].SetXYZT(0., 0., pz1, E1);
8661 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.intermediateVid.size=" << mela_event.
intermediateVid.size() << endl;
8662 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pMothers.size=" << mela_event.
pMothers.size() << endl;
8663 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pDaughters.size=" << mela_event.
pDaughters.size() << endl;
8664 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pAssociated.size=" << mela_event.
pAssociated.size() << endl;
8665 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pStableTops.size=" << mela_event.
pStableTops.size() << endl;
8666 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pStableAntitops.size=" << mela_event.
pStableAntitops.size() << endl;
8667 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pTopDaughters.size=" << mela_event.
pTopDaughters.size() << endl;
8668 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pAntitopDaughters.size=" << mela_event.
pAntitopDaughters.size() << endl;
8669 MELAout <<
"End GetBoostedParticleVectors" << endl;
8682 std::vector<MELAParticle*>* particleList,
8683 std::vector<MELACandidate*>* candList
8687 if (!pDaughters){
MELAerr <<
"TUtil::ConvertVectorFormat: No daughters!" << endl;
return cand; }
8688 else if (pDaughters->size()==0){
MELAerr <<
"TUtil::ConvertVectorFormat: Daughter size==0!" << endl;
return cand; }
8689 else if (pDaughters->size()>4){
MELAerr <<
"TUtil::ConvertVectorFormat: Daughter size " << pDaughters->size() <<
">4 is not supported!" << endl;
return cand; }
8690 if (pMothers && pMothers->size()!=2){
MELAerr <<
"TUtil::ConvertVectorFormat: Mothers momentum size (" << pMothers->size() <<
") has to have had been 2! Continuing by omitting mothers." << endl; }
8693 std::vector<MELAParticle*> daughters;
8694 std::vector<MELAParticle*> aparticles;
8695 std::vector<MELAParticle*> mothers;
8696 for (
auto& spart:(*pDaughters)){
8699 if (particleList) particleList->push_back(onePart);
8700 daughters.push_back(onePart);
8703 for (
auto& spart:(*pAssociated)){
8706 if (particleList) particleList->push_back(onePart);
8707 aparticles.push_back(onePart);
8710 if (pMothers && pMothers->size()==2){
8711 for (
auto& spart:(*pMothers)){
8714 if (particleList) particleList->push_back(onePart);
8715 mothers.push_back(onePart);
8728 if (daughters.size()==1){
8733 else if (daughters.size()==2){
8736 TLorentzVector pH = F1->
p4+F2->
p4;
8747 else if (daughters.size()==3){
8761 TLorentzVector pH = F1->
p4+F2->
p4+gamma->
p4;
8775 TLorentzVector pH(0, 0, 0, 0);
8777 for (
unsigned char ip=0; ip<4; ip++){ pH = pH + (daughters.at(ip))->p4; charge += (daughters.at(ip))->charge(); }
8779 for (
unsigned char ip=0; ip<4; ip++) cand->
addDaughter(daughters.at(ip));
8788 MELAerr <<
"TUtil::ConvertVectorFormat: PDGHelpers::HDecayMode = " <<
PDGHelpers::HDecayMode <<
" is set incorrectly for Ndaughters=" << daughters.size() <<
". There is no automatic mechnism for this scenario. ";
8789 MELAerr <<
"Suggestion: Call PDGHelpers::setCandidateDecayMode(TVar::CandidateDecay_XY) before this call." << endl;
8797 if (!mothers.empty()){
8798 for (
auto& part:mothers) cand->
addMother(part);
8802 if (!aparticles.empty()){
8803 for (
auto& part:aparticles){
8804 const int& partId = part->id;
8813 if (candList && cand) candList->push_back(cand);
8823 std::vector<MELAParticle*>* particleList,
8824 std::vector<MELAThreeBodyDecayCandidate*>* tbdCandList
8828 if (!tbdDaughters){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: No daughters!" << endl;
return cand; }
8829 else if (tbdDaughters->empty()){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: Daughter size==0!" << endl;
return cand; }
8830 else if (!(tbdDaughters->size()==1 || tbdDaughters->size()==3)){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: Daughter size " << tbdDaughters->size() <<
"!=1 or 3 is not supported!" << endl;
return cand; }
8832 if (tbdDaughters->size()==1){
8833 if (abs((tbdDaughters->at(0)).first)==6 || (tbdDaughters->at(0)).first==0){
8835 tbdCandList->push_back(cand);
8838 else if (tbdDaughters->size()==3){
8843 if (Wf->
id<0 || Wfb->
id>0) std::swap(Wf, Wfb);
8845 particleList->push_back(partnerPart);
8846 particleList->push_back(Wf);
8847 particleList->push_back(Wfb);
8850 tbdCandList->push_back(cand);
8856 MELAout <<
"***** TUtil::PrintCandidateSummary *****" << endl;
8857 MELAout <<
"Candidate: " << cand << endl;
8862 MELAout <<
"***** TUtil::PrintCandidateSummary (Simple Event Record) *****" << endl;
8863 MELAout <<
"Candidate: " << cand << endl;
8873 for (
unsigned int ip=0; ip<cand->
pMothers.size(); ip++){
8876 <<
"\t\tV" << ip <<
" (" << part->first <<
") (X,Y,Z,T)=( "
8877 << part->second.X() <<
" , "
8878 << part->second.Y() <<
" , "
8879 << part->second.Z() <<
" , "
8880 << part->second.T() <<
" )" << endl;
8885 for (
unsigned int ip=0; ip<cand->
pDaughters.size(); ip++){
8888 <<
"\t\tDau[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8889 << part->second.X() <<
" , "
8890 << part->second.Y() <<
" , "
8891 << part->second.Z() <<
" , "
8892 << part->second.T() <<
" )" << endl;
8895 for (
unsigned int ip=0; ip<cand->
pAssociated.size(); ip++){
8898 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8899 << part->second.X() <<
" , "
8900 << part->second.Y() <<
" , "
8901 << part->second.Z() <<
" , "
8902 << part->second.T() <<
" )" << endl;
8905 for (
unsigned int ip=0; ip<cand->
pStableTops.size(); ip++){
8908 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8909 << part->second.X() <<
" , "
8910 << part->second.Y() <<
" , "
8911 << part->second.Z() <<
" , "
8912 << part->second.T() <<
" )" << endl;
8918 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8919 << part->second.X() <<
" , "
8920 << part->second.Y() <<
" , "
8921 << part->second.Z() <<
" , "
8922 << part->second.T() <<
" )" << endl;
8926 for (
unsigned int ip=0; ip<cand->
pTopDaughters.size(); ip++){
8927 MELAout <<
"\t\tTop[" << ip <<
"] daughters:" << endl;
8928 for (
unsigned int jp=0; jp<cand->
pTopDaughters.at(ip).size(); jp++){
8931 <<
"\t\t- Top daughter[" << ip << jp <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8932 << part->second.X() <<
" , "
8933 << part->second.Y() <<
" , "
8934 << part->second.Z() <<
" , "
8935 << part->second.T() <<
" )" << endl;
8940 MELAout <<
"\t\tAntitop[" << ip <<
"] daughters:" << endl;
8944 <<
"\t\t- Antitop daughter[" << ip << jp <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8945 << part->second.X() <<
" , "
8946 << part->second.Y() <<
" , "
8947 << part->second.Z() <<
" , "
8948 << part->second.T() <<
" )" << endl;