13 #include "TLorentzRotation.h"
40 TLorentzVector
const nullFourVector(0, 0, 0, 0);
41 TLorentzVector p1hat, p2hat;
42 if (p1==nullFourVector || p2==nullFourVector)
return;
45 TLorentzVector p12=p1+p2;
46 TLorentzVector diffp2p1=p2-p1;
47 double p1sq = p1.M2();
48 double p2sq = p2.M2();
49 double p1p2 = p1.Dot(p2);
50 double m1sq = m1*fabs(m1);
51 double m2sq = m2*fabs(m2);
52 double p12sq = p12.M2();
54 TLorentzVector avec=(p1sq*p2 - p2sq*p1 + p1p2*diffp2p1);
56 double b = (p12sq + m2sq - m1sq) * (pow(p1p2, 2) - p1sq*p2sq);
57 double c = pow((p12sq + m2sq - m1sq), 2)*p1sq/4. - pow((p1sq + p1p2), 2)*m2sq;
58 double eta = (-b - sqrt(fabs(b*b -4.*a*c)))/(2.*a);
59 double xi = (p12sq + m2sq - m1sq - 2.*eta*(p2sq + p1p2))/(2.*(p1sq + p1p2));
61 p1hat = (1.-xi)*p1 + (1.-eta)*p2;
62 p2hat = xi*p1 + eta*p2;
67 double energy, p3, newp3,
ratio;
68 energy = massiveJet.T();
70 newp3 = sqrt(max(pow(energy, 2)-mass*fabs(mass), 0.));
71 ratio = (p3>0. ? (newp3/p3) : 1.);
72 masslessJet.SetXYZT(massiveJet.X()*
ratio, massiveJet.Y()*
ratio, massiveJet.Z()*
ratio, energy);
75 TLorentzVector
const& jet1,
int const& jet1Id,
76 TLorentzVector
const& jet2,
int const& jet2Id,
79 TLorentzVector
const nullFourVector(0, 0, 0, 0);
80 TLorentzVector jet1massless(0, 0, 0, 0), jet2massless(0, 0, 0, 0);
129 pair<TLorentzVector, TLorentzVector> result(jet1massless, jet2massless);
138 daughters.at(1).second = corrWpair.first;
139 daughters.at(2).second = corrWpair.second;
141 TVector3 pW_boost_old = -pW.BoostVector();
149 TVector3 pW_boost_new = pW.BoostVector();
150 for (
unsigned int idau=1; idau<
daughters.size(); idau++){
151 daughters.at(idau).second.Boost(pW_boost_old);
152 daughters.at(idau).second.Boost(pW_boost_new);
157 TLorentzVector masslessRealJet(0, 0, 0, 0);
159 else masslessRealJet = realJet;
160 fakeJet = others + masslessRealJet;
161 fakeJet.SetVect(-fakeJet.Vect());
162 fakeJet.SetE(fakeJet.P());
166 std::pair<TLorentzVector, TLorentzVector>
TUtil::ComplexBoost(TVector3
const& beta, TLorentzVector
const& p4){
170 double b2 = bx*bx + by*by + bz*bz;
171 double gammasqinv = 1.-
b2;
174 double bp = bx*p4.X() + by*p4.Y() + bz*p4.Z();
175 double gammap_real=0;
176 double gammap_imag=0;
177 TLorentzVector p4new_real(0, 0, 0, 0), p4new_imag(0, 0, 0, 0);
179 gamma = 1./sqrt(gammasqinv);
180 if (
b2>0.) gammap_real = (gamma-1.)/
b2;
182 p4new_real.SetX(p4.X() + gammap_real*bp*bx + gamma*bx*p4.T());
183 p4new_real.SetY(p4.Y() + gammap_real*bp*by + gamma*by*p4.T());
184 p4new_real.SetZ(p4.Z() + gammap_real*bp*bz + gamma*bz*p4.T());
185 p4new_real.SetT(gamma*(p4.T() + bp));
187 else if (gammasqinv<0.){
188 gamma = -1./sqrt(-gammasqinv);
190 gammap_real = -1./
b2;
191 gammap_imag = gamma/
b2;
194 p4new_real.SetX(p4.X() + gammap_real*bp*bx);
195 p4new_real.SetY(p4.Y() + gammap_real*bp*by);
196 p4new_real.SetZ(p4.Z() + gammap_real*bp*bz);
198 p4new_imag.SetX(gammap_imag*bp*bx + gamma*bx*p4.T());
199 p4new_imag.SetY(gammap_imag*bp*by + gamma*by*p4.T());
200 p4new_imag.SetZ(gammap_imag*bp*bz + gamma*bz*p4.T());
201 p4new_imag.SetT(gamma*(p4.T() + bp));
204 return (pair<TLorentzVector, TLorentzVector>(p4new_real, p4new_imag));
214 TLorentzVector p4M11,
int Z1_lept1Id,
215 TLorentzVector p4M12,
int Z1_lept2Id,
216 TLorentzVector p4M21,
int Z2_lept1Id,
217 TLorentzVector p4M22,
int Z2_lept2Id
219 TLorentzVector
const nullFourVector(0, 0, 0, 0);
220 if (p4M12==nullFourVector || p4M22==nullFourVector){
222 p4M11 = f13Pair.first;
223 p4M21 = f13Pair.second;
225 else if (p4M11==nullFourVector || p4M21==nullFourVector){
227 p4M12 = f24Pair.first;
228 p4M22 = f24Pair.second;
233 p4M11 = f12Pair.first;
234 p4M12 = f12Pair.second;
235 p4M21 = f34Pair.first;
236 p4M22 = f34Pair.second;
240 TLorentzVector p4Z1 = p4M11 + p4M12;
241 TLorentzVector p4Z2 = p4M21 + p4M22;
248 (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0)
250 ((Z1_lept1Id*Z1_lept2Id>0 || (Z1_lept1Id==0 && Z1_lept2Id==0)) && p4M11.Phi()<=p4M12.Phi())
252 ) swap(p4M11, p4M12);
258 (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0)
260 ((Z2_lept1Id*Z2_lept2Id>0 || (Z2_lept1Id==0 && Z2_lept2Id==0)) && p4M21.Phi()<=p4M22.Phi())
262 ) swap(p4M21, p4M22);
267 TLorentzVector p4H = p4Z1 + p4Z2;
272 TVector3 boostX = -(p4H.BoostVector());
273 TLorentzVector thep4Z1inXFrame(p4Z1);
274 TLorentzVector thep4Z2inXFrame(p4Z2);
275 thep4Z1inXFrame.Boost(boostX);
276 thep4Z2inXFrame.Boost(boostX);
277 TVector3 theZ1X_p3 = thep4Z1inXFrame.Vect();
278 TVector3 theZ2X_p3 = thep4Z2inXFrame.Vect();
279 costhetastar = theZ1X_p3.CosTheta();
281 TVector3 boostV1(0, 0, 0);
282 TVector3 boostV2(0, 0, 0);
284 if (!(fabs(Z1_lept1Id)==21 || fabs(Z1_lept1Id)==22 || fabs(Z1_lept2Id)==21 || fabs(Z1_lept2Id)==22)){
285 boostV1 = -(p4Z1.BoostVector());
286 if (boostV1.Mag()>=1.) {
287 MELAout <<
"Warning: Mela::computeAngles: Z1 boost with beta=1, scaling down" << endl;
288 boostV1*=0.9999/boostV1.Mag();
290 TLorentzVector p4M11_BV1(p4M11);
291 TLorentzVector p4M12_BV1(p4M12);
292 TLorentzVector p4M21_BV1(p4M21);
293 TLorentzVector p4M22_BV1(p4M22);
294 p4M11_BV1.Boost(boostV1);
295 p4M12_BV1.Boost(boostV1);
296 p4M21_BV1.Boost(boostV1);
297 p4M22_BV1.Boost(boostV1);
299 TLorentzVector p4V2_BV1 = p4M21_BV1 + p4M22_BV1;
301 costheta1 = -p4V2_BV1.Vect().Unit().Dot(p4M11_BV1.Vect().Unit());
306 if (!(fabs(Z2_lept1Id)==21 || fabs(Z2_lept1Id)==22 || fabs(Z2_lept2Id)==21 || fabs(Z2_lept2Id)==22)){
307 boostV2 = -(p4Z2.BoostVector());
308 if (boostV2.Mag()>=1.) {
309 MELAout <<
"Warning: Mela::computeAngles: Z2 boost with beta=1, scaling down" << endl;
310 boostV2*=0.9999/boostV2.Mag();
312 TLorentzVector p4M11_BV2(p4M11);
313 TLorentzVector p4M12_BV2(p4M12);
314 TLorentzVector p4M21_BV2(p4M21);
315 TLorentzVector p4M22_BV2(p4M22);
316 p4M11_BV2.Boost(boostV2);
317 p4M12_BV2.Boost(boostV2);
318 p4M21_BV2.Boost(boostV2);
319 p4M22_BV2.Boost(boostV2);
321 TLorentzVector p4V1_BV2 = p4M11_BV2 + p4M12_BV2;
323 costheta2 = -p4V1_BV2.Vect().Unit().Dot(p4M21_BV2.Vect().Unit());
328 TLorentzVector p4M11_BX(p4M11);
329 TLorentzVector p4M12_BX(p4M12);
330 TLorentzVector p4M21_BX(p4M21);
331 TLorentzVector p4M22_BX(p4M22);
333 p4M11_BX.Boost(boostX);
334 p4M12_BX.Boost(boostX);
335 p4M21_BX.Boost(boostX);
336 p4M22_BX.Boost(boostX);
337 TLorentzVector p4V1_BX = p4M11_BX + p4M12_BX;
339 TVector3
const beamAxis(0, 0, 1);
340 TVector3 p3V1_BX = p4V1_BX.Vect().Unit();
341 TVector3 normal1_BX = (p4M11_BX.Vect().Cross(p4M12_BX.Vect())).Unit();
342 TVector3 normal2_BX = (p4M21_BX.Vect().Cross(p4M22_BX.Vect())).Unit();
343 TVector3 normalSC_BX = (beamAxis.Cross(p3V1_BX)).Unit();
347 float tmpSgnPhi = p3V1_BX.Dot(normal1_BX.Cross(normal2_BX));
349 if (fabs(tmpSgnPhi)>0.) sgnPhi = tmpSgnPhi/fabs(tmpSgnPhi);
350 float dot_BX12 = normal1_BX.Dot(normal2_BX);
351 if (fabs(dot_BX12)>=1.) dot_BX12 *= 1./fabs(dot_BX12);
352 Phi = sgnPhi * acos(-1.*dot_BX12);
356 float tmpSgnPhi1 = p3V1_BX.Dot(normal1_BX.Cross(normalSC_BX));
358 if (fabs(tmpSgnPhi1)>0.) sgnPhi1 = tmpSgnPhi1/fabs(tmpSgnPhi1);
359 float dot_BX1SC = normal1_BX.Dot(normalSC_BX);
360 if (fabs(dot_BX1SC)>=1.) dot_BX1SC *= 1./fabs(dot_BX1SC);
361 Phi1 = sgnPhi1 * acos(dot_BX1SC);
364 MELAout <<
"WARNING: NaN in computeAngles: "
365 << costhetastar <<
" "
369 << Phi1 <<
" " << endl;
370 MELAout <<
" boostV1: " <<boostV1.Pt() <<
" " << boostV1.Eta() <<
" " << boostV1.Phi() <<
" " << boostV1.Mag() << endl;
371 MELAout <<
" boostV2: " <<boostV2.Pt() <<
" " << boostV2.Eta() <<
" " << boostV2.Phi() <<
" " << boostV2.Mag() << endl;
381 TLorentzVector p4M11,
int Z1_lept1Id,
382 TLorentzVector p4M12,
int Z1_lept2Id,
383 TLorentzVector p4M21,
int Z2_lept1Id,
384 TLorentzVector p4M22,
int Z2_lept2Id
386 TLorentzVector
const nullFourVector(0, 0, 0, 0);
387 if (p4M12==nullFourVector || p4M22==nullFourVector){
389 p4M11 = f13Pair.first;
390 p4M21 = f13Pair.second;
392 else if (p4M11==nullFourVector || p4M21==nullFourVector){
394 p4M12 = f24Pair.first;
395 p4M22 = f24Pair.second;
400 p4M11 = f12Pair.first;
401 p4M12 = f12Pair.second;
402 p4M21 = f34Pair.first;
403 p4M22 = f34Pair.second;
406 TVector3 LabXaxis(1.0, 0.0, 0.0);
407 TVector3 LabYaxis(0.0, 1.0, 0.0);
408 TVector3 LabZaxis(0.0, 0.0, 1.0);
411 float Ebeam = sqrt(pbeam*pbeam + Mprot*Mprot);
413 TLorentzVector targ(0., 0., -pbeam, Ebeam);
414 TLorentzVector beam(0., 0., pbeam, Ebeam);
417 TLorentzVector p4Z1 = p4M11 + p4M12;
418 TLorentzVector p4Z2 = p4M21 + p4M22;
425 (Z1_lept1Id*Z1_lept2Id<0 && Z1_lept1Id<0)
427 ((Z1_lept1Id*Z1_lept2Id>0 || (Z1_lept1Id==0 && Z1_lept2Id==0)) && p4M11.Phi()<=p4M12.Phi())
429 ) swap(p4M11, p4M12);
435 (Z2_lept1Id*Z2_lept2Id<0 && Z2_lept1Id<0)
437 ((Z2_lept1Id*Z2_lept2Id>0 || (Z2_lept1Id==0 && Z2_lept2Id==0)) && p4M21.Phi()<=p4M22.Phi())
439 ) swap(p4M21, p4M22);
443 TLorentzVector p4H = p4Z1 + p4Z2;
444 TVector3 boostX = -(p4H.BoostVector());
453 TRotation rotationCS;
455 TLorentzVector beaminX(beam);
456 TLorentzVector targinX(targ);
457 targinX.Boost(boostX);
458 beaminX.Boost(boostX);
461 TVector3 beam_targ_bisecinX((beaminX.Vect().Unit() - targinX.Vect().Unit()).Unit());
464 TVector3 newZaxisCS(beam_targ_bisecinX.Unit());
465 TVector3 newYaxisCS(beaminX.Vect().Unit().Cross(newZaxisCS).Unit());
466 TVector3 newXaxisCS(newYaxisCS.Unit().Cross(newZaxisCS).Unit());
467 rotationCS.RotateAxes(newXaxisCS, newYaxisCS, newZaxisCS);
471 TLorentzVector thep4Z1inXFrame_rotCS(p4Z1);
472 TLorentzVector thep4Z2inXFrame_rotCS(p4Z2);
473 thep4Z1inXFrame_rotCS.Transform(rotationCS);
474 thep4Z2inXFrame_rotCS.Transform(rotationCS);
475 thep4Z1inXFrame_rotCS.Boost(boostX);
476 thep4Z2inXFrame_rotCS.Boost(boostX);
477 TVector3 theZ1XrotCS_p3 = TVector3(thep4Z1inXFrame_rotCS.X(), thep4Z1inXFrame_rotCS.Y(), thep4Z1inXFrame_rotCS.Z());
478 costhetastar = theZ1XrotCS_p3.CosTheta();
482 TLorentzVector p4M11_BX_rotCS(p4M11);
483 TLorentzVector p4M12_BX_rotCS(p4M12);
484 TLorentzVector p4M21_BX_rotCS(p4M21);
485 TLorentzVector p4M22_BX_rotCS(p4M22);
486 p4M11_BX_rotCS.Transform(rotationCS);
487 p4M12_BX_rotCS.Transform(rotationCS);
488 p4M21_BX_rotCS.Transform(rotationCS);
489 p4M22_BX_rotCS.Transform(rotationCS);
490 p4M11_BX_rotCS.Boost(boostX);
491 p4M12_BX_rotCS.Boost(boostX);
492 p4M21_BX_rotCS.Boost(boostX);
493 p4M22_BX_rotCS.Boost(boostX);
495 TLorentzVector p4Z1_BX_rotCS = p4M11_BX_rotCS + p4M12_BX_rotCS;
496 TVector3 p3V1_BX_rotCS = (p4Z1_BX_rotCS.Vect()).Unit();
497 TVector3
const beamAxis(0, 0, 1);
498 TVector3 normal1_BX_rotCS = (p4M11_BX_rotCS.Vect().Cross(p4M12_BX_rotCS.Vect())).Unit();
499 TVector3 normal2_BX_rotCS = (p4M21_BX_rotCS.Vect().Cross(p4M22_BX_rotCS.Vect())).Unit();
500 TVector3 normalSC_BX_rotCS = (beamAxis.Cross(p3V1_BX_rotCS)).Unit();
503 float tmpSgnPhi = p3V1_BX_rotCS.Dot(normal1_BX_rotCS.Cross(normal2_BX_rotCS));
505 if (fabs(tmpSgnPhi)>0.) sgnPhi = tmpSgnPhi/fabs(tmpSgnPhi);
506 float dot_BX12 = normal1_BX_rotCS.Dot(normal2_BX_rotCS);
507 if (fabs(dot_BX12)>=1.) dot_BX12 *= 1./fabs(dot_BX12);
508 Phi = sgnPhi * acos(-1.*dot_BX12);
511 float tmpSgnPhi1 = p3V1_BX_rotCS.Dot(normal1_BX_rotCS.Cross(normalSC_BX_rotCS));
513 if (fabs(tmpSgnPhi1)>0.) sgnPhi1 = tmpSgnPhi1/fabs(tmpSgnPhi1);
514 float dot_BX1SC = normal1_BX_rotCS.Dot(normalSC_BX_rotCS);
515 if (fabs(dot_BX1SC)>=1.) dot_BX1SC *= 1./fabs(dot_BX1SC);
516 Phi1 = sgnPhi1 * acos(dot_BX1SC);
519 if (!(fabs(Z1_lept1Id)==21 || fabs(Z1_lept1Id)==22 || fabs(Z1_lept2Id)==21 || fabs(Z1_lept2Id)==22)){
521 TRotation rotationZ1;
522 TVector3 newZaxisZ1(thep4Z1inXFrame_rotCS.Vect().Unit());
523 TVector3 newXaxisZ1(newYaxisCS.Cross(newZaxisZ1).Unit());
524 TVector3 newYaxisZ1(newZaxisZ1.Cross(newXaxisZ1).Unit());
525 rotationZ1.RotateAxes(newXaxisZ1, newYaxisZ1, newZaxisZ1);
528 TLorentzVector thep4Z1inXFrame_rotCS_rotZ1(thep4Z1inXFrame_rotCS);
529 thep4Z1inXFrame_rotCS_rotZ1.Transform(rotationZ1);
530 TVector3 boostZ1inX_rotCS_rotZ1= -(thep4Z1inXFrame_rotCS_rotZ1.BoostVector());
532 TLorentzVector p4M11_BX_rotCS_rotZ1(p4M11_BX_rotCS);
533 TLorentzVector p4M12_BX_rotCS_rotZ1(p4M12_BX_rotCS);
534 TLorentzVector p4M21_BX_rotCS_rotZ1(p4M21_BX_rotCS);
535 TLorentzVector p4M22_BX_rotCS_rotZ1(p4M22_BX_rotCS);
536 p4M11_BX_rotCS_rotZ1.Transform(rotationZ1);
537 p4M12_BX_rotCS_rotZ1.Transform(rotationZ1);
538 p4M21_BX_rotCS_rotZ1.Transform(rotationZ1);
539 p4M22_BX_rotCS_rotZ1.Transform(rotationZ1);
540 p4M11_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
541 p4M12_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
542 p4M21_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
543 p4M22_BX_rotCS_rotZ1.Boost(boostZ1inX_rotCS_rotZ1);
545 TLorentzVector p4V2_BX_rotCS_rotZ1 = p4M21_BX_rotCS_rotZ1 + p4M22_BX_rotCS_rotZ1;
547 costheta1 = -p4V2_BX_rotCS_rotZ1.Vect().Dot(p4M11_BX_rotCS_rotZ1.Vect())/p4V2_BX_rotCS_rotZ1.Vect().Mag()/p4M11_BX_rotCS_rotZ1.Vect().Mag();
552 if (!(fabs(Z2_lept1Id)==21 || fabs(Z2_lept1Id)==22 || fabs(Z2_lept2Id)==21 || fabs(Z2_lept2Id)==22)){
554 TRotation rotationZ2;
555 TVector3 newZaxisZ2(thep4Z2inXFrame_rotCS.Vect().Unit());
556 TVector3 newXaxisZ2(newYaxisCS.Cross(newZaxisZ2).Unit());
557 TVector3 newYaxisZ2(newZaxisZ2.Cross(newXaxisZ2).Unit());
558 rotationZ2.RotateAxes(newXaxisZ2, newYaxisZ2, newZaxisZ2);
561 TLorentzVector thep4Z2inXFrame_rotCS_rotZ2(thep4Z2inXFrame_rotCS);
562 thep4Z2inXFrame_rotCS_rotZ2.Transform(rotationZ2);
563 TVector3 boostZ2inX_rotCS_rotZ2= -(thep4Z2inXFrame_rotCS_rotZ2.BoostVector());
565 TLorentzVector p4M11_BX_rotCS_rotZ2(p4M11_BX_rotCS);
566 TLorentzVector p4M12_BX_rotCS_rotZ2(p4M12_BX_rotCS);
567 TLorentzVector p4M21_BX_rotCS_rotZ2(p4M21_BX_rotCS);
568 TLorentzVector p4M22_BX_rotCS_rotZ2(p4M22_BX_rotCS);
569 p4M11_BX_rotCS_rotZ2.Transform(rotationZ2);
570 p4M12_BX_rotCS_rotZ2.Transform(rotationZ2);
571 p4M21_BX_rotCS_rotZ2.Transform(rotationZ2);
572 p4M22_BX_rotCS_rotZ2.Transform(rotationZ2);
573 p4M11_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
574 p4M12_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
575 p4M21_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
576 p4M22_BX_rotCS_rotZ2.Boost(boostZ2inX_rotCS_rotZ2);
579 TLorentzVector p4V1_BX_rotCS_rotZ2= p4M11_BX_rotCS_rotZ2 + p4M12_BX_rotCS_rotZ2;
581 costheta2 = -p4V1_BX_rotCS_rotZ2.Vect().Dot(p4M21_BX_rotCS_rotZ2.Vect())/p4V1_BX_rotCS_rotZ2.Vect().Mag()/p4M21_BX_rotCS_rotZ2.Vect().Mag();
586 MELAout <<
"WARNING: NaN in computeAngles: "
587 << costhetastar <<
" "
591 << Phi1 <<
" " << endl;
604 TLorentzVector p4M11,
int Z1_lept1Id,
605 TLorentzVector p4M12,
int Z1_lept2Id,
606 TLorentzVector p4M21,
int Z2_lept1Id,
607 TLorentzVector p4M22,
int Z2_lept2Id,
608 TLorentzVector jet1,
int jet1Id,
609 TLorentzVector jet2,
int jet2Id,
610 TLorentzVector* injet1,
int injet1Id,
611 TLorentzVector* injet2,
int injet2Id
613 TLorentzVector
const nullFourVector(0, 0, 0, 0);
615 TLorentzVector jet1massless, jet2massless;
617 jet1massless = jetPair.first;
618 jet2massless = jetPair.second;
621 TLorentzVector p4Z1 = p4M11 + p4M12;
622 TLorentzVector p4Z2 = p4M21 + p4M22;
623 TLorentzVector pH = p4Z1+p4Z2;
625 if (jet1massless.Z() < jet2massless.Z()) { swap(jet1massless, jet2massless); swap(jet1Id, jet2Id); }
631 TLorentzRotation movingframe;
632 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
633 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
634 movingframe.Boost(-pHJJ_perp.BoostVector());
635 pHJJ.Boost(-pHJJ_perp.BoostVector());
636 movingframe.Boost(-pHJJ.BoostVector());
637 pHJJ.Boost(-pHJJ.BoostVector());
639 TLorentzVector P1(0, 0, pHJJ.T()/2, pHJJ.T()/2);
640 TLorentzVector P2(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
642 P1.Transform(movingframe.Inverse());
643 P2.Transform(movingframe.Inverse());
644 pHJJ.Transform(movingframe.Inverse());
647 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
650 if (P1.Z() < P2.Z()){
652 swap(injet1Id, injet2Id);
655 int diff1Id = jet1Id-injet1Id;
656 int diff2Id = jet2Id-injet2Id;
659 (diff1Id==0 && diff2Id==0 && !(injet1Id==21 || injet2Id==21))
661 ((fabs(diff1Id)==1 || fabs(diff1Id)==3 || fabs(diff1Id)==5) && (fabs(diff2Id)==1 || fabs(diff2Id)==3 || fabs(diff2Id)==5))
664 int diff12Id = jet1Id-injet2Id;
665 int diff21Id = jet2Id-injet1Id;
667 ((diff12Id==0 || diff21Id==0) && !(injet1Id==21 || injet2Id==21))
669 ((fabs(diff12Id)==1 || fabs(diff12Id)==3 || fabs(diff12Id)==5) || (fabs(diff21Id)==1 || fabs(diff21Id)==3 || fabs(diff21Id)==5))
672 swap(injet1Id, injet2Id);
677 TLorentzRotation ZZframe;
678 ZZframe.Boost(-pH.BoostVector());
679 P1.Transform(ZZframe);
680 P2.Transform(ZZframe);
681 p4Z1.Transform(ZZframe);
682 p4Z2.Transform(ZZframe);
683 jet1massless.Transform(ZZframe);
684 jet2massless.Transform(ZZframe);
686 TLorentzVector fermion1, fermion2, antifermion1, antifermion2;
689 if (injet1 && injet1Id<0){
691 antifermion1 = jet1massless;
694 fermion1 = jet1massless;
697 if (injet2 && injet2Id<0){
699 antifermion2 = jet2massless;
702 fermion2 = jet2massless;
707 TLorentzVector V1 = fermion1 + antifermion1;
708 TLorentzVector V2 = fermion2 + antifermion2;
709 TVector3 normvec1 = fermion1.Vect().Cross(antifermion1.Vect()).Unit();
710 TVector3 normvec2 = fermion2.Vect().Cross(antifermion2.Vect()).Unit();
711 TVector3 normvec3 = p4Z2.Vect().Cross(V1.Vect()).Unit();
712 double cosPhi = normvec1.Dot(normvec2);
713 double sgnPhi = normvec1.Cross(normvec2).Dot(V1.Vect());
714 if (fabs(sgnPhi)>0.) sgnPhi = sgnPhi/fabs(sgnPhi);
715 double cosPhi1 = normvec1.Dot(normvec3);
716 double sgnPhi1 = normvec1.Cross(normvec3).Dot(V1.Vect());
717 if (fabs(sgnPhi1)>0.) sgnPhi1 = sgnPhi1/fabs(sgnPhi1);
718 if (fabs(cosPhi)>1) cosPhi *= 1./fabs(cosPhi);
719 if (fabs(cosPhi1)>1) cosPhi1 *= 1./fabs(cosPhi1);
720 Phi = acos(-cosPhi)*sgnPhi;
721 Phi1 = acos(cosPhi1)*sgnPhi1;
722 costhetastar = V1.Vect().Unit().Dot(p4Z2.Vect().Unit());
727 costheta1 = V1.Vect().Unit().Dot(fermion1.Vect().Unit());
728 costheta2 = V2.Vect().Unit().Dot(fermion2.Vect().Unit());
732 float& costheta1_real,
float& costheta1_imag,
733 float& costheta2_real,
float& costheta2_imag,
738 TLorentzVector p4M11,
int Z1_lept1Id,
739 TLorentzVector p4M12,
int Z1_lept2Id,
740 TLorentzVector p4M21,
int Z2_lept1Id,
741 TLorentzVector p4M22,
int Z2_lept2Id,
742 TLorentzVector jet1,
int jet1Id,
743 TLorentzVector jet2,
int jet2Id,
744 TLorentzVector* injet1,
int injet1Id,
745 TLorentzVector* injet2,
int injet2Id
747 TLorentzVector
const nullFourVector(0, 0, 0, 0);
749 TLorentzVector jet1massless, jet2massless;
751 jet1massless = jetPair.first;
752 jet2massless = jetPair.second;
755 TLorentzVector p4Z1 = p4M11 + p4M12;
756 TLorentzVector p4Z2 = p4M21 + p4M22;
757 TLorentzVector pH = p4Z1+p4Z2;
759 if (jet1massless.Z() < jet2massless.Z()) { swap(jet1massless, jet2massless); swap(jet1Id, jet2Id); }
765 TLorentzRotation movingframe;
766 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
767 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
768 movingframe.Boost(-pHJJ_perp.BoostVector());
769 pHJJ.Boost(-pHJJ_perp.BoostVector());
770 movingframe.Boost(-pHJJ.BoostVector());
771 pHJJ.Boost(-pHJJ.BoostVector());
773 TLorentzVector P1(0, 0, pHJJ.T()/2, pHJJ.T()/2);
774 TLorentzVector P2(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
776 P1.Transform(movingframe.Inverse());
777 P2.Transform(movingframe.Inverse());
778 pHJJ.Transform(movingframe.Inverse());
781 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
784 if (P1.Z() < P2.Z()){
786 swap(injet1Id, injet2Id);
789 int diff1Id = jet1Id-injet1Id;
790 int diff2Id = jet2Id-injet2Id;
793 (diff1Id==0 && diff2Id==0 && !(injet1Id==21 || injet2Id==21))
795 ((fabs(diff1Id)==1 || fabs(diff1Id)==3 || fabs(diff1Id)==5) && (fabs(diff2Id)==1 || fabs(diff2Id)==3 || fabs(diff2Id)==5))
798 int diff12Id = jet1Id-injet2Id;
799 int diff21Id = jet2Id-injet1Id;
801 ((diff12Id==0 || diff21Id==0) && !(injet1Id==21 || injet2Id==21))
803 ((fabs(diff12Id)==1 || fabs(diff12Id)==3 || fabs(diff12Id)==5) || (fabs(diff21Id)==1 || fabs(diff21Id)==3 || fabs(diff21Id)==5))
806 swap(injet1Id, injet2Id);
811 TLorentzRotation ZZframe;
812 ZZframe.Boost(-pH.BoostVector());
813 P1.Transform(ZZframe);
814 P2.Transform(ZZframe);
815 p4Z1.Transform(ZZframe);
816 p4Z2.Transform(ZZframe);
817 jet1massless.Transform(ZZframe);
818 jet2massless.Transform(ZZframe);
820 TLorentzVector fermion1, fermion2, antifermion1, antifermion2;
823 if (injet1 && injet1Id<0){
825 antifermion1 = jet1massless;
828 fermion1 = jet1massless;
831 if (injet2 && injet2Id<0){
833 antifermion2 = jet2massless;
836 fermion2 = jet2massless;
841 TLorentzVector V1 = fermion1 + antifermion1;
842 TLorentzVector V2 = fermion2 + antifermion2;
843 TVector3 normvec1 = fermion1.Vect().Cross(antifermion1.Vect()).Unit();
844 TVector3 normvec2 = fermion2.Vect().Cross(antifermion2.Vect()).Unit();
845 TVector3 normvec3 = p4Z2.Vect().Cross(V1.Vect()).Unit();
846 double cosPhi = normvec1.Dot(normvec2);
847 double sgnPhi = normvec1.Cross(normvec2).Dot(V1.Vect());
848 if (fabs(sgnPhi)>0.) sgnPhi = sgnPhi/fabs(sgnPhi);
849 double cosPhi1 = normvec1.Dot(normvec3);
850 double sgnPhi1 = normvec1.Cross(normvec3).Dot(V1.Vect());
851 if (fabs(sgnPhi1)>0.) sgnPhi1 = sgnPhi1/fabs(sgnPhi1);
852 if (fabs(cosPhi)>1) cosPhi *= 1./fabs(cosPhi);
853 if (fabs(cosPhi1)>1) cosPhi1 *= 1./fabs(cosPhi1);
854 Phi = acos(-cosPhi)*sgnPhi;
855 Phi1 = acos(cosPhi1)*sgnPhi1;
856 costhetastar = V1.Vect().Unit().Dot(p4Z2.Vect().Unit());
864 pair<TLorentzVector, TLorentzVector> fermion1_BV1 =
TUtil::ComplexBoost(V1.BoostVector(), fermion1);
865 costheta1_real = -(V2_BV1.first.Vect().Unit().Dot(fermion1_BV1.first.Vect().Unit()) - V2_BV1.second.Vect().Unit().Dot(fermion1_BV1.second.Vect().Unit()));
866 costheta1_imag = -(V2_BV1.first.Vect().Unit().Dot(fermion1_BV1.second.Vect().Unit()) + V2_BV1.second.Vect().Unit().Dot(fermion1_BV1.first.Vect().Unit()));
869 pair<TLorentzVector, TLorentzVector> fermion2_BV2 =
TUtil::ComplexBoost(V2.BoostVector(), fermion2);
870 costheta2_real = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()) - V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()));
871 costheta2_imag = -(V1_BV2.first.Vect().Unit().Dot(fermion2_BV2.second.Vect().Unit()) + V1_BV2.second.Vect().Unit().Dot(fermion2_BV2.first.Vect().Unit()));
881 TLorentzVector p4M11,
int Z1_lept1Id,
882 TLorentzVector p4M12,
int Z1_lept2Id,
883 TLorentzVector p4M21,
int Z2_lept1Id,
884 TLorentzVector p4M22,
int Z2_lept2Id,
885 TLorentzVector jet1,
int jet1Id,
886 TLorentzVector jet2,
int jet2Id,
887 TLorentzVector* injet1,
int injet1Id,
888 TLorentzVector* injet2,
int injet2Id
890 TLorentzVector
const nullFourVector(0, 0, 0, 0);
892 TLorentzVector jet1massless, jet2massless;
894 jet1massless = jetPair.first;
895 jet2massless = jetPair.second;
898 TLorentzVector p4Z1 = p4M11 + p4M12;
899 TLorentzVector p4Z2 = p4M21 + p4M22;
900 TLorentzVector pH = p4Z1 + p4Z2;
904 (jet1Id*jet2Id<0 && jet1Id<0)
906 ((jet1Id*jet2Id>0 || (jet1Id==0 && jet2Id==0)) && jet1massless.Phi()<=jet2massless.Phi())
908 swap(jet1massless, jet2massless);
909 swap(jet1Id, jet2Id);
917 TLorentzRotation movingframe;
918 TLorentzVector pHJJ = pH+jet1massless+jet2massless;
919 TLorentzVector pHJJ_perp(pHJJ.X(), pHJJ.Y(), 0, pHJJ.T());
920 movingframe.Boost(-pHJJ_perp.BoostVector());
921 pHJJ.Boost(-pHJJ_perp.BoostVector());
922 movingframe.Boost(-pHJJ.BoostVector());
923 pHJJ.Boost(-pHJJ.BoostVector());
925 TLorentzVector P1(0, 0, -pHJJ.T()/2, pHJJ.T()/2);
926 TLorentzVector P2(0, 0, pHJJ.T()/2, pHJJ.T()/2);
928 P1.Transform(movingframe.Inverse());
929 P2.Transform(movingframe.Inverse());
930 pHJJ.Transform(movingframe.Inverse());
933 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pHJJ.P())<pHJJ.P()*1e-4){
938 (injet1Id*injet2Id<0 && injet1Id>0)
940 (injet1Id*injet2Id>0 && P1.Z()>=P2.Z())
943 swap(injet1Id, injet2Id);
948 TLorentzRotation ZZframe;
949 TVector3
const beamAxis(0, 0, 1);
950 if (p4Z1==nullFourVector || p4Z2==nullFourVector){
951 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
952 if (pNewAxis != nullFourVector.Vect()){
953 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
954 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
956 P1.Transform(ZZframe);
957 P2.Transform(ZZframe);
958 jet1massless = -jet1massless; jet1massless.Transform(ZZframe); jet1massless = -jet1massless;
959 jet2massless = -jet2massless; jet2massless.Transform(ZZframe); jet2massless = -jet2massless;
963 ZZframe.Boost(-pH.BoostVector());
964 p4Z1.Boost(-pH.BoostVector());
965 p4Z2.Boost(-pH.BoostVector());
966 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
967 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
968 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
969 P1.Transform(ZZframe);
970 P2.Transform(ZZframe);
971 jet1massless = -jet1massless; jet1massless.Transform(ZZframe); jet1massless = -jet1massless;
972 jet2massless = -jet2massless; jet2massless.Transform(ZZframe); jet2massless = -jet2massless;
987 m2 = (jet1massless + jet2massless).M();
1011 TLorentzVector p4M11,
int Z1_lept1Id,
1012 TLorentzVector p4M12,
int Z1_lept2Id,
1013 TLorentzVector p4M21,
int Z2_lept1Id,
1014 TLorentzVector p4M22,
int Z2_lept2Id,
1015 TLorentzVector b,
int bId,
1016 TLorentzVector Wplusf,
int WplusfId,
1017 TLorentzVector Wplusfb,
int WplusfbId,
1018 TLorentzVector bbar,
int bbarId,
1019 TLorentzVector Wminusf,
int WminusfId,
1020 TLorentzVector Wminusfb,
int WminusfbId,
1021 TLorentzVector* injet1,
int injet1Id,
1022 TLorentzVector* injet2,
int injet2Id
1024 TLorentzVector
const nullFourVector(0, 0, 0, 0);
1025 TVector3
const beamAxis(0, 0, 1);
1026 TVector3
const xAxis(1, 0, 0);
1029 hs=hincoming=hTT=PhiTT=Phi1
1030 =hbb=hWW=Phibb=Phi1bb
1032 =hWminusf=PhiWminusf=0;
1038 float m1_dummy=0, m2_dummy=0;
1040 hs, hincoming, hTT, PhiTT, Phi1,
1055 WplusfId!=-9000 && WplusfbId!=-9000
1058 (WplusfId*WplusfbId<0 && WplusfId<0)
1060 ((WplusfId*WplusfbId>0 || (WplusfId==0 && WplusfbId==0)) && Wplusf.Phi()<=Wplusfb.Phi())
1063 swap(Wplusf, Wplusfb);
1064 swap(WplusfId, WplusfbId);
1068 WminusfId!=-9000 && WminusfbId!=-9000
1071 (WminusfId*WminusfbId<0 && WminusfId<0)
1073 ((WminusfId*WminusfbId>0 || (WminusfId==0 && WminusfbId==0)) && Wminusf.Phi()<=Wminusfb.Phi())
1076 swap(Wminusf, Wminusfb);
1077 swap(WminusfId, WminusfbId);
1081 if (bId!=-9000 && WplusfId!=-9000 && WplusfbId!=-9000){
1083 tmp_daus.reserve(3);
1084 tmp_daus.emplace_back(bId, b);
1085 tmp_daus.emplace_back(WplusfId, Wplusf);
1086 tmp_daus.emplace_back(WplusfbId, Wplusfb);
1088 b = tmp_daus.at(0).second;
1089 Wplusf = tmp_daus.at(1).second;
1090 Wplusfb = tmp_daus.at(2).second;
1092 else if (WplusfId!=-9000 && WplusfbId!=-9000){
1094 Wplusf = Wffb.first;
1095 Wplusfb = Wffb.second;
1098 if (bbarId!=-9000 && WminusfId!=-9000 && WminusfbId!=-9000){
1100 tmp_daus.reserve(3);
1101 tmp_daus.emplace_back(bbarId, bbar);
1102 tmp_daus.emplace_back(WminusfId, Wminusf);
1103 tmp_daus.emplace_back(WminusfbId, Wminusfb);
1105 bbar = tmp_daus.at(0).second;
1106 Wminusf = tmp_daus.at(1).second;
1107 Wminusfb = tmp_daus.at(2).second;
1109 else if (WminusfId!=-9000 && WminusfbId!=-9000){
1111 Wminusf = Wffb.first;
1112 Wminusfb = Wffb.second;
1117 TLorentzVector p4Z1 = p4M11 + p4M12;
1118 TLorentzVector p4Z2 = p4M21 + p4M22;
1120 TLorentzVector pH = p4Z1 + p4Z2;
1122 TLorentzVector pWplus = Wplusf + Wplusfb;
1123 TLorentzVector pTop = b + pWplus;
1125 TLorentzVector pWminus = Wminusf + Wminusfb;
1126 TLorentzVector pATop = bbar + pWminus;
1128 TLorentzVector pTT = pTop + pATop;
1129 TLorentzVector pTTH = pTT + pH;
1136 TLorentzRotation movingframe;
1137 TLorentzVector pTTH_perp(pTTH.X(), pTTH.Y(), 0, pTTH.T());
1138 movingframe.Boost(-pTTH_perp.BoostVector());
1139 pTTH.Boost(-pTTH_perp.BoostVector());
1140 movingframe.Boost(-pTTH.BoostVector());
1141 pTTH.Boost(-pTTH.BoostVector());
1143 TLorentzVector P1(0, 0, -pTTH.T()/2, pTTH.T()/2);
1144 TLorentzVector P2(0, 0, pTTH.T()/2, pTTH.T()/2);
1146 P1.Transform(movingframe.Inverse());
1147 P2.Transform(movingframe.Inverse());
1148 pTTH.Transform(movingframe.Inverse());
1151 if (injet1 && injet2 && fabs((*injet1+*injet2).P()-pTTH.P())<pTTH.P()*1e-4){
1156 (injet1Id*injet2Id<0 && injet1Id>0)
1158 (injet1Id*injet2Id>0 && P1.Z()>=P2.Z())
1161 swap(injet1Id, injet2Id);
1166 TLorentzRotation ZZframe;
1167 bool applyZZframe=
false;
1168 if (p4Z1==nullFourVector || p4Z2==nullFourVector){
1169 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
1170 if (pNewAxis != nullFourVector.Vect()){
1171 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
1172 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
1177 TVector3 pHboost = pH.BoostVector();
1178 ZZframe.Boost(-pHboost);
1179 p4Z1.Boost(-pHboost);
1180 p4Z2.Boost(-pHboost);
1181 TVector3 pNewAxis = (p4Z2-p4Z1).Vect().Unit();
1182 TVector3 pNewAxisPerp = pNewAxis.Cross(beamAxis);
1183 ZZframe.Rotate(acos(pNewAxis.Dot(beamAxis)), pNewAxisPerp);
1185 p4Z1.Boost(pHboost);
1186 p4Z2.Boost(pHboost);
1190 P1.Transform(ZZframe); P2.Transform(ZZframe);
1191 p4Z1 = -p4Z1; p4Z1.Transform(ZZframe); p4Z1 = -p4Z1;
1192 p4Z2 = -p4Z2; p4Z2.Transform(ZZframe); p4Z2 = -p4Z2;
1193 b = -b; b.Transform(ZZframe); b = -b;
1194 Wplusf = -Wplusf; Wplusf.Transform(ZZframe); Wplusf = -Wplusf;
1195 Wplusfb = -Wplusfb; Wplusfb.Transform(ZZframe); Wplusfb = -Wplusfb;
1196 bbar = -bbar; bbar.Transform(ZZframe); bbar = -bbar;
1197 Wminusf = -Wminusf; Wminusf.Transform(ZZframe); Wminusf = -Wminusf;
1198 Wminusfb = -Wminusfb; Wminusfb.Transform(ZZframe); Wminusfb = -Wminusfb;
1202 pWplus = Wplusf + Wplusfb;
1204 pWminus = Wminusf + Wminusfb;
1205 pATop = bbar + pWminus;
1224 TVector3 pHboost = pH.BoostVector();
1232 Wplusf.Boost(-pHboost);
1233 Wplusfb.Boost(-pHboost);
1234 bbar.Boost(-pHboost);
1235 Wminusf.Boost(-pHboost);
1236 Wminusfb.Boost(-pHboost);
1241 pWplus = Wplusf + Wplusfb;
1243 pWminus = Wminusf + Wminusfb;
1244 pATop = bbar + pWminus;
1251 TLorentzRotation TTframe;
1254 TVector3 nNewZAxis = (-P1-P2-pTop-pATop).Vect().Unit();
1255 if (nNewZAxis != nullFourVector.Vect()){
1256 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1257 TTframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1260 P1.Transform(TTframe);
1261 P2.Transform(TTframe);
1262 b.Transform(TTframe);
1263 Wplusf.Transform(TTframe);
1264 Wplusfb.Transform(TTframe);
1265 bbar.Transform(TTframe);
1266 Wminusf.Transform(TTframe);
1267 Wminusfb.Transform(TTframe);
1270 pWplus = Wplusf + Wplusfb;
1272 pWminus = Wminusf + Wminusfb;
1273 pATop = bbar + pWminus;
1279 int id_dummy_Wplus=24;
1280 int id_dummy_Wminus=-24;
1281 if (WplusfId!=-9000 && WplusfbId!=-9000) id_dummy_Wplus=-9000;
1282 if (WminusfId!=-9000 && WminusfbId!=-9000) id_dummy_Wminus=-9000;
1291 pWplus, id_dummy_Wplus,
1292 pWminus, id_dummy_Wminus
1297 if (WplusfId!=-9000 && WplusfbId!=-9000){
1298 TLorentzVector pATop_tmp = pATop;
1299 TVector3 pWplus_boost = (Wplusf+Wplusfb).BoostVector();
1300 b.Boost(-pWplus_boost);
1301 Wplusf.Boost(-pWplus_boost);
1302 Wplusfb.Boost(-pWplus_boost);
1303 pATop_tmp.Boost(-pWplus_boost);
1305 TLorentzRotation Wplusframe;
1308 TVector3 nNewZAxis = (b+Wplusf+Wplusfb-pATop_tmp).Vect().Unit();
1309 if (nNewZAxis != nullFourVector.Vect()){
1310 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1311 Wplusframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1314 b.Transform(Wplusframe);
1315 Wplusf.Transform(Wplusframe);
1316 Wplusfb.Transform(Wplusframe);
1318 TVector3 n3_Wb = -b.Vect().Unit();
1319 TVector3 nW = ((Wplusf-Wplusfb).Vect().Cross(n3_Wb)).Unit();
1320 TVector3 nS = (beamAxis.Cross(n3_Wb)).Unit();
1323 hWplusf = n3_Wb.Dot(Wplusf.Vect().Unit());
1325 float tmpSgnPhiWplusf = n3_Wb.Dot(nW.Cross(nS));
1326 float sgnPhiWplusf = 0;
1327 if (fabs(tmpSgnPhiWplusf)>0.) sgnPhiWplusf = tmpSgnPhiWplusf/fabs(tmpSgnPhiWplusf);
1328 float dot_nWnS = nW.Dot(nS);
1329 if (fabs(dot_nWnS)>=1.) dot_nWnS *= 1./fabs(dot_nWnS);
1330 PhiWplusf = sgnPhiWplusf * acos(dot_nWnS);
1339 if (WminusfId!=-9000 && WminusfbId!=-9000){
1340 TLorentzVector pTop_tmp = pTop;
1341 TVector3 pWminus_boost = (Wminusf+Wminusfb).BoostVector();
1342 bbar.Boost(-pWminus_boost);
1343 Wminusf.Boost(-pWminus_boost);
1344 Wminusfb.Boost(-pWminus_boost);
1345 pTop_tmp.Boost(-pWminus_boost);
1347 TLorentzRotation Wminusframe;
1350 TVector3 nNewZAxis = (pTop_tmp-bbar-Wminusf-Wminusfb).Vect().Unit();
1351 if (nNewZAxis != nullFourVector.Vect()){
1352 TVector3 nNewZAxisPerp = nNewZAxis.Cross(beamAxis);
1353 Wminusframe.Rotate(acos(nNewZAxis.Dot(beamAxis)), nNewZAxisPerp);
1356 bbar.Transform(Wminusframe);
1357 Wminusf.Transform(Wminusframe);
1358 Wminusfb.Transform(Wminusframe);
1360 TVector3 n3_Wb = -bbar.Vect().Unit();
1361 TVector3 nW = ((Wminusf-Wminusfb).Vect().Cross(n3_Wb)).Unit();
1362 TVector3 nS = (beamAxis.Cross(n3_Wb)).Unit();
1365 hWminusf = n3_Wb.Dot(Wminusf.Vect().Unit());
1367 float tmpSgnPhiWminusf = n3_Wb.Dot(nW.Cross(nS));
1368 float sgnPhiWminusf = 0;
1369 if (fabs(tmpSgnPhiWminusf)>0.) sgnPhiWminusf = tmpSgnPhiWminusf/fabs(tmpSgnPhiWminusf);
1370 float dot_nWnS = nW.Dot(nS);
1371 if (fabs(dot_nWnS)>=1.) dot_nWnS *= 1./fabs(dot_nWnS);
1372 PhiWminusf = sgnPhiWminusf * acos(dot_nWnS);
1388 const double GeV=1./100.;
1389 double ext_mZ_jhu = ext_mZ*GeV;
1390 double ext_mW_jhu = ext_mW*GeV;
1391 double ext_Gf_jhu = ext_Gf/pow(GeV, 2);
1395 if (ext_ewscheme<-1 || ext_ewscheme>3) ext_ewscheme=3;
1446 const int ipartabs = abs(ipart);
1447 bool runcoupling_mcfm=
false;
1448 bool runcoupling_jhugen=
false;
1454 else if (ipartabs==6){
masses_mcfm_.mt=inmass; runcoupling_mcfm=
true; }
1463 else if (ipartabs==23){
masses_mcfm_.zmass=inmass;
ewinput_.zmass_inp = inmass; runcoupling_mcfm=
true; }
1464 else if (ipartabs==24){
masses_mcfm_.wmass=inmass;
ewinput_.wmass_inp = inmass; runcoupling_mcfm=
true; }
1471 (ipartabs>=11 && ipartabs<=16)
1473 ipartabs==23 || ipartabs==24 || ipartabs==25
1475 ipartabs==32 || ipartabs==34
1477 runcoupling_jhugen=(
1482 const double GeV=1./100.;
1483 double jinmass = inmass*GeV;
1492 const int ipartabs = abs(ipart);
1502 const double GeV=1./100.;
1503 double jinwidth = inwidth*GeV;
1507 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){
1530 const int ipartabs = abs(ipart);
1551 (ipartabs>=11 && ipartabs<=16)
1553 ipartabs==23 || ipartabs==24 || ipartabs==25
1555 ipartabs==32 || ipartabs==34
1557 const double GeV=1./100.;
1560 double outmass = joutmass/GeV;
1567 const int ipartabs = abs(ipart);
1577 const double GeV=1./100.;
1580 double outwidth = joutwidth/GeV;
1585 if (part==
nullptr)
return GetMass(-9000);
1603 TLorentzVector
const nullFourVector(0, 0, 0, 0);
1612 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5];
1613 Q = fabs(pTotal.M());
1616 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5]+
p[6]+
p[7];
1617 Q = fabs(pTotal.M());
1620 TLorentzVector qH =
p[2]+
p[3]+
p[4]+
p[5];
1621 TLorentzVector qJJ =
p[6]+
p[7];
1622 Q = fabs(qH.M()+qJJ.M());
1625 TLorentzVector qH =
p[2]+
p[3]+
p[4]+
p[5];
1626 Q = fabs(qH.M()+
p[6].M()+
p[7].M());
1629 for (
int c=2; c<
mxpart; c++)
Q +=
p[c].Pt();
1633 for (
int c=6; c<
mxpart; c++)
Q = std::max(
Q,
p[c].Pt());
1638 for (
int c=7; c<
mxpart; c++){
if (
p[c].Pt()>0.)
Q = std::min(
Q,
p[c].Pt()); }
1658 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5];
1659 Q = fabs(pTotal.M());
1701 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5];
1702 Q = fabs(pTotal.M());
1708 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5]+
p[6]+
p[7];
1709 Q = fabs(pTotal.M());
1715 MELAerr <<
"Scaling " << scheme <<
" fails for production " <<
production <<
", defaulting to dynamic scheme m3456 " << endl;
1716 TLorentzVector pTotal =
p[2]+
p[3]+
p[4]+
p[5];
1717 Q = fabs(pTotal.M());
1721 void TUtil::SetAlphaS(
double& Q_ren,
double& Q_fac,
double multiplier_ren,
double multiplier_fac,
int mynloop,
int mynflav,
string mypartons){
1722 bool hasReset=
false;
1723 if (multiplier_ren<=0. || multiplier_fac<=0.){
1724 MELAerr <<
"TUtil::SetAlphaS: Invalid scale multipliers" << endl;
1727 if (Q_ren<=1. || Q_fac<=1. || mynloop<=0 || mypartons.compare(
"Default")==0){
1728 if (Q_ren<0.)
MELAout <<
"TUtil::SetAlphaS: Invalid QCD scale for alpha_s, setting to mH/2..." << endl;
1729 if (Q_fac<0.)
MELAout <<
"TUtil::SetAlphaS: Invalid factorization scale, setting to mH/2..." << endl;
1736 Q_ren *= multiplier_ren;
1737 Q_fac *= multiplier_fac;
1741 bool nflav_is_same = (
nflav_.nflav == mynflav);
1742 if (!nflav_is_same)
MELAout <<
"TUtil::SetAlphaS: nflav=" <<
nflav_.nflav <<
" is the only one supported." << endl;
1744 scale_.musq = Q_ren*Q_ren;
1747 MELAout <<
"TUtil::SetAlphaS: Only nloop=1 is supported!" << endl;
1775 const double GeV=1./100.;
1776 double muren_jhu =
scale_.scale*GeV;
1777 double mufac_jhu =
facscale_.facscale*GeV;
1795 double alphasVal=0, alphasmzVal=0;
1798 if (alphasmz_!=0) *alphasmz_ = alphasmzVal;
1809 unsigned int ndau = mela_event.
pDaughters.size();
1810 int* pId =
new int[ndau];
1811 for (
unsigned int ip=0; ip<ndau; ip++){
1821 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::MCFM_chooser: isWW=" << (
int)isWW <<
", isZZ=" << (
int)isZZ <<
", hasZZ4fInterf=" << (
int)hasZZ4fInterf << endl;
1883 (ndau>=4 || ndau==2)
2210 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;
2215 MELAerr <<
"TUtil::MCFM_chooser: Can't identify (process, production) = (" <<
process <<
", " <<
production <<
")" << endl;
2216 MELAerr <<
"TUtil::MCFM_chooser: ndau: " << ndau <<
'\t';
2217 MELAerr <<
"TUtil::MCFM_chooser: isZZ: " << isZZ <<
'\t';
2218 MELAerr <<
"TUtil::MCFM_chooser: isWW: " << isWW <<
'\t';
2219 MELAerr <<
"TUtil::MCFM_chooser: isGG: " << isGG <<
'\t';
2231 std::vector<int>* partOrder, std::vector<int>* apartOrder
2244 TString strplabel[
mxpart];
2245 for (
int ip=0; ip<
mxpart; ip++) strplabel[ip]=
" ";
2248 unsigned int ndau = mela_event.
pDaughters.size();
2249 if (ndau<1)
return false;
2250 unsigned int napart = mela_event.
pAssociated.size();
2265 bool hasZ1 = (isZZ || isZG || isZJJ);
2271 vector<TString> strcfgs;
2272 if (isGG) strcfgs.push_back(TString(
"GG"));
2273 if (isZG) strcfgs.push_back(TString(
"ZG"));
2274 if (isZZ) strcfgs.push_back(TString(
"ZZ"));
2275 if (isZJJ) strcfgs.push_back(TString(
"ZJJ"));
2276 if (isWW) strcfgs.push_back(TString(
"WW"));
2277 if (!strcfgs.empty()){
2279 for (
unsigned int istr=0; istr<strcfgs.size(); istr++){
2280 if (istr==0)
MELAout << strcfgs.at(istr);
2281 else MELAout <<
", " << strcfgs.at(istr);
2285 else MELAout <<
"has no valid decay!";
2290 int* pApartOrder = 0;
2293 pApartOrder =
new int[napart];
2294 pApartId =
new int[napart];
2295 for (
unsigned int ip=0; ip<napart; ip++){
2300 int* pOrder =
new int[ndau];
2301 int* pZOrder =
new int[ndau];
2302 int* pWOrder =
new int[ndau];
2303 int* pId =
new int[ndau];
2304 for (
unsigned int ip=0; ip<ndau; ip++){
2329 bool useQQVVQQstrong =
2335 bool useQQVVQQboth =
2341 bool useQQVVQQany = useQQVVQQ || useQQVVQQstrong || useQQVVQQboth;
2348 !(isWW || isZZ || isZJJ || isZG || isGG)
2354 if (isWW && !hasZ1 && !hasZ2){
2358 swap(pZOrder[0], pZOrder[2]);
2367 if (!hasZ1 && hasZ2){
2368 swap(pZOrder[0], pZOrder[2]);
2369 swap(pZOrder[1], pZOrder[3]);
2373 if (hasZ1)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Found a Z1";
2374 if (hasZ2)
MELAout <<
" and a Z2";
2375 if (hasZ1 || hasZ2)
MELAout <<
" in WW." << endl;
2378 if (isZZ && !hasW1 && !hasW2){
2379 swap(pWOrder[0], pWOrder[2]);
2387 if (V1id<0 || V2id>0){
2388 swap(pWOrder[0], pWOrder[2]);
2389 swap(pWOrder[1], pWOrder[3]);
2395 if (hasW1)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Found a W1(" << V1id <<
")";
2396 if (hasW2)
MELAout <<
" and a W2(" << V2id <<
")";
2397 if (hasW1 || hasW2)
MELAout <<
" in ZZ." << endl;
2409 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Setting up mother labels for MCFM:";
2410 for (
int ip=0; ip<min(2, (
int)mela_event.
pMothers.size()); ip++){
2411 const int& idmot = mela_event.
pMothers.at(ip).first;
2442 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2443 sprintf(
runstring_.runstring, strrun.c_str());
2470 else result =
false;
2482 else result =
false;
2494 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2495 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2496 sprintf(
runstring_.runstring, strrun.c_str());
2501 for (
unsigned int ix=0; ix<napart; ix++){
2503 for (
unsigned int iy=ix+1; iy<napart; iy++){
2507 int i1=ix;
int i2=iy;
2508 if (pApartId[ix]<0){ swap(i1, i2); }
2509 int* tmpOrder =
new int[napart];
2510 tmpOrder[0] = i1; tmpOrder[1] = i2;
2512 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2513 delete[] pApartOrder; pApartOrder = tmpOrder;
2518 int i1=ix;
int i2=iy;
2519 if (pApartId[ix]<0){ swap(i1, i2); }
2520 int* tmpOrder =
new int[napart];
2521 tmpOrder[0] = i1; tmpOrder[1] = i2;
2523 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2524 delete[] pApartOrder; pApartOrder = tmpOrder;
2528 if (hasZll || hasZnn)
break;
2531 unsigned int iout=6;
2532 unsigned int jout=7;
2534 int qqvvqq_apartordering[2]={ -1, -1 };
2537 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2538 qqvvqq_apartordering
2540 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2541 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]){
2543 swap(pApartOrder[0], pApartOrder[1]);
2547 strplabel[iout]=
"el";
2548 strplabel[jout]=
"ea";
2551 strplabel[iout]=
"nl";
2552 strplabel[jout]=
"na";
2554 else result =
false;
2562 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2563 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2564 sprintf(
runstring_.runstring, strrun.c_str());
2567 bool hasWplus=
false;
2568 bool hasWminus=
false;
2569 for (
unsigned int ix=0; ix<napart; ix++){
2571 for (
unsigned int iy=ix+1; iy<napart; iy++){
2583 int i1=ix;
int i2=iy;
2584 if (pApartId[ix]<0){ swap(i1, i2); }
2585 int* tmpOrder =
new int[napart];
2586 tmpOrder[0] = i1; tmpOrder[1] = i2;
2588 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2589 delete[] pApartOrder; pApartOrder = tmpOrder;
2602 int i1=ix;
int i2=iy;
2603 if (pApartId[ix]<0){ swap(i1, i2); }
2604 int* tmpOrder =
new int[napart];
2605 tmpOrder[0] = i1; tmpOrder[1] = i2;
2607 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2608 delete[] pApartOrder; pApartOrder = tmpOrder;
2613 if (hasWplus || hasWminus)
break;
2615 unsigned int iout=6;
2616 unsigned int jout=7;
2618 int qqvvqq_apartordering[2]={ -1, -1 };
2621 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2622 qqvvqq_apartordering
2624 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2625 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]){
2627 swap(pApartOrder[0], pApartOrder[1]);
2631 strplabel[iout]=
"nl";
2632 strplabel[jout]=
"ea";
2634 else if (hasWminus){
2635 strplabel[iout]=
"el";
2636 strplabel[jout]=
"na";
2638 else result =
false;
2646 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2647 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2648 sprintf(
runstring_.runstring, strrun.c_str());
2654 for (
unsigned int ix=0; ix<napart; ix++){
2656 for (
unsigned int iy=ix+1; iy<napart; iy++){
2662 int i1=ix;
int i2=iy;
2663 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2664 int* tmpOrder =
new int[napart];
2665 tmpOrder[0] = i1; tmpOrder[1] = i2;
2667 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2668 delete[] pApartOrder; pApartOrder = tmpOrder;
2673 int i1=ix;
int i2=iy;
2674 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2675 int* tmpOrder =
new int[napart];
2676 tmpOrder[0] = i1; tmpOrder[1] = i2;
2678 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2679 delete[] pApartOrder; pApartOrder = tmpOrder;
2684 int i1=ix;
int i2=iy;
2685 int* tmpOrder =
new int[napart];
2686 tmpOrder[0] = i1; tmpOrder[1] = i2;
2688 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2689 delete[] pApartOrder; pApartOrder = tmpOrder;
2693 if (hasZuu || hasZdd || hasZjj)
break;
2697 int qqvvqq_apartordering[2]={ -1, -1 };
2700 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2701 qqvvqq_apartordering
2703 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2704 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2706 if (hasZuu || hasZdd){
2714 else result =
false;
2722 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2723 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2724 sprintf(
runstring_.runstring, strrun.c_str());
2727 bool hasWplus=
false;
2728 bool hasWminus=
false;
2730 for (
unsigned int ix=0; ix<napart; ix++){
2732 for (
unsigned int iy=ix+1; iy<napart; iy++){
2750 int i1=ix;
int i2=iy;
2751 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2752 int* tmpOrder =
new int[napart];
2753 tmpOrder[0] = i1; tmpOrder[1] = i2;
2755 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2756 delete[] pApartOrder; pApartOrder = tmpOrder;
2769 int i1=ix;
int i2=iy;
2770 if (pApartId[ix]<0 || pApartId[iy]>0){ swap(i1, i2); }
2771 int* tmpOrder =
new int[napart];
2772 tmpOrder[0] = i1; tmpOrder[1] = i2;
2774 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2775 delete[] pApartOrder; pApartOrder = tmpOrder;
2780 int i1=ix;
int i2=iy;
2781 int* tmpOrder =
new int[napart];
2782 tmpOrder[0] = i1; tmpOrder[1] = i2;
2784 for (
unsigned int ipart=0; ipart<napart; ipart++){
if (ctr<napart && pApartOrder[ipart]!=i1 && pApartOrder[ipart]!=i2){ tmpOrder[ctr] = pApartOrder[ipart]; ctr++; } }
2785 delete[] pApartOrder; pApartOrder = tmpOrder;
2790 if (hasWplus || hasWminus || hasWjj)
break;
2793 int qqvvqq_apartordering[2]={ -1, -1 };
2796 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2797 qqvvqq_apartordering
2799 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2800 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2802 if (hasWplus || hasWminus){
2810 else result =
false;
2818 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2819 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2820 sprintf(
runstring_.runstring, strrun.c_str());
2824 int qqvvqq_apartordering[2]={ -1, -1 };
2827 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2828 qqvvqq_apartordering
2830 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2831 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2842 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2843 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2844 sprintf(
runstring_.runstring, strrun.c_str());
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]);
2866 if (isWW && (!hasZ1 || !hasZ2)) strrun +=
"_ww";
2867 else if (isZZ && (!hasW1 || !hasW2)) strrun +=
"_zz";
2868 sprintf(
runstring_.runstring, strrun.c_str());
2872 int qqvvqq_apartordering[2]={ -1, -1 };
2875 pApartId[pApartOrder[0]], pApartId[pApartOrder[1]],
2876 qqvvqq_apartordering
2878 if (qqvvqq_apartordering[0]==-1 || qqvvqq_apartordering[1]==-1) result=
false;
2879 else if (qqvvqq_apartordering[0]>qqvvqq_apartordering[1]) swap(pApartOrder[0], pApartOrder[1]);
2900 else if (useQQVVQQany){
2922 else if (useQQVVQQany){
2938 if (jetid==6) jetid = 2;
2948 else if (useQQVVQQany){
2962 double yy_up = pow(gV_up, 2) + pow(gA_up, 2);
2963 double yy_dn = pow(gV_dn, 2) + pow(gA_dn, 2);
2964 double xx_up = gV_up*gA_up;
2965 double xx_dn = gV_dn*gA_dn;
2966 double yy = (2.*yy_up+3.*yy_dn);
2967 double xx = (2.*xx_up+3.*xx_dn);
2968 double discriminant = pow(yy, 2)-4.*pow(xx, 2);
2969 double gVsq = (yy+sqrt(fabs(discriminant)))/2.;
2970 double gAsq = pow(xx, 2)/gVsq;
2971 double gV=-sqrt(gVsq);
2972 double gA=-sqrt(gAsq);
2982 else if (useQQVVQQany){
2998 else if (useQQVVQQany){
3054 if (jetid==6) jetid = 2;
3074 double yy_up = pow(gV_up, 2) + pow(gA_up, 2);
3075 double yy_dn = pow(gV_dn, 2) + pow(gA_dn, 2);
3076 double xx_up = gV_up*gA_up;
3077 double xx_dn = gV_dn*gA_dn;
3078 double yy = (2.*yy_up+3.*yy_dn);
3079 double xx = (2.*xx_up+3.*xx_dn);
3080 double discriminant = pow(yy, 2)-4.*pow(xx, 2);
3081 double gVsq = (yy+sqrt(fabs(discriminant)))/2.;
3082 double gAsq = pow(xx, 2)/gVsq;
3083 double gV=-sqrt(gVsq);
3084 double gA=-sqrt(gAsq);
3106 else if (useQQVVQQany){
3117 int* ordering=
nullptr;
3145 else if (
process ==
TVar::bkgZZ){ ordering = pZOrder;
if (!hasZ1 || !hasZ2) result =
false; }
3146 else if (
process ==
TVar::bkgWW){ ordering = pWOrder;
if (!hasW1 || !hasW2) result =
false; }
3149 else{ ordering = pZOrder;
if (!hasZ1 || !hasZ2) result =
false; }
3166 if (!hasZ1 || !hasZ2 || !hasW1 || !hasW2) result =
false;
3170 if (!hasW1 || !hasW2) result =
false;
3202 if (!hasZ1 || !hasZ2 || !hasW1 || !hasW2) result =
false;
3206 if (!hasW1 || !hasW2) result =
false;
3212 if (partOrder && ordering) {
for (
unsigned int ip=0; ip<ndau; ip++) partOrder->push_back(ordering[ip]); }
3213 if (apartOrder && pApartOrder) {
for (
unsigned int ip=0; ip<napart; ip++) apartOrder->push_back(pApartOrder[ip]); }
3214 for (
int ip=0; ip<
mxpart; ip++) sprintf((
plabel_.plabel)[ip], strplabel[ip].Data());
3217 MELAout <<
"TUtil::MCFM_SetupParticleCouplings: Summary (result=" << (
int)result <<
"):\n";
3221 if (hasZ1)
MELAout <<
"\tProcess found a Z1." << endl;
3222 if (hasZ2)
MELAout <<
"\tProcess found a Z2." << endl;
3223 if (hasW1)
MELAout <<
"\tProcess found a W1." << endl;
3224 if (hasW2)
MELAout <<
"\tProcess found a W2." << endl;
3230 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
3232 if (!partOrder->empty())
MELAout <<
"\tpartOrder[" << partOrder->size() <<
"]:\n";
3233 for (
unsigned int ip=0; ip<partOrder->size(); ip++)
MELAout <<
"\t[" << ip <<
"] -> " << partOrder->at(ip) << endl;
3236 if (!apartOrder->empty())
MELAout <<
"\tapartOrder[" << apartOrder->size() <<
"]:\n";
3237 for (
unsigned int ip=0; ip<apartOrder->size(); ip++)
MELAout <<
"\t[" << ip <<
"] -> " << apartOrder->at(ip) << endl;
3245 if (pApartOrder)
delete[] pApartOrder;
3246 if (pApartId)
delete[] pApartId;
3251 if (useQJ)
return TString(
"qj");
3252 else return TString(
"pp");
3255 else if (pid==1)
return TString(
"dq");
3256 else if (pid==2)
return TString(
"uq");
3257 else if (pid==3)
return TString(
"sq");
3258 else if (pid==4)
return TString(
"cq");
3259 else if (pid==5)
return TString(
"bq");
3260 else if (pid==-1)
return TString(
"da");
3261 else if (pid==-2)
return TString(
"ua");
3262 else if (pid==-3)
return TString(
"sa");
3263 else if (pid==-4)
return TString(
"ca");
3264 else if (pid==-5)
return TString(
"ba");
3265 else if (abs(pid)==6)
return TString(
"qj");
3266 else if (pid==11)
return TString(
"el");
3267 else if (pid==13)
return TString(
"ml");
3268 else if (pid==15)
return TString(
"tl");
3269 else if (pid==-11)
return TString(
"ea");
3270 else if (pid==-13)
return TString(
"ma");
3271 else if (pid==-15)
return TString(
"ta");
3272 else if (std::abs(pid)>=12 && std::abs(pid)<=16){
3273 if (!useExtendedConventions){
3274 if (pid>0)
return TString(
"nl");
3275 else return TString(
"na");
3278 if (pid==12)
return TString(
"ne");
3279 else if (pid==14)
return TString(
"nm");
3280 else if (pid==16)
return TString(
"nt");
3281 else if (pid==-12)
return TString(
"ke");
3282 else if (pid==-14)
return TString(
"km");
3283 else return TString(
"kt");
3286 else return TString(
" ");
3290 double minZmassSqr=10*10;
3294 (s[2][3]<minZmassSqr || s[4][5]<minZmassSqr)
3306 (-s[5-1][1-1]<
cutoff_.cutoff)
3307 || (-s[5-1][2-1]<
cutoff_.cutoff)
3308 || (-s[4-1][1-1]<
cutoff_.cutoff)
3309 || (-s[4-1][2-1]<
cutoff_.cutoff)
3310 || (-s[3-1][1-1]<
cutoff_.cutoff)
3311 || (-s[3-1][2-1]<
cutoff_.cutoff)
3312 || (+s[5-1][4-1]<
cutoff_.cutoff)
3313 || (+s[5-1][3-1]<
cutoff_.cutoff)
3314 || (+s[4-1][3-1]<
cutoff_.cutoff)
3322 (-s[5-1][1-1]<
cutoff_.cutoff)
3323 || (-s[5-1][2-1]<
cutoff_.cutoff)
3324 || (-s[6-1][1-1]<
cutoff_.cutoff)
3325 || (-s[6-1][2-1]<
cutoff_.cutoff)
3326 || (+s[6-1][5-1]<
cutoff_.cutoff)
3338 const double GeV = 1./100.;
3339 collider_sqrts *= GeV;
3341 char path_pdf_c[200];
3342 sprintf(path_pdf_c,
"%s", pathtoPDFSet);
3343 int pathpdfLength = strlen(path_pdf_c);
3348 const double GeV = 1./100.;
3354 int iAllow = (doAllow ? 1 : 0);
3404 for (
int im=0; im<2; im++){
3532 for (
int im=0; im<2; im++){
3642 for (
int im=0; im<2; im++){
3718 for (
int im=0; im<2; im++){
3771 for (
int im=0; im<2; im++){
3826 for (
int im=0; im<2; im++){
3902 for (
int im=0; im<2; im++){
3955 for (
int im=0; im<2; im++){
3995 for (
int im=0; im<2; im++){
4027 for (
int im=0; im<2; im++){
4079 const double GeV = 1./100.;
4080 int iWWcoupl = (useWWcoupl ? 1 : 0);
4106 const double& EBEAM,
4112 int nRequested_AssociatedJets=0;
4113 int nRequested_AssociatedLeptons=0;
4114 int AssociationVCompatibility=0;
4118 nRequested_AssociatedJets = 1;
4119 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::SumMatrixElementPDF: Requesting " << nRequested_AssociatedJets <<
" jets" << endl;
4132 nRequested_AssociatedJets = 2;
4133 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::SumMatrixElementPDF: Requesting " << nRequested_AssociatedJets <<
" jets" << endl;
4141 nRequested_AssociatedLeptons = 2;
4147 ) AssociationVCompatibility=24;
4152 ) AssociationVCompatibility=23;
4165 vector<int> partOrder;
4166 vector<int> apartOrder;
4173 if (partOrder.size()!=mela_event.
pDaughters.size()){
4175 MELAerr <<
"TUtil::SumMatrixElementPDF: Ordering size " << partOrder.size() <<
" and number of daughter particles " << mela_event.
pDaughters.size() <<
" are not the same!" << endl;
4180 if (apartOrder.size()!=mela_event.
pAssociated.size()){
4182 MELAerr <<
"TUtil::SumMatrixElementPDF: Ordering size " << apartOrder.size() <<
" and number of associated particles " << mela_event.
pAssociated.size() <<
" are not the same!" << endl;
4189 int NPart=
npart_.npart+2;
4190 double p4[4][
mxpart]={ { 0 } };
4191 double p4_tmp[4][
mxpart]={ { 0 } };
4192 int id[
mxpart];
for (
int ipar=0; ipar<
mxpart; ipar++)
id[ipar]=-9000;
4194 double msq_tmp[
nmsq][
nmsq]={ { 0 } };
4195 int channeltoggle=0;
4197 TLorentzVector MomStore[
mxpart];
4198 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
4202 for (
int ipar=0; ipar<2; ipar++){
4203 if (mela_event.
pMothers.at(ipar).second.T()>0.){
4204 p4[0][ipar] = -mela_event.
pMothers.at(ipar).second.X();
4205 p4[1][ipar] = -mela_event.
pMothers.at(ipar).second.Y();
4206 p4[2][ipar] = -mela_event.
pMothers.at(ipar).second.Z();
4207 p4[3][ipar] = -mela_event.
pMothers.at(ipar).second.T();
4208 MomStore[ipar] = mela_event.
pMothers.at(ipar).second;
4209 id[ipar] = mela_event.
pMothers.at(ipar).first;
4212 p4[0][ipar] = mela_event.
pMothers.at(ipar).second.X();
4213 p4[1][ipar] = mela_event.
pMothers.at(ipar).second.Y();
4214 p4[2][ipar] = mela_event.
pMothers.at(ipar).second.Z();
4215 p4[3][ipar] = mela_event.
pMothers.at(ipar).second.T();
4216 MomStore[ipar] = -mela_event.
pMothers.at(ipar).second;
4217 id[ipar] = mela_event.
pMothers.at(ipar).first;
4227 for (
int ix=0; ix<(
int)partOrder.size(); ix++){
4228 int ipar = min((
int)
mxpart, min(NPart, 2))+ix;
4230 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(partOrder.at(ix)).second);
4231 p4[0][ipar] = momTmp->X();
4232 p4[1][ipar] = momTmp->Y();
4233 p4[2][ipar] = momTmp->Z();
4234 p4[3][ipar] = momTmp->T();
4235 MomStore[ipar]=*momTmp;
4236 id[ipar] = mela_event.
pDaughters.at(partOrder.at(ix)).first;
4238 for (
int ix=0; ix<(
int)apartOrder.size(); ix++){
4239 int ipar = min((
int)
mxpart, min(NPart, (
int)(partOrder.size()+2)))+ix;
4241 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(apartOrder.at(ix)).second);
4242 p4[0][ipar] = momTmp->X();
4243 p4[1][ipar] = momTmp->Y();
4244 p4[2][ipar] = momTmp->Z();
4245 p4[3][ipar] = momTmp->T();
4246 MomStore[ipar]=*momTmp;
4247 id[ipar] = mela_event.
pAssociated.at(apartOrder.at(ix)).first;
4250 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; }
4252 double defaultRenScale =
scale_.scale;
4253 double defaultFacScale =
facscale_.facscale;
4255 int defaultNflav =
nflav_.nflav;
4256 string defaultPdflabel =
pdlabel_.pdlabel;
4260 double alphasVal, alphasmzVal;
4270 <<
"TUtil::SumMatrixElementPDF: Set AlphaS:\n"
4271 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
4273 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
4291 else channeltoggle=2;
4298 for (
int iquark=-5; iquark<=5; iquark++){
4299 for (
int jquark=-5; jquark<=5; jquark++){
4309 ) msqjk =
msq[3][7]+
msq[7][3];
4312 ) msqjk =
msq[jquark+5][-jquark+5];
4315 ) msqjk =
msq[-iquark+5][iquark+5];
4316 else msqjk +=
msq[jquark+5][iquark+5];
4320 SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4333 MELAout <<
"\tTUtil::SumMatrixElementPDF: ZZGG && ZZ/WW/ZZ+WW MEs using ZZ (runstring: " <<
runstring_.runstring <<
")" << endl;
4334 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;
4342 MELAout <<
"\tTUtil::SumMatrixElementPDF: ZZGG && ZZ/WW/ZZ+WW MEs using WW (runstring: " <<
runstring_.runstring <<
")" << endl;
4343 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;
4349 for (
int iquark=-5; iquark<=5; iquark++){
4350 for (
int jquark=-5; jquark<=5; jquark++){
4351 if (iquark==0 && jquark==0)
continue;
4352 msq[jquark+5][iquark+5] = 0;
4366 SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4411 for (
unsigned int ix=0; ix<4; ix++){
4412 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4413 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4416 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; } }
4419 for (
unsigned int ix=0; ix<4; ix++){
4420 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4421 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4423 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4424 for (
int ip=0; ip<
mxpart; ip++){
4425 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4426 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4430 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4432 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4433 MELAout <<
"\tMEsq initial:" << endl;
4434 for (
int iquark=-5; iquark<=5; iquark++){
4435 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4438 MELAout <<
"\tMEsq added:" << endl;
4439 for (
int iquark=-5; iquark<=5; iquark++){
4440 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4444 for (
int iquark=-5; iquark<=5; iquark++){
4445 for (
int jquark=-5; jquark<=5; jquark++){
4446 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4457 for (
unsigned int ix=0; ix<4; ix++){
4458 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4459 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4462 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; } }
4465 for (
int ip=0; ip<
mxpart; ip++){
4466 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4467 else sprintf((
plabel_.plabel)[ip], gglabel.Data());
4470 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]); }
4472 MELAout <<
"TUtil::SumMatrixElementPDF: Subtracting qqb/qbq->gg double-counted contribution:\n";
4474 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4476 for (
int iquark=-5; iquark<=5; iquark++){
4477 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4483 for (
unsigned int ix=0; ix<4; ix++){
4484 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4485 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4487 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4488 for (
int ip=0; ip<
mxpart; ip++){
4489 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4490 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4494 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4496 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4497 MELAout <<
"\tMEsq initial:" << endl;
4498 for (
int iquark=-5; iquark<=5; iquark++){
4499 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4502 MELAout <<
"\tMEsq added:" << endl;
4503 for (
int iquark=-5; iquark<=5; iquark++){
4504 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4508 for (
int iquark=-5; iquark<=5; iquark++){
4509 for (
int jquark=-5; jquark<=5; jquark++){
4510 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4530 for (
unsigned int ix=0; ix<4; ix++){
4531 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4532 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4535 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; } }
4538 for (
unsigned int ix=0; ix<4; ix++){
4539 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4540 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4542 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4543 for (
int ip=0; ip<
mxpart; ip++){
4544 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4545 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4549 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4551 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4552 MELAout <<
"\tMEsq initial:" << endl;
4553 for (
int iquark=-5; iquark<=5; iquark++){
4554 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4557 MELAout <<
"\tMEsq added:" << endl;
4558 for (
int iquark=-5; iquark<=5; iquark++){
4559 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4563 for (
int iquark=-5; iquark<=5; iquark++){
4564 for (
int jquark=-5; jquark<=5; jquark++){
4565 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4575 for (
unsigned int ix=0; ix<4; ix++){
4576 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4577 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4580 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; } }
4583 for (
int ip=0; ip<
mxpart; ip++){
4584 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4585 else sprintf((
plabel_.plabel)[ip], gglabel.Data());
4588 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]); }
4590 MELAout <<
"TUtil::SumMatrixElementPDF: Subtracting qqb/qbq->gg double-counted contribution:\n";
4592 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4594 for (
int iquark=-5; iquark<=5; iquark++){
4595 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4601 for (
unsigned int ix=0; ix<4; ix++){
4602 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4603 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4605 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
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], strlabel[ip-6].Data());
4612 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4614 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4615 MELAout <<
"\tMEsq initial:" << endl;
4616 for (
int iquark=-5; iquark<=5; iquark++){
4617 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4620 MELAout <<
"\tMEsq added:" << endl;
4621 for (
int iquark=-5; iquark<=5; iquark++){
4622 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4626 for (
int iquark=-5; iquark<=5; iquark++){
4627 for (
int jquark=-5; jquark<=5; jquark++){
4628 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4648 for (
unsigned int ix=0; ix<4; ix++){
4649 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4650 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4653 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; } }
4656 for (
unsigned int ix=0; ix<4; ix++){
4657 for (
int ip=0; ip<
mxpart; ip++) p4_tmp[ix][ip]=p4[ix][ip];
4658 swap(p4_tmp[ix][partOrder.size()+2], p4_tmp[ix][partOrder.size()+3]);
4660 TString strlabel[2] ={ (
plabel_.plabel)[7], (
plabel_.plabel)[6] };
for (
unsigned int il=0; il<2; il++) strlabel[il].Resize(2);
4661 for (
int ip=0; ip<
mxpart; ip++){
4662 if (ip!=6 && ip!=7) sprintf((
plabel_.plabel)[ip], (
plabel_.plabel)[ip]);
4663 else sprintf((
plabel_.plabel)[ip], strlabel[ip-6].Data());
4667 MELAout <<
"TUtil::SumMatrixElementPDF: Adding missing contributions:\n";
4669 for (
int ip=0; ip<
mxpart; ip++)
MELAout <<
"\t[" << ip <<
"]=" << (
plabel_.plabel)[ip] << endl;
4670 MELAout <<
"\tMEsq initial:" << endl;
4671 for (
int iquark=-5; iquark<=5; iquark++){
4672 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4675 MELAout <<
"\tMEsq added:" << endl;
4676 for (
int iquark=-5; iquark<=5; iquark++){
4677 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << msq_tmp[jquark+5][iquark+5] <<
'\t';
4681 for (
int iquark=-5; iquark<=5; iquark++){
4682 for (
int jquark=-5; jquark<=5; jquark++){
4683 if (
msq[jquark+5][iquark+5]==0.)
msq[jquark+5][iquark+5] = msq_tmp[jquark+5][iquark+5];
4694 msqjk =
SumMEPDF(MomStore[0], MomStore[1],
msq, RcdME, EBEAM, verbosity);
4711 if (msqjk != msqjk){
4713 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;
4714 MELAout <<
"TUtil::SumMatrixElementPDF: The number of ME instances: " << nInstances <<
". MEsq[ip][jp] = " << endl;
4715 for (
int iquark=-5; iquark<=5; iquark++){
4716 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4722 MELAout <<
"TUtil::SumMatrixElementPDF: The number of ME instances: " << nInstances <<
". MEsq[ip][jp] = " << endl;
4723 for (
int iquark=-5; iquark<=5; iquark++){
4724 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
4730 <<
"TUtil::SumMatrixElementPDF: Reset AlphaS:\n"
4731 <<
"\tBefore reset, alphas scale: " <<
scale_.scale
4732 <<
", PDF scale: " <<
facscale_.facscale
4734 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
4736 <<
"TUtil::SumMatrixElementPDF: Reset AlphaS result:\n"
4737 <<
"\tAfter reset, alphas scale: " <<
scale_.scale
4738 <<
", PDF scale: " <<
facscale_.facscale
4741 if (verbosity>=
TVar::DEBUG)
MELAout <<
"End TUtil::SumMatrixElementPDF(" << msqjk <<
")" << endl;
4749 const double& EBEAM,
4752 const double GeV=1./100.;
4755 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::JHUGenMatEl: Non-JHUGen MEs are not supported" << endl;
return MatElSq; }
4790 if (isSpinZero)
MELAout <<
"0";
4791 else if (isSpinOne)
MELAout <<
"1";
4792 else if (isSpinTwo)
MELAout <<
"2";
4798 int MYIDUP_tmp[4]={ 0 };
4799 vector<pair<int, int>> idarray[2];
4800 int MYIDUP[4]={ 0 };
4801 int idfirst[2]={ 0 };
4802 int idsecond[2]={ 0 };
4803 double p4[6][4]={ { 0 } };
4804 TLorentzVector MomStore[
mxpart];
4805 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
4825 for (
int ipar=0; ipar<2; ipar++){
4826 if (mela_event.
pMothers.at(ipar).second.T()>0.){
4827 p4[ipar][0] = -mela_event.
pMothers.at(ipar).second.T()*GeV;
4828 p4[ipar][1] = -mela_event.
pMothers.at(ipar).second.X()*GeV;
4829 p4[ipar][2] = -mela_event.
pMothers.at(ipar).second.Y()*GeV;
4830 p4[ipar][3] = -mela_event.
pMothers.at(ipar).second.Z()*GeV;
4831 MomStore[ipar] = mela_event.
pMothers.at(ipar).second;
4834 p4[ipar][0] = mela_event.
pMothers.at(ipar).second.T()*GeV;
4835 p4[ipar][1] = mela_event.
pMothers.at(ipar).second.X()*GeV;
4836 p4[ipar][2] = mela_event.
pMothers.at(ipar).second.Y()*GeV;
4837 p4[ipar][3] = mela_event.
pMothers.at(ipar).second.Z()*GeV;
4838 MomStore[ipar] = -mela_event.
pMothers.at(ipar).second;
4846 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
4847 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
4848 int* idtmp = &(mela_event.
pDaughters.at(ipar).first);
4850 int arrindex = 2*ipar;
4852 else MYIDUP_tmp[arrindex] = 0;
4853 MYIDUP_tmp[arrindex+1] = -9000;
4854 p4[arrindex+2][0] = momTmp->T()*GeV;
4855 p4[arrindex+2][1] = momTmp->X()*GeV;
4856 p4[arrindex+2][2] = momTmp->Y()*GeV;
4857 p4[arrindex+2][3] = momTmp->Z()*GeV;
4858 MomStore[arrindex+2] = *momTmp;
4859 if (verbosity >=
TVar::DEBUG)
MELAout <<
"MYIDUP_tmp[" << arrindex <<
"(" << ipar <<
")" <<
"]=" << MYIDUP_tmp[arrindex] << endl;
4863 for (
unsigned int ipar=0; ipar<4; ipar++){
4865 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
4866 int* idtmp = &(mela_event.
pDaughters.at(ipar).first);
4869 else MYIDUP_tmp[ipar] = 0;
4870 p4[ipar+2][0] = momTmp->T()*GeV;
4871 p4[ipar+2][1] = momTmp->X()*GeV;
4872 p4[ipar+2][2] = momTmp->Y()*GeV;
4873 p4[ipar+2][3] = momTmp->Z()*GeV;
4874 MomStore[ipar+2] = *momTmp;
4876 else MYIDUP_tmp[ipar] = -9000;
4878 if (verbosity >=
TVar::DEBUG)
MELAout <<
"MYIDUP_tmp[" << ipar <<
"]=" << MYIDUP_tmp[ipar] << endl;
4883 double defaultRenScale =
scale_.scale;
4884 double defaultFacScale =
facscale_.facscale;
4886 int defaultNflav =
nflav_.nflav;
4887 string defaultPdflabel =
pdlabel_.pdlabel;
4891 double alphasVal, alphasmzVal;
4901 <<
"TUtil::JHUGenMatEl: Set AlphaS:\n"
4902 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
4904 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
4908 for (
int iv=0; iv<2; iv++){
4909 if (MYIDUP_tmp[2*iv+0]!=0 && MYIDUP_tmp[2*iv+1]!=0){
4912 TMath::Sign(1, MYIDUP_tmp[2*iv+0])==TMath::Sign(1, -MYIDUP_tmp[2*iv+1])
4914 (MYIDUP_tmp[2*iv+0]==-9000 || MYIDUP_tmp[2*iv+1]==-9000)
4915 ) idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], MYIDUP_tmp[2*iv+1]));
4917 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]));
4918 else idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], -MYIDUP_tmp[2*iv+1]));
4920 else if (MYIDUP_tmp[2*iv+0]!=0){
4923 int id_dn = TMath::Sign(1, -MYIDUP_tmp[2*iv+0]);
4924 int id_st = TMath::Sign(3, -MYIDUP_tmp[2*iv+0]);
4925 int id_bt = TMath::Sign(5, -MYIDUP_tmp[2*iv+0]);
4926 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_dn));
4927 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_st));
4928 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_bt));
4931 int id_up = TMath::Sign(2, -MYIDUP_tmp[2*iv+0]);
4932 int id_ch = TMath::Sign(4, -MYIDUP_tmp[2*iv+0]);
4933 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_up));
4934 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], id_ch));
4938 idarray[iv].push_back(pair<int, int>(MYIDUP_tmp[2*iv+0], -MYIDUP_tmp[2*iv+0]));
4941 else if (MYIDUP_tmp[2*iv+1]!=0){
4944 int id_dn = TMath::Sign(1, -MYIDUP_tmp[2*iv+1]);
4945 int id_st = TMath::Sign(3, -MYIDUP_tmp[2*iv+1]);
4946 int id_bt = TMath::Sign(5, -MYIDUP_tmp[2*iv+1]);
4947 idarray[iv].push_back(pair<int, int>(id_dn, MYIDUP_tmp[2*iv+1]));
4948 idarray[iv].push_back(pair<int, int>(id_st, MYIDUP_tmp[2*iv+1]));
4949 idarray[iv].push_back(pair<int, int>(id_bt, MYIDUP_tmp[2*iv+1]));
4952 int id_up = TMath::Sign(2, -MYIDUP_tmp[2*iv+1]);
4953 int id_ch = TMath::Sign(4, -MYIDUP_tmp[2*iv+1]);
4954 idarray[iv].push_back(pair<int, int>(id_up, MYIDUP_tmp[2*iv+1]));
4955 idarray[iv].push_back(pair<int, int>(id_ch, MYIDUP_tmp[2*iv+1]));
4959 idarray[iv].push_back(pair<int, int>(-MYIDUP_tmp[2*iv+1], MYIDUP_tmp[2*iv+1]));
4964 for (
int iqu=1; iqu<=2; iqu++){
4966 for (
int iqd=1; iqd<=3; iqd++){
4967 int id_dn = iqd*2-1;
4968 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)); }
4969 else{ idarray[iv].push_back(pair<int, int>(id_dn, -id_up)); idarray[iv].push_back(pair<int, int>(-id_up, id_dn)); }
4974 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)); }
4980 const int InvalidMode=-1, WWMode=0, ZZMode=1, ZgMode=5, ggMode=7;
4981 int VVMode=InvalidMode;
4990 double aL1=0, aL2=0, aR1=0, aR2=0;
4993 for (
unsigned int v1=0; v1<idarray[0].size(); v1++){
4994 for (
unsigned int v2=0; v2<idarray[1].size(); v2++){
4996 if (idarray[0].at(v1).first!=-9000) MYIDUP[0] =
convertLHEreverse(&(idarray[0].at(v1).first));
4998 if (idarray[1].at(v2).first!=-9000) MYIDUP[2] =
convertLHEreverse(&(idarray[1].at(v2).first));
5002 if (verbosity>=
TVar::DEBUG){
for (
unsigned int idau=0; idau<4; idau++)
MELAout <<
"MYIDUP[" << idau <<
"]=" << MYIDUP[idau] << endl; }
5005 for (
int ip=0; ip<2; ip++){ idfirst[ip]=MYIDUP[ip]; idsecond[ip]=MYIDUP[ip+2]; }
5010 MELAout <<
"TUtil::JHUGenMatEl: M_V=" << mv/GeV <<
", Ga_V=" << gv/GeV << endl;
5012 MELAout <<
"TUtil::JHUGenMatEl: M_Vprime=" << mv/GeV <<
", Ga_Vprime=" << gv/GeV << endl;
5016 double aLRtmp[4]={ 0 };
5018 if (idarray[0].size()>1){
5019 aL1 = sqrt(pow(aL1, 2)+pow(aLRtmp[0], 2));
5020 aR1 = sqrt(pow(aR1, 2)+pow(aLRtmp[1], 2));
5026 if (idarray[1].size()>1){
5027 aL2 = sqrt(pow(aL2, 2)+pow(aLRtmp[2], 2));
5028 aR2 = sqrt(pow(aR2, 2)+pow(aLRtmp[3], 2));
5053 if (verbosity >=
TVar::DEBUG)
MELAout <<
"=====\nTUtil::JHUGenMatEl: Instance MatElTmp = " << MatElTmp <<
"\n=====" << endl;
5054 MatElSq += MatElTmp;
5055 if (MatElTmp>0.) nNonZero++;
5058 if (nNonZero>0) MatElSq /= ((double)nNonZero);
5060 MELAout <<
"TUtil::JHUGenMatEl: Number of matrix element instances computed: " << nNonZero << endl;
5061 MELAout <<
"TUtil::JHUGenMatEl: MatElSq after division = " << MatElSq << endl;
5068 <<
"TUtil::JHUGenMatEl: aL1, aR1, aL2, aR2: "
5069 << aL1 <<
", " << aR1 <<
", " << aL2 <<
", " << aR2
5073 int GeVexponent_MEsq = 4-((
int)mela_event.
pDaughters.size())*2;
5075 double constant = pow(GeV, -GeVexponent_MEsq);
5076 MatElSq *= constant;
5079 if (mela_event.
pMothers.at(0).first!=0 && mela_event.
pMothers.at(1).first!=0){
5081 if (abs(mela_event.
pMothers.at(0).first)<6) ix=mela_event.
pMothers.at(0).first + 5;
5082 if (abs(mela_event.
pMothers.at(1).first)<6) iy=mela_event.
pMothers.at(1).first + 5;
5083 msq[iy][ix]=MatElSq;
5088 for (
int ix=0; ix<5; ix++){
msq[ix][10-ix]=MatElSq;
msq[10-ix][ix]=MatElSq; }
5093 double fx_dummy[
nmsq]={ 0 }; fx_dummy[5]=1.;
5099 if (verbosity >=
TVar::DEBUG)
MELAout <<
"TUtil::JHUGenMatEl: Final MatElSq = " << MatElSq << endl;
5104 <<
"TUtil::JHUGenMatEl: Reset AlphaS:\n"
5105 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
5107 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
5111 <<
"TUtil::JHUGenMatEl: Reset AlphaS result:\n"
5112 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5121 const double& EBEAM,
5124 const double GeV=1./100.;
5125 double sum_msqjk = 0;
5133 double MatElsq[
nmsq][
nmsq]={ { 0 } };
5134 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
5137 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::HJJMatEl: Non-JHUGen MEs are not supported" << endl;
return sum_msqjk; }
5141 int nRequested_AssociatedJets=2;
5152 if (mela_event.
pAssociated.size()==0){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::HJJMatEl: Number of associated particles is 0!" << endl;
return sum_msqjk; }
5154 int MYIDUP_tmp[4]={ 0 };
5155 double p4[5][4]={ { 0 } };
5156 double pOneJet[4][4] ={ { 0 } };
5157 TLorentzVector MomStore[
mxpart];
5158 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
5164 for (
int ipar=0; ipar<2; ipar++){
5165 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
5166 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
5168 else MYIDUP_tmp[ipar] = 0;
5169 if (momTmp->T()>0.){
5170 p4[ipar][0] = momTmp->T()*GeV;
5171 p4[ipar][1] = momTmp->X()*GeV;
5172 p4[ipar][2] = momTmp->Y()*GeV;
5173 p4[ipar][3] = momTmp->Z()*GeV;
5174 MomStore[ipar] = (*momTmp);
5177 p4[ipar][0] = -momTmp->T()*GeV;
5178 p4[ipar][1] = -momTmp->X()*GeV;
5179 p4[ipar][2] = -momTmp->Y()*GeV;
5180 p4[ipar][3] = -momTmp->Z()*GeV;
5181 MomStore[ipar] = -(*momTmp);
5182 MYIDUP_tmp[ipar] = -MYIDUP_tmp[ipar];
5185 for (
unsigned int ipar=0; ipar<2; ipar++){
5187 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(ipar).second);
5188 int* idtmp = &(mela_event.
pAssociated.at(ipar).first);
5190 else MYIDUP_tmp[ipar+2] = 0;
5191 p4[ipar+2][0] = momTmp->T()*GeV;
5192 p4[ipar+2][1] = momTmp->X()*GeV;
5193 p4[ipar+2][2] = momTmp->Y()*GeV;
5194 p4[ipar+2][3] = momTmp->Z()*GeV;
5195 MomStore[ipar+6] = (*momTmp);
5197 else MYIDUP_tmp[ipar+2] = -9000;
5199 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
5200 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
5201 p4[4][0] += momTmp->T()*GeV;
5202 p4[4][1] += momTmp->X()*GeV;
5203 p4[4][2] += momTmp->Y()*GeV;
5204 p4[4][3] += momTmp->Z()*GeV;
5205 MomStore[5] = MomStore[5] + (*momTmp);
5208 for (
unsigned int i = 0; i < 4; i++){
5209 if (i<2){
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[i][j]; }
5210 else if (i==2){
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[4][j]; }
5211 else{
for (
unsigned int j = 0; j < 4; j++) pOneJet[i][j] = p4[2][j]; }
5214 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;
5215 MELAout <<
"One-jet momenta: " << endl;
5216 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;
5219 double defaultRenScale =
scale_.scale;
5220 double defaultFacScale =
facscale_.facscale;
5222 int defaultNflav =
nflav_.nflav;
5223 string defaultPdflabel =
pdlabel_.pdlabel;
5227 double alphasVal, alphasmzVal;
5237 <<
"TUtil::HJJMatEl: Set AlphaS:\n"
5238 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
5240 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
5244 bool partonIsUnknown[4];
5245 for (
unsigned int ip=0; ip<4; ip++) partonIsUnknown[ip] = (MYIDUP_tmp[ip]==0);
5251 for (
int isel=-5; isel<=5; isel++){
5252 if (!partonIsUnknown[0] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
continue;
5253 for (
int jsel=-5; jsel<=5; jsel++){
5254 if (!partonIsUnknown[1] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
continue;
5256 if (isel!=0 && jsel!=0) rsel=0;
5257 else if (isel==0) rsel=jsel;
5259 if (!partonIsUnknown[2] && !((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
continue;
5260 MatElsq[jsel+5][isel+5] = MatElsq_tmp[jsel+5][isel+5];
5261 if (verbosity >=
TVar::DEBUG)
MELAout <<
"Channel (isel, jsel)=" << isel <<
", " << jsel << endl;
5267 int nijchannels = ijsel.size();
5268 for (
int ic=0; ic<nijchannels; ic++){
5270 int isel = ijsel[ic][0];
5271 int jsel = ijsel[ic][1];
5272 int code = ijsel[ic][2];
5282 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
5284 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
5287 if (isel==0 && jsel==0){
5298 MatElsq[jsel+5][isel+5] += msq_tmp;
5299 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5307 MatElsq[jsel+5][isel+5] += msq_tmp;
5308 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5319 MatElsq[jsel+5][isel+5] += msq_tmp;
5320 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5324 else if (isel==0 || jsel==0){
5326 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
5328 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && ssel==0) || MYIDUP_tmp[3]==ssel))
5331 MatElsq[jsel+5][isel+5] += msq_tmp;
5332 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5335 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && ssel==0) || MYIDUP_tmp[2]==ssel))
5337 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && rsel==0) || MYIDUP_tmp[3]==rsel))
5340 MatElsq[jsel+5][isel+5] += msq_tmp;
5341 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5344 else if ((isel>0 && jsel<0) || (isel<0 && jsel>0)){
5345 if (code==1 && isel==-jsel){
5353 MatElsq[jsel+5][isel+5] += msq_tmp;
5354 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5357 else if (code==3 && isel==-jsel){
5358 if (abs(isel)!=1){ rsel=1; ssel=-1; }
5359 else{ rsel=2; ssel=-2; }
5362 (partonIsUnknown[2] && partonIsUnknown[3])
5364 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]<0)
5366 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]>0 && partonIsUnknown[3])
5368 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]>0)
5371 MatElsq[jsel+5][isel+5] += msq_tmp;
5372 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5375 (partonIsUnknown[2] && partonIsUnknown[3])
5377 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]>0)
5379 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]<0 && partonIsUnknown[3])
5381 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]<0)
5384 MatElsq[jsel+5][isel+5] += msq_tmp;
5385 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5390 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5392 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5395 MatElsq[jsel+5][isel+5] += msq_tmp;
5396 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5399 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5401 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5404 MatElsq[jsel+5][isel+5] += msq_tmp;
5405 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5411 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5413 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5416 MatElsq[jsel+5][isel+5] += msq_tmp;
5417 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5422 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5424 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5427 MatElsq[jsel+5][isel+5] += msq_tmp;
5428 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5433 if (isel==jsel)
continue;
5434 isel = ijsel[ic][1];
5435 jsel = ijsel[ic][0];
5443 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
5445 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
5448 if (isel==0 || jsel==0){
5450 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && rsel==0) || MYIDUP_tmp[2]==rsel))
5452 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && ssel==0) || MYIDUP_tmp[3]==ssel))
5455 MatElsq[jsel+5][isel+5] += msq_tmp;
5456 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5459 (partonIsUnknown[2] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[2]) && ssel==0) || MYIDUP_tmp[2]==ssel))
5461 (partonIsUnknown[3] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[3]) && rsel==0) || MYIDUP_tmp[3]==rsel))
5464 MatElsq[jsel+5][isel+5] += msq_tmp;
5465 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5468 else if ((isel>0 && jsel<0) || (isel<0 && jsel>0)){
5469 if (code==1 && isel==-jsel){
5477 MatElsq[jsel+5][isel+5] += msq_tmp;
5478 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5481 else if(code==3 && isel==-jsel){
5482 if (abs(isel)!=1){ rsel=1; ssel=-1; }
5483 else{ rsel=2; ssel=-2; }
5486 (partonIsUnknown[2] && partonIsUnknown[3])
5488 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]<0)
5490 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]>0 && partonIsUnknown[3])
5492 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]>0)
5495 MatElsq[jsel+5][isel+5] += msq_tmp;
5496 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5499 (partonIsUnknown[2] && partonIsUnknown[3])
5501 (partonIsUnknown[2] && std::abs(MYIDUP_tmp[3])!=std::abs(isel) && MYIDUP_tmp[3]>0)
5503 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]<0 && partonIsUnknown[3])
5505 (std::abs(MYIDUP_tmp[2])!=std::abs(isel) && MYIDUP_tmp[2]==-MYIDUP_tmp[3] && MYIDUP_tmp[2]<0)
5508 MatElsq[jsel+5][isel+5] += msq_tmp;
5509 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5514 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5516 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5519 MatElsq[jsel+5][isel+5] += msq_tmp;
5520 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5523 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5525 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5528 MatElsq[jsel+5][isel+5] += msq_tmp;
5529 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5535 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
5537 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
5540 MatElsq[jsel+5][isel+5] += msq_tmp;
5541 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
5546 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
5548 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
5551 MatElsq[jsel+5][isel+5] += msq_tmp;
5552 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
5559 int isel, jsel, rsel, ssel;
5568 double ckm_takeout=1;
5569 double msq_uu_zz_ijrs1234[2]={ 0 };
5570 double msq_uu_zz_ijrs1243[2]={ 0 };
5571 double msq_dd_zz_ijrs1234[2]={ 0 };
5572 double msq_dd_zz_ijrs1243[2]={ 0 };
5573 double msq_ubarubar_zz_ijrs1234[2]={ 0 };
5574 double msq_ubarubar_zz_ijrs1243[2]={ 0 };
5575 double msq_dbardbar_zz_ijrs1234[2]={ 0 };
5576 double msq_dbardbar_zz_ijrs1243[2]={ 0 };
5577 double msq_uu_zzid_ijrs1234[2]={ 0 };
5578 double msq_uu_zzid_ijrs1243[2]={ 0 };
5579 double msq_dd_zzid_ijrs1234[2]={ 0 };
5580 double msq_dd_zzid_ijrs1243[2]={ 0 };
5581 double msq_ubarubar_zzid_ijrs1234[2]={ 0 };
5582 double msq_ubarubar_zzid_ijrs1243[2]={ 0 };
5583 double msq_dbardbar_zzid_ijrs1234[2]={ 0 };
5584 double msq_dbardbar_zzid_ijrs1243[2]={ 0 };
5585 double msq_udbar_zz_ijrs1234[2]={ 0 };
5586 double msq_udbar_zz_ijrs1243[2]={ 0 };
5587 double msq_dubar_zz_ijrs1234[2]={ 0 };
5588 double msq_dubar_zz_ijrs1243[2]={ 0 };
5590 double msq_uubar_zz_ijrs1234[2]={ 0 };
5591 double msq_uubar_zz_ijrs1243[2]={ 0 };
5592 double msq_ddbar_zz_ijrs1234[2]={ 0 };
5593 double msq_ddbar_zz_ijrs1243[2]={ 0 };
5594 double msq_uubar_ww_ijrs1234[2]={ 0 };
5595 double msq_uubar_ww_ijrs1243[2]={ 0 };
5596 double msq_ddbar_ww_ijrs1234[2]={ 0 };
5597 double msq_ddbar_ww_ijrs1243[2]={ 0 };
5599 double msq_ud_wwonly_ijrs1234[2]={ 0 };
5600 double msq_ud_wwonly_ijrs1243[2]={ 0 };
5601 double msq_ubardbar_wwonly_ijrs1234[2]={ 0 };
5602 double msq_ubardbar_wwonly_ijrs1243[2]={ 0 };
5616 rsel=isel; ssel=jsel;
5620 msq_uu_zz_ijrs1234[1] = msq_uu_zz_ijrs1234[0];
5621 msq_uu_zz_ijrs1243[1] = msq_uu_zz_ijrs1243[0];
5624 rsel=isel; ssel=jsel;
5628 msq_uu_zzid_ijrs1234[1] = msq_uu_zzid_ijrs1234[0];
5629 msq_uu_zzid_ijrs1243[1] = msq_uu_zzid_ijrs1243[0];
5642 rsel=isel; ssel=jsel;
5646 msq_dd_zz_ijrs1234[1] = msq_dd_zz_ijrs1234[0];
5647 msq_dd_zz_ijrs1243[1] = msq_dd_zz_ijrs1243[0];
5650 rsel=isel; ssel=jsel;
5654 msq_dd_zzid_ijrs1234[1] = msq_dd_zzid_ijrs1234[0];
5655 msq_dd_zzid_ijrs1243[1] = msq_dd_zzid_ijrs1243[0];
5668 rsel=isel; ssel=jsel;
5672 msq_ubarubar_zz_ijrs1234[1] = msq_ubarubar_zz_ijrs1234[0];
5673 msq_ubarubar_zz_ijrs1243[1] = msq_ubarubar_zz_ijrs1243[0];
5676 rsel=isel; ssel=jsel;
5680 msq_ubarubar_zzid_ijrs1234[1] = msq_ubarubar_zzid_ijrs1234[0];
5681 msq_ubarubar_zzid_ijrs1243[1] = msq_ubarubar_zzid_ijrs1243[0];
5694 rsel=isel; ssel=jsel;
5698 msq_dbardbar_zz_ijrs1234[1] = msq_dbardbar_zz_ijrs1234[0];
5699 msq_dbardbar_zz_ijrs1243[1] = msq_dbardbar_zz_ijrs1243[0];
5702 rsel=isel; ssel=jsel;
5706 msq_dbardbar_zzid_ijrs1234[1] = msq_dbardbar_zzid_ijrs1234[0];
5707 msq_dbardbar_zzid_ijrs1243[1] = msq_dbardbar_zzid_ijrs1243[0];
5712 (partonIsUnknown[0] && partonIsUnknown[1])
5726 (partonIsUnknown[2] && partonIsUnknown[3])
5740 rsel=isel; ssel=jsel;
5749 (partonIsUnknown[0] && partonIsUnknown[1])
5763 (partonIsUnknown[2] && partonIsUnknown[3])
5777 rsel=isel; ssel=jsel;
5785 (partonIsUnknown[0] && partonIsUnknown[1])
5795 (partonIsUnknown[2] && partonIsUnknown[3])
5803 rsel=isel; ssel=jsel;
5810 (partonIsUnknown[2] && partonIsUnknown[3])
5829 if (abs(isel)>=6)
break;
5838 ckm_takeout = pow(ckm_ir*ckm_js, 2);
5839 for (
unsigned int iswap=0; iswap<2; iswap++){
5840 msq_uubar_ww_ijrs1234[iswap] /= ckm_takeout;
5841 msq_uubar_ww_ijrs1243[iswap] /= ckm_takeout;
5848 (partonIsUnknown[0] && partonIsUnknown[1])
5858 (partonIsUnknown[2] && partonIsUnknown[3])
5866 rsel=isel; ssel=jsel;
5873 (partonIsUnknown[2] && partonIsUnknown[3])
5892 if (abs(isel)>=6)
break;
5901 ckm_takeout = pow(ckm_ir*ckm_js, 2);
5902 for (
unsigned int iswap=0; iswap<2; iswap++){
5903 msq_ddbar_ww_ijrs1234[iswap] /= ckm_takeout;
5904 msq_ddbar_ww_ijrs1243[iswap] /= ckm_takeout;
5912 (partonIsUnknown[0] && partonIsUnknown[1])
5926 (partonIsUnknown[2] && partonIsUnknown[3])
5942 const unsigned int nCombinations = 6;
5943 int ijrs_arr[nCombinations][4] ={
5951 for (
unsigned int ijrs=0; ijrs<nCombinations; ijrs++){
5952 isel = ijrs_arr[ijrs][0];
5953 jsel = ijrs_arr[ijrs][1];
5954 rsel = ijrs_arr[ijrs][2];
5955 ssel = ijrs_arr[ijrs][3];
5958 if (ckm_ir!=0. && ckm_js!=0.)
break;
5960 if (ckm_ir!=0. && ckm_js!=0.){
5965 ckm_takeout = pow(ckm_ir*ckm_js, 2);
5966 for (
unsigned int iswap=0; iswap<2; iswap++){
5967 msq_ud_wwonly_ijrs1234[iswap] /= ckm_takeout;
5968 msq_ud_wwonly_ijrs1243[iswap] /= ckm_takeout;
5975 (partonIsUnknown[0] && partonIsUnknown[1])
5989 (partonIsUnknown[2] && partonIsUnknown[3])
6005 const unsigned int nCombinations = 6;
6006 int ijrs_arr[nCombinations][4] ={
6014 for (
unsigned int ijrs=0; ijrs<nCombinations; ijrs++){
6015 isel = ijrs_arr[ijrs][0];
6016 jsel = ijrs_arr[ijrs][1];
6017 rsel = ijrs_arr[ijrs][2];
6018 ssel = ijrs_arr[ijrs][3];
6021 if (ckm_ir!=0. && ckm_js!=0.)
break;
6023 if (ckm_ir!=0. && ckm_js!=0.){
6028 ckm_takeout = pow(ckm_ir*ckm_js, 2);
6029 for (
unsigned int iswap=0; iswap<2; iswap++){
6030 msq_ubardbar_wwonly_ijrs1234[iswap] /= ckm_takeout;
6031 msq_ubardbar_wwonly_ijrs1243[iswap] /= ckm_takeout;
6037 MELAout <<
"TUtil::HJJMatEl: The pre-computed MEs:" << endl;
6038 for (
unsigned int iswap=0; iswap<2; iswap++){
6040 <<
"\tmsq_uu_zz_ijrs1234[" << iswap <<
"] = " << msq_uu_zz_ijrs1234[iswap] <<
'\n'
6041 <<
"\tmsq_uu_zz_ijrs1243[" << iswap <<
"] = " << msq_uu_zz_ijrs1243[iswap] <<
'\n'
6042 <<
"\tmsq_dd_zz_ijrs1234[" << iswap <<
"] = " << msq_dd_zz_ijrs1234[iswap] <<
'\n'
6043 <<
"\tmsq_dd_zz_ijrs1243[" << iswap <<
"] = " << msq_dd_zz_ijrs1243[iswap] <<
'\n'
6044 <<
"\tmsq_ubarubar_zz_ijrs1234[" << iswap <<
"] = " << msq_ubarubar_zz_ijrs1234[iswap] <<
'\n'
6045 <<
"\tmsq_ubarubar_zz_ijrs1243[" << iswap <<
"] = " << msq_ubarubar_zz_ijrs1243[iswap] <<
'\n'
6046 <<
"\tmsq_dbardbar_zz_ijrs1234[" << iswap <<
"] = " << msq_dbardbar_zz_ijrs1234[iswap] <<
'\n'
6047 <<
"\tmsq_dbardbar_zz_ijrs1243[" << iswap <<
"] = " << msq_dbardbar_zz_ijrs1243[iswap] <<
'\n'
6048 <<
"\tmsq_uu_zzid_ijrs1234[" << iswap <<
"] = " << msq_uu_zzid_ijrs1234[iswap] <<
'\n'
6049 <<
"\tmsq_uu_zzid_ijrs1243[" << iswap <<
"] = " << msq_uu_zzid_ijrs1243[iswap] <<
'\n'
6050 <<
"\tmsq_dd_zzid_ijrs1234[" << iswap <<
"] = " << msq_dd_zzid_ijrs1234[iswap] <<
'\n'
6051 <<
"\tmsq_dd_zzid_ijrs1243[" << iswap <<
"] = " << msq_dd_zzid_ijrs1243[iswap] <<
'\n'
6052 <<
"\tmsq_ubarubar_zzid_ijrs1234[" << iswap <<
"] = " << msq_ubarubar_zzid_ijrs1234[iswap] <<
'\n'
6053 <<
"\tmsq_ubarubar_zzid_ijrs1243[" << iswap <<
"] = " << msq_ubarubar_zzid_ijrs1243[iswap] <<
'\n'
6054 <<
"\tmsq_dbardbar_zzid_ijrs1234[" << iswap <<
"] = " << msq_dbardbar_zzid_ijrs1234[iswap] <<
'\n'
6055 <<
"\tmsq_dbardbar_zzid_ijrs1243[" << iswap <<
"] = " << msq_dbardbar_zzid_ijrs1243[iswap] <<
'\n'
6056 <<
"\tmsq_udbar_zz_ijrs1234[" << iswap <<
"] = " << msq_udbar_zz_ijrs1234[iswap] <<
'\n'
6057 <<
"\tmsq_udbar_zz_ijrs1243[" << iswap <<
"] = " << msq_udbar_zz_ijrs1243[iswap] <<
'\n'
6058 <<
"\tmsq_dubar_zz_ijrs1234[" << iswap <<
"] = " << msq_dubar_zz_ijrs1234[iswap] <<
'\n'
6059 <<
"\tmsq_dubar_zz_ijrs1243[" << iswap <<
"] = " << msq_dubar_zz_ijrs1243[iswap] <<
'\n'
6060 <<
"\tmsq_uubar_zz_ijrs1234[" << iswap <<
"] = " << msq_uubar_zz_ijrs1234[iswap] <<
'\n'
6061 <<
"\tmsq_uubar_zz_ijrs1243[" << iswap <<
"] = " << msq_uubar_zz_ijrs1243[iswap] <<
'\n'
6062 <<
"\tmsq_ddbar_zz_ijrs1234[" << iswap <<
"] = " << msq_ddbar_zz_ijrs1234[iswap] <<
'\n'
6063 <<
"\tmsq_ddbar_zz_ijrs1243[" << iswap <<
"] = " << msq_ddbar_zz_ijrs1243[iswap] <<
'\n'
6064 <<
"\tmsq_uubar_ww_ijrs1234[" << iswap <<
"] = " << msq_uubar_ww_ijrs1234[iswap] <<
'\n'
6065 <<
"\tmsq_uubar_ww_ijrs1243[" << iswap <<
"] = " << msq_uubar_ww_ijrs1243[iswap] <<
'\n'
6066 <<
"\tmsq_ddbar_ww_ijrs1234[" << iswap <<
"] = " << msq_ddbar_ww_ijrs1234[iswap] <<
'\n'
6067 <<
"\tmsq_ddbar_ww_ijrs1243[" << iswap <<
"] = " << msq_ddbar_ww_ijrs1243[iswap] <<
'\n'
6068 <<
"\tmsq_ud_wwonly_ijrs1234[" << iswap <<
"] = " << msq_ud_wwonly_ijrs1234[iswap] <<
'\n'
6069 <<
"\tmsq_ud_wwonly_ijrs1243[" << iswap <<
"] = " << msq_ud_wwonly_ijrs1243[iswap] <<
'\n'
6070 <<
"\tmsq_ubardbar_wwonly_ijrs1234[" << iswap <<
"] = " << msq_ubardbar_wwonly_ijrs1234[iswap] <<
'\n'
6071 <<
"\tmsq_ubardbar_wwonly_ijrs1243[" << iswap <<
"] = " << msq_ubardbar_wwonly_ijrs1243[iswap]
6077 int nijchannels = ijsel.size();
6079 for (
int ic=0; ic<nijchannels; ic++){
6081 isel = ijsel[ic][0];
6082 jsel = ijsel[ic][1];
6083 int code = ijsel[ic][2];
6084 bool ijselIsUpType[2];
6085 bool ijselIsDownType[2];
6086 bool ijselIsParticle[2];
6087 bool ijselIsAntiparticle[2];
6101 ijselIsParticle[0] = (isel>0); ijselIsAntiparticle[0] = (isel<0);
6102 ijselIsParticle[1] = (jsel>0); ijselIsAntiparticle[1] = (jsel<0);
6105 (partonIsUnknown[0] || MYIDUP_tmp[0]==isel)
6107 (partonIsUnknown[1] || MYIDUP_tmp[1]==jsel)
6112 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6114 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6117 if (ijselIsUpType[0] && ijselIsUpType[1]){
6118 if (ijselIsParticle[0] && ijselIsParticle[1]){
6119 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1234[0];
6120 else msq_tmp = msq_uu_zzid_ijrs1234[0];
6122 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6123 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1234[0];
6124 else msq_tmp = msq_ubarubar_zzid_ijrs1234[0];
6126 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_zz_ijrs1234[0];
6128 else if (ijselIsDownType[0] && ijselIsDownType[1]){
6129 if (ijselIsParticle[0] && ijselIsParticle[1]){
6130 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1234[0];
6131 else msq_tmp = msq_dd_zzid_ijrs1234[0];
6133 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6134 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1234[0];
6135 else msq_tmp = msq_dbardbar_zzid_ijrs1234[0];
6137 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_zz_ijrs1234[0];
6139 else if (ijselIsUpType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_udbar_zz_ijrs1234[0];
6140 else if (ijselIsDownType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_dubar_zz_ijrs1234[0];
6141 MatElsq[jsel+5][isel+5] += msq_tmp;
6142 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6147 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6149 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6152 if (ijselIsUpType[0] && ijselIsUpType[1]){
6153 if (ijselIsParticle[0] && ijselIsParticle[1]){
6154 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1243[0];
6155 else msq_tmp = msq_uu_zzid_ijrs1243[0];
6157 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6158 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1243[0];
6159 else msq_tmp = msq_ubarubar_zzid_ijrs1243[0];
6161 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_zz_ijrs1243[0];
6163 else if (ijselIsDownType[0] && ijselIsDownType[1]){
6164 if (ijselIsParticle[0] && ijselIsParticle[1]){
6165 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1243[0];
6166 else msq_tmp = msq_dd_zzid_ijrs1243[0];
6168 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]){
6169 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1243[0];
6170 else msq_tmp = msq_dbardbar_zzid_ijrs1243[0];
6172 else if (ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_zz_ijrs1243[0];
6174 else if (ijselIsUpType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_udbar_zz_ijrs1243[0];
6175 else if (ijselIsDownType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_dubar_zz_ijrs1243[0];
6176 MatElsq[jsel+5][isel+5] += msq_tmp;
6177 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6181 vector<int> possible_rsel;
6182 vector<int> possible_ssel;
6183 vector<double> possible_Vsqir;
6184 vector<double> possible_Vsqjs;
6185 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6186 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6187 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6188 rsel=possible_rsel.at(ix);
6190 possible_Vsqir.push_back(ckmval);
6192 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6193 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6194 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6195 ssel=possible_ssel.at(iy);
6197 possible_Vsqjs.push_back(ckmval);
6201 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6202 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6203 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6204 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6205 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6207 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6209 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6212 if (ijselIsUpType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_ww_ijrs1234[0];
6213 else if (ijselIsDownType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_ww_ijrs1234[0];
6215 MatElsq[jsel+5][isel+5] += msq_tmp;
6216 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6221 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6223 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6226 if (ijselIsUpType[0] && ijselIsUpType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_uubar_ww_ijrs1243[0];
6227 else if (ijselIsDownType[0] && ijselIsDownType[1] && ijselIsParticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ddbar_ww_ijrs1243[0];
6229 MatElsq[jsel+5][isel+5] += msq_tmp;
6230 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6236 vector<int> possible_rsel;
6237 vector<int> possible_ssel;
6238 vector<double> possible_Vsqir;
6239 vector<double> possible_Vsqjs;
6240 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6241 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6242 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6243 rsel=possible_rsel.at(ix);
6245 possible_Vsqir.push_back(ckmval);
6247 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6248 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6249 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6250 ssel=possible_ssel.at(iy);
6252 possible_Vsqjs.push_back(ckmval);
6256 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6257 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6258 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6259 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6260 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6262 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6264 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6266 if (rsel!=jsel && ssel!=isel){
6268 if (ijselIsUpType[0] && ijselIsDownType[1]){
6269 if (ijselIsParticle[0] && ijselIsParticle[1]) msq_tmp = msq_ud_wwonly_ijrs1234[0];
6270 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ubardbar_wwonly_ijrs1234[0];
6273 MatElsq[jsel+5][isel+5] += msq_tmp;
6274 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6278 MatElsq[jsel+5][isel+5] += msq_tmp;
6279 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6285 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6287 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6289 if (rsel!=jsel && ssel!=isel){
6291 if (ijselIsUpType[0] && ijselIsDownType[1]){
6292 if (ijselIsParticle[0] && ijselIsParticle[1]) msq_tmp = msq_ud_wwonly_ijrs1243[0];
6293 else if (ijselIsAntiparticle[0] && ijselIsAntiparticle[1]) msq_tmp = msq_ubardbar_wwonly_ijrs1243[0];
6296 MatElsq[jsel+5][isel+5] += msq_tmp;
6297 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6301 MatElsq[jsel+5][isel+5] += msq_tmp;
6302 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6310 if (isel==jsel)
continue;
6311 isel = ijsel[ic][1];
6312 jsel = ijsel[ic][0];
6324 ijselIsParticle[0] = (isel>0); ijselIsAntiparticle[0] = (isel<0);
6325 ijselIsParticle[1] = (jsel>0); ijselIsAntiparticle[1] = (jsel<0);
6328 (partonIsUnknown[0] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[0]) && isel==0) || MYIDUP_tmp[0]==isel))
6330 (partonIsUnknown[1] || ((
PDGHelpers::isAGluon(MYIDUP_tmp[1]) && jsel==0) || MYIDUP_tmp[1]==jsel))
6336 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6338 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6341 if (ijselIsUpType[1] && ijselIsUpType[0]){
6342 if (ijselIsParticle[1] && ijselIsParticle[0]){
6343 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1234[1];
6344 else msq_tmp = msq_uu_zzid_ijrs1234[1];
6346 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6347 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1234[1];
6348 else msq_tmp = msq_ubarubar_zzid_ijrs1234[1];
6350 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_zz_ijrs1234[1];
6352 else if (ijselIsDownType[1] && ijselIsDownType[0]){
6353 if (ijselIsParticle[1] && ijselIsParticle[0]){
6354 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1234[1];
6355 else msq_tmp = msq_dd_zzid_ijrs1234[1];
6357 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6358 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1234[1];
6359 else msq_tmp = msq_dbardbar_zzid_ijrs1234[1];
6361 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_zz_ijrs1234[1];
6363 else if (ijselIsUpType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_udbar_zz_ijrs1234[1];
6364 else if (ijselIsDownType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_dubar_zz_ijrs1234[1];
6365 MatElsq[jsel+5][isel+5] += msq_tmp;
6366 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6371 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6373 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6376 if (ijselIsUpType[1] && ijselIsUpType[0]){
6377 if (ijselIsParticle[1] && ijselIsParticle[0]){
6378 if (isel!=jsel) msq_tmp = msq_uu_zz_ijrs1243[1];
6379 else msq_tmp = msq_uu_zzid_ijrs1243[1];
6381 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6382 if (isel!=jsel) msq_tmp = msq_ubarubar_zz_ijrs1243[1];
6383 else msq_tmp = msq_ubarubar_zzid_ijrs1243[1];
6385 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_zz_ijrs1243[1];
6387 else if (ijselIsDownType[1] && ijselIsDownType[0]){
6388 if (ijselIsParticle[1] && ijselIsParticle[0]){
6389 if (isel!=jsel) msq_tmp = msq_dd_zz_ijrs1243[1];
6390 else msq_tmp = msq_dd_zzid_ijrs1243[1];
6392 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]){
6393 if (isel!=jsel) msq_tmp = msq_dbardbar_zz_ijrs1243[1];
6394 else msq_tmp = msq_dbardbar_zzid_ijrs1243[1];
6396 else if (ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_zz_ijrs1243[1];
6398 else if (ijselIsUpType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_udbar_zz_ijrs1243[1];
6399 else if (ijselIsDownType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_dubar_zz_ijrs1243[1];
6400 MatElsq[jsel+5][isel+5] += msq_tmp;
6401 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6405 vector<int> possible_rsel;
6406 vector<int> possible_ssel;
6407 vector<double> possible_Vsqir;
6408 vector<double> possible_Vsqjs;
6409 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6410 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6411 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6412 rsel=possible_rsel.at(ix);
6414 possible_Vsqir.push_back(ckmval);
6416 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6417 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6418 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6419 ssel=possible_ssel.at(iy);
6421 possible_Vsqjs.push_back(ckmval);
6425 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6426 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6427 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6428 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6429 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6431 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6433 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6436 if (ijselIsUpType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_ww_ijrs1234[1];
6437 else if (ijselIsDownType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_ww_ijrs1234[1];
6439 MatElsq[jsel+5][isel+5] += msq_tmp;
6440 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6445 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6447 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6450 if (ijselIsUpType[1] && ijselIsUpType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_uubar_ww_ijrs1243[1];
6451 else if (ijselIsDownType[1] && ijselIsDownType[0] && ijselIsParticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ddbar_ww_ijrs1243[1];
6453 MatElsq[jsel+5][isel+5] += msq_tmp;
6454 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6460 vector<int> possible_rsel;
6461 vector<int> possible_ssel;
6462 vector<double> possible_Vsqir;
6463 vector<double> possible_Vsqjs;
6464 if (ijselIsUpType[0]){ possible_rsel.push_back(1); possible_rsel.push_back(3); possible_rsel.push_back(5); }
6465 else if (ijselIsDownType[0]){ possible_rsel.push_back(2); possible_rsel.push_back(4); }
6466 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6467 rsel=possible_rsel.at(ix);
6469 possible_Vsqir.push_back(ckmval);
6471 if (ijselIsUpType[1]){ possible_ssel.push_back(1); possible_ssel.push_back(3); possible_ssel.push_back(5); }
6472 else if (ijselIsDownType[1]){ possible_ssel.push_back(2); possible_ssel.push_back(4); }
6473 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6474 ssel=possible_ssel.at(iy);
6476 possible_Vsqjs.push_back(ckmval);
6480 for (
unsigned int ix=0; ix<possible_rsel.size(); ix++){
6481 rsel=possible_rsel.at(ix)*TMath::Sign(1, isel);
6482 for (
unsigned int iy=0; iy<possible_ssel.size(); iy++){
6483 ssel=possible_ssel.at(iy)*TMath::Sign(1, jsel);
6484 double ckmval = possible_Vsqir.at(ix)*possible_Vsqjs.at(iy);
6486 (partonIsUnknown[2] || MYIDUP_tmp[2]==rsel)
6488 (partonIsUnknown[3] || MYIDUP_tmp[3]==ssel)
6490 if (rsel!=jsel && ssel!=isel){
6492 if (ijselIsUpType[1] && ijselIsDownType[0]){
6493 if (ijselIsParticle[1] && ijselIsParticle[0]) msq_tmp = msq_ud_wwonly_ijrs1234[1];
6494 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ubardbar_wwonly_ijrs1234[1];
6497 MatElsq[jsel+5][isel+5] += msq_tmp;
6498 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6502 MatElsq[jsel+5][isel+5] += msq_tmp;
6503 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << rsel <<
", " << ssel <<
'\t' << msq_tmp << endl;
6509 (partonIsUnknown[2] || MYIDUP_tmp[2]==ssel)
6511 (partonIsUnknown[3] || MYIDUP_tmp[3]==rsel)
6513 if (rsel!=jsel && ssel!=isel){
6515 if (ijselIsUpType[1] && ijselIsDownType[0]){
6516 if (ijselIsParticle[1] && ijselIsParticle[0]) msq_tmp = msq_ud_wwonly_ijrs1243[1];
6517 else if (ijselIsAntiparticle[1] && ijselIsAntiparticle[0]) msq_tmp = msq_ubardbar_wwonly_ijrs1243[1];
6520 MatElsq[jsel+5][isel+5] += msq_tmp;
6521 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6525 MatElsq[jsel+5][isel+5] += msq_tmp;
6526 if (verbosity >=
TVar::DEBUG_VERBOSE)
MELAout <<
"Channel (isel, jsel, rsel, ssel)=" << isel <<
", " << jsel <<
", " << ssel <<
", " << rsel <<
'\t' << msq_tmp << endl;
6537 int GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets)*2;
6538 double constant = pow(GeV, -GeVexponent_MEsq);
6539 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
6545 for (
int ii = 0; ii <
nmsq; ii++){
for (
int jj = 0; jj <
nmsq; jj++)
MELAout << MatElsq[jj][ii] <<
'\t';
MELAout << endl; }
6547 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
6564 <<
"TUtil::HJJMatEl: Reset AlphaS:\n"
6565 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
6567 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
6571 <<
"TUtil::HJJMatEl: Reset AlphaS result:\n"
6572 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
6581 const double& EBEAM,
6582 bool includeHiggsDecay,
6585 const double GeV=1./100.;
6586 double sum_msqjk = 0;
6594 double MatElsq[
nmsq][
nmsq]={ { 0 } };
6596 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::VHiggsMatEl: Non-JHUGen MEs are not supported" << endl;
return sum_msqjk; }
6599 int nRequested_AssociatedJets=0;
6600 int nRequested_AssociatedLeptons=0;
6601 int nRequested_AssociatedPhotons=0;
6602 int AssociationVCompatibility=0;
6606 nRequested_AssociatedJets=2;
6610 nRequested_AssociatedLeptons=2;
6614 nRequested_AssociatedPhotons=1;
6632 MELAerr <<
"TUtil::VHiggsMatEl: Number of associated particles (" << mela_event.
pAssociated.size() <<
") is less than ";
6634 else MELAerr << nRequested_AssociatedPhotons;
6640 int MYIDUP_prod[4]={ 0 };
6641 int MYIDUP_dec[2]={ -9000, -9000 };
6642 double p4[9][4] ={ { 0 } };
6643 double helicities[9] ={ 0 };
6644 int vh_ids[9] ={ 0 };
6645 TLorentzVector MomStore[
mxpart];
6646 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
6654 for (
int ipar=0; ipar<2; ipar++){
6655 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
6656 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
6658 else MYIDUP_prod[ipar] = 0;
6659 if (momTmp->T()>0.){
6660 p4[ipar][0] = momTmp->T()*GeV;
6661 p4[ipar][1] = momTmp->X()*GeV;
6662 p4[ipar][2] = momTmp->Y()*GeV;
6663 p4[ipar][3] = momTmp->Z()*GeV;
6664 MomStore[ipar] = (*momTmp);
6667 p4[ipar][0] = -momTmp->T()*GeV;
6668 p4[ipar][1] = -momTmp->X()*GeV;
6669 p4[ipar][2] = -momTmp->Y()*GeV;
6670 p4[ipar][3] = -momTmp->Z()*GeV;
6671 MomStore[ipar] = -(*momTmp);
6672 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
6677 TLorentzVector* momTmp = &(mela_event.
pAssociated.at(ipar).second);
6678 int* idtmp = &(mela_event.
pAssociated.at(ipar).first);
6680 else MYIDUP_prod[ipar+2] = 0;
6681 p4[ipar+5][0] = momTmp->T()*GeV;
6682 p4[ipar+5][1] = momTmp->X()*GeV;
6683 p4[ipar+5][2] = momTmp->Y()*GeV;
6684 p4[ipar+5][3] = momTmp->Z()*GeV;
6685 MomStore[ipar+6] = (*momTmp);
6695 for (
int iv=0; iv<2; iv++){
6698 else MYIDUP_dec[iv] = 0;
6701 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
6702 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
6704 p4[ipar+7][0] = momTmp->T()*GeV;
6705 p4[ipar+7][1] = momTmp->X()*GeV;
6706 p4[ipar+7][2] = momTmp->Y()*GeV;
6707 p4[ipar+7][3] = momTmp->Z()*GeV;
6708 MomStore[5] = (*momTmp);
6711 p4[ipar+7][0] = momTmp->T()*GeV;
6712 p4[ipar+7][1] = momTmp->X()*GeV;
6713 p4[ipar+7][2] = momTmp->Y()*GeV;
6714 p4[ipar+7][3] = momTmp->Z()*GeV;
6715 MomStore[2*ipar+2] = (*momTmp);
6720 p4[7][0] += momTmp->T()*GeV;
6721 p4[7][1] += momTmp->X()*GeV;
6722 p4[7][2] += momTmp->Y()*GeV;
6723 p4[7][3] += momTmp->Z()*GeV;
6726 p4[8][0] = momTmp->T()*GeV;
6727 p4[8][1] = momTmp->X()*GeV;
6728 p4[8][2] = momTmp->Y()*GeV;
6729 p4[8][3] = momTmp->Z()*GeV;
6731 MomStore[ipar+2] = (*momTmp);
6735 p4[7][0] += momTmp->T()*GeV;
6736 p4[7][1] += momTmp->X()*GeV;
6737 p4[7][2] += momTmp->Y()*GeV;
6738 p4[7][3] += momTmp->Z()*GeV;
6741 p4[8][0] += momTmp->T()*GeV;
6742 p4[8][1] += momTmp->X()*GeV;
6743 p4[8][2] += momTmp->Y()*GeV;
6744 p4[8][3] += momTmp->Z()*GeV;
6746 MomStore[ipar+2] = (*momTmp);
6749 p4[7][0] += momTmp->T()*GeV;
6750 p4[7][1] += momTmp->X()*GeV;
6751 p4[7][2] += momTmp->Y()*GeV;
6752 p4[7][3] += momTmp->Z()*GeV;
6753 MomStore[5] = MomStore[5] + (*momTmp);
6756 for (
int ix=0; ix<4; ix++){
6757 p4[3][ix] = p4[5][ix] + p4[6][ix];
6758 p4[4][ix] = p4[7][ix] + p4[8][ix];
6759 p4[2][ix] = p4[3][ix] + p4[4][ix];
6766 else vh_ids[3] = 22;
6769 vh_ids[7] = 5; helicities[7] = 1;
6770 vh_ids[8] = -5; helicities[8] = 1;
6772 if (includeHiggsDecay && MYIDUP_dec[0]!=-9000 && MYIDUP_dec[1]!=-9000 && MYIDUP_dec[0]==-MYIDUP_dec[1]){
6777 else if (verbosity>=
TVar::INFO && includeHiggsDecay)
MELAerr <<
"TUtil::VHiggsMatEl: includeHiggsDecay=true is not supported for the present decay mode." << endl;
6780 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";
6783 double defaultRenScale =
scale_.scale;
6784 double defaultFacScale =
facscale_.facscale;
6786 int defaultNflav =
nflav_.nflav;
6787 string defaultPdflabel =
pdlabel_.pdlabel;
6791 double alphasVal, alphasmzVal;
6801 <<
"TUtil::VHiggsMatEl: Set AlphaS:\n"
6802 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
6804 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
6808 bool partonIsKnown[4];
6809 for (
unsigned int ip=0; ip<4; ip++) partonIsKnown[ip] = (MYIDUP_prod[ip]!=0);
6812 const double allowed_helicities[2] ={ -1, 1 };
6814 vector<pair<int, int>> Hffparticles;
6819 else Hffparticles.push_back(pair<int, int>(-MYIDUP_dec[1], MYIDUP_dec[1]));
6822 double Hffscalesum=0;
6823 for (
int hquark=1; hquark<=5; hquark++){
6825 Hffscalesum += Hffmass;
6827 Hffparticles.push_back(pair<int, int>(hquark, -hquark));
6828 Hffparticles.push_back(pair<int, int>(-hquark, hquark));
6829 Hffscale /= Hffmass;
6832 Hffscale *= Hffscalesum;
6836 MELAout <<
"TUtil::VHiggsMatEl: Outgoing H-> f fbar particles to compute for the ME template:" << endl;
6837 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++)
MELAout <<
"\t - (id8, id9) = (" << Hffparticles.at(ihf).first <<
", " << Hffparticles.at(ihf).second <<
")" << endl;
6838 MELAout <<
"TUtil::VHiggsMatEl: ME scale for the H-> f fbar particles: " << Hffscale << endl;
6843 vector<pair<int, int>> incomingPartons;
6844 if (partonIsKnown[0] && partonIsKnown[1]) incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], MYIDUP_prod[1]));
6846 else if (!partonIsKnown[0] && !partonIsKnown[1]){
6848 incomingPartons.push_back(pair<int, int>(1, -2));
6849 incomingPartons.push_back(pair<int, int>(-2, 1));
6850 incomingPartons.push_back(pair<int, int>(2, -1));
6851 incomingPartons.push_back(pair<int, int>(-1, 2));
6853 else if (!partonIsKnown[1] && partonIsKnown[0]){
6856 for (
int iqf=1; iqf<=5; iqf++){
6857 if (iqf%2==0)
continue;
6858 int jqf = -TMath::Sign(iqf, MYIDUP_prod[0]);
6860 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], jqf));
break; }
6864 for (
int iqf=2; iqf<=4; iqf++){
6865 if (iqf%2==1)
continue;
6866 int jqf = -TMath::Sign(iqf, MYIDUP_prod[0]);
6868 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], jqf));
break; }
6875 for (
int iqf=1; iqf<=5; iqf++){
6876 if (iqf%2==0)
continue;
6877 int jqf = -TMath::Sign(iqf, MYIDUP_prod[1]);
6879 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[1]));
break; }
6883 for (
int iqf=2; iqf<=4; iqf++){
6884 if (iqf%2==1)
continue;
6885 int jqf = -TMath::Sign(iqf, MYIDUP_prod[1]);
6887 if (ckm_test!=0.){ incomingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[1]));
break; }
6892 MELAout <<
"TUtil::VHiggsMatEl: Incoming partons to compute for the ME template:" << endl;
6893 for (
auto& inpart:incomingPartons)
MELAout <<
"\t - (id1, id2) = (" << inpart.first <<
", " << inpart.second <<
")" << endl;
6897 vector<pair<int, int>> outgoingPartons;
6898 if (partonIsKnown[2] && partonIsKnown[3]) outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], MYIDUP_prod[3]));
6900 else if (!partonIsKnown[2] && !partonIsKnown[3]){
6902 outgoingPartons.push_back(pair<int, int>(1, -2));
6903 outgoingPartons.push_back(pair<int, int>(-2, 1));
6904 outgoingPartons.push_back(pair<int, int>(2, -1));
6905 outgoingPartons.push_back(pair<int, int>(-1, 2));
6907 else if (!partonIsKnown[3] && partonIsKnown[2]){
6910 for (
int iqf=1; iqf<=5; iqf++){
6911 if (iqf%2==0)
continue;
6912 int jqf = -TMath::Sign(iqf, MYIDUP_prod[2]);
6914 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], jqf));
break; }
6918 for (
int iqf=2; iqf<=4; iqf++){
6919 if (iqf%2==1)
continue;
6920 int jqf = -TMath::Sign(iqf, MYIDUP_prod[2]);
6922 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], jqf));
break; }
6929 for (
int iqf=1; iqf<=5; iqf++){
6930 if (iqf%2==0)
continue;
6931 int jqf = -TMath::Sign(iqf, MYIDUP_prod[3]);
6933 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[3]));
break; }
6937 for (
int iqf=2; iqf<=4; iqf++){
6938 if (iqf%2==1)
continue;
6939 int jqf = -TMath::Sign(iqf, MYIDUP_prod[3]);
6941 if (ckm_test!=0.){ outgoingPartons.push_back(pair<int, int>(jqf, MYIDUP_prod[3]));
break; }
6946 MELAout <<
"TUtil::VHiggsMatEl: Outgoing particles to compute for the ME template:" << endl;
6947 for (
auto& outpart:outgoingPartons)
MELAout <<
"\t - (id6, id7) = (" << outpart.first <<
", " << outpart.second <<
")" << endl;
6950 for (
auto& inpart:incomingPartons){
6951 vh_ids[0] = inpart.first;
6952 vh_ids[1] = inpart.second;
6955 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tIncoming " << vh_ids[0] <<
"," << vh_ids[1] <<
" -> " << vh_ids[2] << endl;
6958 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tNeed to divide the ME by |VCKM_incoming|**2 = " << Vckmsq_in << endl;
6960 for (
auto& outpart:outgoingPartons){
6961 vh_ids[5] = outpart.first;
6962 vh_ids[6] = outpart.second;
6964 if (vh_ids[2]!=vh_ids[3])
continue;
6965 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tOutgoing " << vh_ids[3] <<
" -> " << vh_ids[5] <<
"," << vh_ids[6] << endl;
6969 for (
int h01 = 0; h01 < 2; h01++){
6970 helicities[0] = allowed_helicities[h01];
6971 helicities[1] = -helicities[0];
6972 for (
int h56 = 0; h56 < 2; h56++){
6973 helicities[5] = allowed_helicities[h56];
6974 helicities[6] = -helicities[5];
6979 for (
int h78=0; h78<2; h78++){
6980 helicities[7]=allowed_helicities[h78];
6981 helicities[8]=allowed_helicities[h78];
6982 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++){
6983 vh_ids[7]=Hffparticles.at(ihf).first;
6984 vh_ids[8]=Hffparticles.at(ihf).second;
6985 double msq_inst_LR=0;
6987 msq_inst += msq_inst_LR;
6990 msq_inst *= Hffscale;
6997 double scalesum_out=0;
6998 if (!(partonIsKnown[2] && partonIsKnown[3])){
7000 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tDividing ME by |VCKM_outgoing|**2 = " << Vckmsq_out << endl;
7002 for (
int outgoing1=1; outgoing1<=
nf; outgoing1++){
7003 if (partonIsKnown[2] && outgoing1!=abs(vh_ids[5]))
continue;
7004 if (outgoing1%2!=abs(vh_ids[5])%2 || outgoing1==6)
continue;
7005 int iout = outgoing1 * TMath::Sign(1, vh_ids[5]);
7006 for (
int outgoing2=1; outgoing2<=
nf; outgoing2++){
7007 if (partonIsKnown[3] && outgoing2!=abs(vh_ids[6]))
continue;
7008 if (outgoing2%2!=abs(vh_ids[6])%2 || outgoing2==6)
continue;
7009 int jout = outgoing2 * TMath::Sign(1, vh_ids[6]);
7014 else scalesum_out = 1;
7015 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for outgoing particles: " << scalesum_out << endl;
7018 if (!partonIsKnown[0] || !partonIsKnown[1]){
7019 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tDividing ME by |VCKM_incoming|**2 = " << Vckmsq_in << endl;
7024 for (
int incoming1=1; incoming1<=
nf; incoming1++){
7025 if (partonIsKnown[0] && incoming1!=abs(vh_ids[0]))
continue;
7026 if (incoming1%2!=abs(vh_ids[0])%2 || incoming1==6)
continue;
7027 int iin = incoming1 * TMath::Sign(1, vh_ids[0]);
7028 for (
int incoming2=1; incoming2<=
nf; incoming2++){
7029 if (partonIsKnown[1] && incoming2!=abs(vh_ids[1]))
continue;
7030 if (incoming2%2!=abs(vh_ids[1])%2 || incoming2==6)
continue;
7031 int jin = incoming2 * TMath::Sign(1, vh_ids[1]);
7032 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\t\tiin,jin = " << iin <<
"," << jin << endl;
7034 if (!partonIsKnown[0] || !partonIsKnown[1]){
7036 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for incoming particles: " << scale_in << endl;
7038 MatElsq[jin+5][iin+5] +=
msq * 0.25 * scale_in *scalesum_out;
7047 vector<pair<int, int>> incomingPartons;
7048 if (partonIsKnown[0] && partonIsKnown[1]) incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], MYIDUP_prod[1]));
7049 else if (!partonIsKnown[0] && !partonIsKnown[1]){
7051 incomingPartons.push_back(pair<int, int>(1, -1));
7052 incomingPartons.push_back(pair<int, int>(-1, 1));
7053 incomingPartons.push_back(pair<int, int>(2, -2));
7054 incomingPartons.push_back(pair<int, int>(-2, 2));
7056 else if (!partonIsKnown[1] && partonIsKnown[0])
7057 incomingPartons.push_back(pair<int, int>(MYIDUP_prod[0], -MYIDUP_prod[0]));
7059 incomingPartons.push_back(pair<int, int>(-MYIDUP_prod[1], MYIDUP_prod[1]));
7061 MELAout <<
"TUtil::VHiggsMatEl: Incoming partons to compute for the ME template:" << endl;
7062 for (
auto& inpart:incomingPartons)
MELAout <<
"\t - (id1, id2) = (" << inpart.first <<
", " << inpart.second <<
")" << endl;
7066 vector<pair<int, int>> outgoingPartons;
7067 if ((partonIsKnown[2] && partonIsKnown[3]) ||
production ==
TVar::GammaH) outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], MYIDUP_prod[3]));
7068 else if (!partonIsKnown[2] && !partonIsKnown[3]){
7070 outgoingPartons.push_back(pair<int, int>(1, -1));
7071 outgoingPartons.push_back(pair<int, int>(-1, 1));
7072 outgoingPartons.push_back(pair<int, int>(2, -2));
7073 outgoingPartons.push_back(pair<int, int>(-2, 2));
7075 else if (!partonIsKnown[3] && partonIsKnown[2])
7076 outgoingPartons.push_back(pair<int, int>(MYIDUP_prod[2], -MYIDUP_prod[2]));
7078 outgoingPartons.push_back(pair<int, int>(-MYIDUP_prod[3], MYIDUP_prod[3]));
7080 MELAout <<
"TUtil::VHiggsMatEl: Outgoing particles to compute for the ME template:" << endl;
7081 for (
unsigned int op=0; op<outgoingPartons.size(); op++)
MELAout <<
"\t - (id6, id7) = (" << outgoingPartons.at(op).first <<
", " << outgoingPartons.at(op).second <<
")" << endl;
7084 for (
auto& inpart:incomingPartons){
7085 vh_ids[0] = inpart.first;
7086 vh_ids[1] = inpart.second;
7089 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\tIncoming " << vh_ids[0] <<
"," << vh_ids[1] <<
" -> " << vh_ids[2] << endl;
7093 for (
auto& outpart:outgoingPartons){
7094 vh_ids[5] = outpart.first;
7095 vh_ids[6] = outpart.second;
7099 if (vh_ids[2]!=vh_ids[3])
continue;
7101 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tOutgoing " << vh_ids[3] <<
" -> " << vh_ids[5] <<
"," << vh_ids[6] << endl;
7105 for (
int h01 = 0; h01 < 2; h01++){
7106 helicities[0] = allowed_helicities[h01];
7107 helicities[1] = -helicities[0];
7108 for (
int h56 = 0; h56 < 2; h56++){
7109 helicities[5] = allowed_helicities[h56];
7110 helicities[6] = -helicities[5];
7115 for (
int h78=0; h78<2; h78++){
7116 helicities[7]=allowed_helicities[h78];
7117 helicities[8]=allowed_helicities[h78];
7118 for (
unsigned int ihf=0; ihf<Hffparticles.size(); ihf++){
7119 vh_ids[7]=Hffparticles.at(ihf).first;
7120 vh_ids[8]=Hffparticles.at(ihf).second;
7121 double msq_inst_LR=0;
7123 msq_inst += msq_inst_LR;
7126 msq_inst *= Hffscale;
7138 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\tScale for outgoing particles: " << scale_out << endl;
7141 for (
int incoming1=1; incoming1<=
nf; incoming1++){
7142 if (partonIsKnown[0] && incoming1!=abs(vh_ids[0]))
continue;
7143 if (incoming1%2!=abs(vh_ids[0])%2 || incoming1==6)
continue;
7144 int iin = incoming1 * TMath::Sign(1, vh_ids[0]);
7146 if (verbosity>=
TVar::DEBUG)
MELAout <<
"\t\t\tiin,jin = " << iin <<
"," << jin << endl;
7147 MatElsq[jin+5][iin+5] +=
msq * 0.25 * scale_out;
7154 int GeVexponent_MEsq;
7155 if (HDKon==0) GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets+nRequested_AssociatedLeptons+nRequested_AssociatedPhotons)*2-4;
7156 else GeVexponent_MEsq = 4-(2+nRequested_AssociatedJets+nRequested_AssociatedLeptons+nRequested_AssociatedPhotons)*2;
7157 double constant = pow(GeV, -GeVexponent_MEsq);
7158 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7161 for (
int ii = 0; ii <
nmsq; ii++){
for (
int jj = 0; jj <
nmsq; jj++)
MELAout << MatElsq[jj][ii] <<
'\t';
MELAout << endl; }
7163 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7173 <<
"TUtil::VHiggsMatEl: Reset AlphaS:\n"
7174 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7176 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7180 <<
"TUtil::VHiggsMatEl: Reset AlphaS result:\n"
7181 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7190 const double& EBEAM,
7191 int topDecay,
int topProcess,
7194 const double GeV=1./100.;
7195 double sum_msqjk = 0;
7196 double MatElsq[
nmsq][
nmsq]={ { 0 } };
7198 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::TTHiggsMatEl: Non-JHUGen MEs are not supported." << endl;
return sum_msqjk; }
7202 int nRequested_Tops=1;
7203 int nRequested_Antitops=1;
7220 <<
"TUtil::TTHiggsMatEl: Number of stable tops (" << mela_event.
pStableTops.size() <<
")"
7221 <<
" and number of stable antitops (" << mela_event.
pStableAntitops.size() <<
")"
7222 <<
" in ttH process are not 1!" << endl;
7227 <<
"TUtil::TTHiggsMatEl: Number of set of top daughters (" << mela_event.
pTopDaughters.size() <<
")"
7228 <<
" and number of set of antitop daughters (" << mela_event.
pAntitopDaughters.size() <<
")"
7229 <<
" in ttH process are not 1!" << endl;
7234 <<
"TUtil::TTHiggsMatEl: Number of top daughters (" << mela_event.
pTopDaughters.at(0).size() <<
")"
7235 <<
" and number of antitop daughters (" << mela_event.
pAntitopDaughters.at(0).size() <<
")"
7236 <<
" in ttH process are not 3!" << endl;
7242 bool isUnknown[2]={
true,
true };
7246 for (
unsigned int itd=0; itd<mela_event.
pTopDaughters.at(0).size(); itd++) topDaughters.push_back(mela_event.
pTopDaughters.at(0).at(itd));
7250 for (
unsigned int itop=0; itop<mela_event.
pStableTops.size(); itop++) topDaughters.push_back(mela_event.
pStableTops.at(itop));
7254 for (
unsigned int itd=0; itd<topDaughters.size(); itd++){
if (!
PDGHelpers::isAnUnknownJet(topDaughters.at(itd).first)){ isUnknown[0]=
false;
break; } }
7255 for (
unsigned int itd=0; itd<antitopDaughters.size(); itd++){
if (!
PDGHelpers::isAnUnknownJet(antitopDaughters.at(itd).first)){ isUnknown[1]=
false;
break; } }
7262 double p4[13][4]={ { 0 } };
7263 double MYIDUP_prod[2]={ 0 };
7264 TLorentzVector MomStore[
mxpart];
7265 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
7266 for (
int ipar=0; ipar<2; ipar++){
7267 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
7268 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
7270 else MYIDUP_prod[ipar] = 0;
7271 if (momTmp->T()>0.){
7272 p4[ipar][0] = -momTmp->T()*GeV;
7273 p4[ipar][1] = -momTmp->X()*GeV;
7274 p4[ipar][2] = -momTmp->Y()*GeV;
7275 p4[ipar][3] = -momTmp->Z()*GeV;
7276 MomStore[ipar] = (*momTmp);
7279 p4[ipar][0] = momTmp->T()*GeV;
7280 p4[ipar][1] = momTmp->X()*GeV;
7281 p4[ipar][2] = momTmp->Y()*GeV;
7282 p4[ipar][3] = momTmp->Z()*GeV;
7283 MomStore[ipar] = -(*momTmp);
7284 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
7290 const unsigned int t_pos=4;
7291 const unsigned int b_pos=9;
7292 const unsigned int Wp_pos=10;
7293 const unsigned int Wpf_pos=12;
7294 const unsigned int Wpfb_pos=11;
7295 for (
unsigned int ipar=0; ipar<topDaughters.size(); ipar++){
7296 TLorentzVector* momTmp = &(topDaughters.at(ipar).second);
7297 if (topDaughters.size()==1){
7298 p4[t_pos][0] = momTmp->T()*GeV;
7299 p4[t_pos][1] = momTmp->X()*GeV;
7300 p4[t_pos][2] = momTmp->Y()*GeV;
7301 p4[t_pos][3] = momTmp->Z()*GeV;
7305 p4[b_pos][0] = momTmp->T()*GeV;
7306 p4[b_pos][1] = momTmp->X()*GeV;
7307 p4[b_pos][2] = momTmp->Y()*GeV;
7308 p4[b_pos][3] = momTmp->Z()*GeV;
7311 p4[Wpfb_pos][0] = momTmp->T()*GeV;
7312 p4[Wpfb_pos][1] = momTmp->X()*GeV;
7313 p4[Wpfb_pos][2] = momTmp->Y()*GeV;
7314 p4[Wpfb_pos][3] = momTmp->Z()*GeV;
7315 p4[Wp_pos][0] += p4[Wpfb_pos][0];
7316 p4[Wp_pos][1] += p4[Wpfb_pos][1];
7317 p4[Wp_pos][2] += p4[Wpfb_pos][2];
7318 p4[Wp_pos][3] += p4[Wpfb_pos][3];
7321 p4[Wpf_pos][0] = momTmp->T()*GeV;
7322 p4[Wpf_pos][1] = momTmp->X()*GeV;
7323 p4[Wpf_pos][2] = momTmp->Y()*GeV;
7324 p4[Wpf_pos][3] = momTmp->Z()*GeV;
7325 p4[Wp_pos][0] += p4[Wpf_pos][0];
7326 p4[Wp_pos][1] += p4[Wpf_pos][1];
7327 p4[Wp_pos][2] += p4[Wpf_pos][2];
7328 p4[Wp_pos][3] += p4[Wpf_pos][3];
7330 MomStore[6] = MomStore[6] + (*momTmp);
7332 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]; } }
7336 const unsigned int tb_pos=3;
7337 const unsigned int bb_pos=5;
7338 const unsigned int Wm_pos=6;
7339 const unsigned int Wmf_pos=7;
7340 const unsigned int Wmfb_pos=8;
7341 for (
unsigned int ipar=0; ipar<antitopDaughters.size(); ipar++){
7342 TLorentzVector* momTmp = &(antitopDaughters.at(ipar).second);
7343 if (antitopDaughters.size()==1){
7344 p4[tb_pos][0] = momTmp->T()*GeV;
7345 p4[tb_pos][1] = momTmp->X()*GeV;
7346 p4[tb_pos][2] = momTmp->Y()*GeV;
7347 p4[tb_pos][3] = momTmp->Z()*GeV;
7351 p4[bb_pos][0] = momTmp->T()*GeV;
7352 p4[bb_pos][1] = momTmp->X()*GeV;
7353 p4[bb_pos][2] = momTmp->Y()*GeV;
7354 p4[bb_pos][3] = momTmp->Z()*GeV;
7357 p4[Wmf_pos][0] = momTmp->T()*GeV;
7358 p4[Wmf_pos][1] = momTmp->X()*GeV;
7359 p4[Wmf_pos][2] = momTmp->Y()*GeV;
7360 p4[Wmf_pos][3] = momTmp->Z()*GeV;
7361 p4[Wm_pos][0] += p4[Wmf_pos][0];
7362 p4[Wm_pos][1] += p4[Wmf_pos][1];
7363 p4[Wm_pos][2] += p4[Wmf_pos][2];
7364 p4[Wm_pos][3] += p4[Wmf_pos][3];
7367 p4[Wmfb_pos][0] = momTmp->T()*GeV;
7368 p4[Wmfb_pos][1] = momTmp->X()*GeV;
7369 p4[Wmfb_pos][2] = momTmp->Y()*GeV;
7370 p4[Wmfb_pos][3] = momTmp->Z()*GeV;
7371 p4[Wm_pos][0] += p4[Wmfb_pos][0];
7372 p4[Wm_pos][1] += p4[Wmfb_pos][1];
7373 p4[Wm_pos][2] += p4[Wmfb_pos][2];
7374 p4[Wm_pos][3] += p4[Wmfb_pos][3];
7376 MomStore[7] = MomStore[7] + (*momTmp);
7378 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]; } }
7380 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
7381 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
7382 p4[2][0] += momTmp->T()*GeV;
7383 p4[2][1] += momTmp->X()*GeV;
7384 p4[2][2] += momTmp->Y()*GeV;
7385 p4[2][3] += momTmp->Z()*GeV;
7386 MomStore[5] = MomStore[5] + (*momTmp);
7390 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; }
7393 double defaultRenScale =
scale_.scale;
7394 double defaultFacScale =
facscale_.facscale;
7396 int defaultNflav =
nflav_.nflav;
7397 string defaultPdflabel =
pdlabel_.pdlabel;
7401 double alphasVal, alphasmzVal;
7411 <<
"TUtil::TTHiggsMatEl: Set AlphaS:\n"
7412 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
7414 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7419 for (
unsigned int ib1=0; ib1<3; ib1++){
7420 if (topDaughters.size()==1 && ib1!=0)
continue;
7422 unsigned int b1index;
7423 if (ib1==0) b1index = b_pos;
7424 else if (ib1==1) b1index = Wpf_pos;
7425 else b1index = Wpfb_pos;
7426 for (
unsigned int if1=0; if1<3; if1++){
7427 if (topDaughters.size()==1 && if1!=0)
continue;
7429 unsigned int f1index;
7430 if (if1==0) f1index = Wpf_pos;
7431 else if (if1==2) f1index = b_pos;
7432 else f1index = Wpfb_pos;
7433 for (
unsigned int ifb1=0; ifb1<3; ifb1++){
7434 if (topDaughters.size()==1 && ifb1!=0)
continue;
7436 unsigned int fb1index;
7437 if (ifb1==2) fb1index = Wpf_pos;
7438 else if (ifb1==1) fb1index = b_pos;
7439 else fb1index = Wpfb_pos;
7441 if (b1index==f1index || b1index==fb1index || f1index==fb1index)
continue;
7443 for (
unsigned int ib2=0; ib2<3; ib2++){
7444 if (antitopDaughters.size()==1 && ib2!=0)
continue;
7446 unsigned int b2index;
7447 if (ib2==0) b2index = bb_pos;
7448 else if (ib2==1) b2index = Wmf_pos;
7449 else b2index = Wmfb_pos;
7450 for (
unsigned int if2=0; if2<3; if2++){
7451 if (antitopDaughters.size()==1 && if2!=0)
continue;
7453 unsigned int f2index;
7454 if (if2==0) f2index = Wmf_pos;
7455 else if (if2==2) f2index = bb_pos;
7456 else f2index = Wmfb_pos;
7457 for (
unsigned int ifb2=0; ifb2<3; ifb2++){
7458 if (antitopDaughters.size()==1 && ifb2!=0)
continue;
7460 unsigned int fb2index;
7461 if (ifb2==2) fb2index = Wmf_pos;
7462 else if (ifb2==1) fb2index = bb_pos;
7463 else fb2index = Wmfb_pos;
7465 if (b2index==f2index || b2index==fb2index || f2index==fb2index)
continue;
7467 double p4_current[13][4]={ { 0 } };
7468 for (
unsigned int ix=0; ix<4; ix++){
7469 for (
unsigned int ip=0; ip<=2; ip++) p4_current[ip][ix] = p4[ip][ix];
7470 p4_current[t_pos][ix] = p4[t_pos][ix];
7471 p4_current[b1index][ix] = p4[b_pos][ix];
7472 p4_current[f1index][ix] = p4[Wpf_pos][ix];
7473 p4_current[fb1index][ix] = p4[Wpfb_pos][ix];
7474 p4_current[Wp_pos][ix] = p4_current[Wpf_pos][ix] + p4_current[Wpfb_pos][ix];
7476 p4_current[tb_pos][ix] = p4[tb_pos][ix];
7477 p4_current[b2index][ix] = p4[bb_pos][ix];
7478 p4_current[f2index][ix] = p4[Wmf_pos][ix];
7479 p4_current[fb2index][ix] = p4[Wmfb_pos][ix];
7480 p4_current[Wm_pos][ix] = p4_current[Wmf_pos][ix] + p4_current[Wmfb_pos][ix];
7484 <<
"TUtil::TTHiggsMatEl: Unswapped instance for "
7485 <<
"b(" << b_pos <<
") -> " << b1index <<
", "
7486 <<
"Wpf(" << Wpf_pos <<
") -> " << f1index <<
", "
7487 <<
"Wpfb(" << Wpfb_pos <<
") -> " << fb1index <<
", "
7488 <<
"bb(" << bb_pos <<
") -> " << b2index <<
", "
7489 <<
"Wmf(" << Wmf_pos <<
") -> " << f2index <<
", "
7490 <<
"Wmfb(" << Wmfb_pos <<
") -> " << fb2index << endl;
7491 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; }
7495 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
7496 double MatElsq_tmp_swap[
nmsq][
nmsq]={ { 0 } };
7498 if (isUnknown[0] && isUnknown[1]){
7499 for (
unsigned int ix=0; ix<4; ix++){
7500 swap(p4_current[t_pos][ix], p4_current[tb_pos][ix]);
7501 swap(p4_current[b_pos][ix], p4_current[bb_pos][ix]);
7502 swap(p4_current[Wp_pos][ix], p4_current[Wm_pos][ix]);
7503 swap(p4_current[Wpf_pos][ix], p4_current[Wmf_pos][ix]);
7504 swap(p4_current[Wpfb_pos][ix], p4_current[Wmfb_pos][ix]);
7507 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.; }
7509 for (
int ix=0; ix<11; ix++){
for (
int iy=0; iy<11; iy++) MatElsq[iy][ix] += MatElsq_tmp[iy][ix]; }
7512 <<
"TUtil::TTHiggsMatEl: Swapped instance for "
7513 <<
"b(" << b_pos <<
") -> " << b1index <<
", "
7514 <<
"Wpf(" << Wpf_pos <<
") -> " << f1index <<
", "
7515 <<
"Wpfb(" << Wpfb_pos <<
") -> " << fb1index <<
", "
7516 <<
"bb(" << bb_pos <<
") -> " << b2index <<
", "
7517 <<
"Wmf(" << Wmf_pos <<
") -> " << f2index <<
", "
7518 <<
"Wmfb(" << Wmfb_pos <<
") -> " << fb2index << endl;
7519 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; }
7529 int defaultTopDecay=-1;
7532 int GeVexponent_MEsq;
7533 if (topDecay>0) GeVexponent_MEsq = 4-(1+3*(nRequested_Tops+nRequested_Antitops))*2;
7534 else GeVexponent_MEsq = 4-(1+nRequested_Tops+nRequested_Antitops)*2;
7535 double constant = pow(GeV, -GeVexponent_MEsq);
7536 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7538 MELAout <<
"TUtil::TTHiggsMatEl: MEsq[ip][jp] = " << endl;
7539 for (
int iquark=-5; iquark<=5; iquark++){
7540 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << MatElsq[jquark+5][iquark+5] <<
'\t';
7544 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7548 <<
"TUtil::TTHiggsMatEl: Reset AlphaS:\n"
7549 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7551 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7555 <<
"TUtil::TTHiggsMatEl: Reset AlphaS result:\n"
7556 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7565 const double& EBEAM,
7569 const double GeV=1./100.;
7570 double sum_msqjk = 0;
7571 double MatElsq[
nmsq][
nmsq]={ { 0 } };
7572 double MatElsq_tmp[
nmsq][
nmsq]={ { 0 } };
7574 if (matrixElement!=
TVar::JHUGen){
if (verbosity>=
TVar::ERROR)
MELAerr <<
"TUtil::BBHiggsMatEl: Non-JHUGen MEs are not supported." << endl;
return sum_msqjk; }
7578 int nRequested_AssociatedJets=2;
7591 <<
"TUtil::BBHiggsMatEl: Number of stable bs (" << mela_event.
pAssociated.size() <<
")"
7592 <<
" in bbH process is not 2!" << endl;
7598 bool isUnknown[2]; isUnknown[0]=
false; isUnknown[1]=
false;
7601 else antitopDaughters.push_back(mela_event.
pAssociated.at(0));
7603 else topDaughters.push_back(mela_event.
pAssociated.at(1));
7610 double p4[13][4]={ { 0 } };
7611 double MYIDUP_prod[2]={ 0 };
7612 TLorentzVector MomStore[
mxpart];
7613 for (
int i = 0; i <
mxpart; i++) MomStore[i].SetXYZT(0, 0, 0, 0);
7614 for (
int ipar=0; ipar<2; ipar++){
7615 TLorentzVector* momTmp = &(mela_event.
pMothers.at(ipar).second);
7616 int* idtmp = &(mela_event.
pMothers.at(ipar).first);
7618 else MYIDUP_prod[ipar] = 0;
7619 if (momTmp->T()>0.){
7620 p4[ipar][0] = -momTmp->T()*GeV;
7621 p4[ipar][1] = -momTmp->X()*GeV;
7622 p4[ipar][2] = -momTmp->Y()*GeV;
7623 p4[ipar][3] = -momTmp->Z()*GeV;
7624 MomStore[ipar] = (*momTmp);
7627 p4[ipar][0] = momTmp->T()*GeV;
7628 p4[ipar][1] = momTmp->X()*GeV;
7629 p4[ipar][2] = momTmp->Y()*GeV;
7630 p4[ipar][3] = momTmp->Z()*GeV;
7631 MomStore[ipar] = -(*momTmp);
7632 MYIDUP_prod[ipar] = -MYIDUP_prod[ipar];
7637 for (
unsigned int ipar=0; ipar<topDaughters.size(); ipar++){
7638 TLorentzVector* momTmp = &(topDaughters.at(ipar).second);
7639 p4[4][0] = momTmp->T()*GeV;
7640 p4[4][1] = momTmp->X()*GeV;
7641 p4[4][2] = momTmp->Y()*GeV;
7642 p4[4][3] = momTmp->Z()*GeV;
7643 MomStore[6] = MomStore[6] + (*momTmp);
7647 for (
unsigned int ipar=0; ipar<antitopDaughters.size(); ipar++){
7648 TLorentzVector* momTmp = &(antitopDaughters.at(ipar).second);
7649 p4[3][0] = momTmp->T()*GeV;
7650 p4[3][1] = momTmp->X()*GeV;
7651 p4[3][2] = momTmp->Y()*GeV;
7652 p4[3][3] = momTmp->Z()*GeV;
7653 MomStore[7] = MomStore[7] + (*momTmp);
7656 for (
unsigned int ipar=0; ipar<mela_event.
pDaughters.size(); ipar++){
7657 TLorentzVector* momTmp = &(mela_event.
pDaughters.at(ipar).second);
7658 p4[2][0] += momTmp->T()*GeV;
7659 p4[2][1] += momTmp->X()*GeV;
7660 p4[2][2] += momTmp->Y()*GeV;
7661 p4[2][3] += momTmp->Z()*GeV;
7662 MomStore[5] = MomStore[5] + (*momTmp);
7666 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; }
7669 double defaultRenScale =
scale_.scale;
7670 double defaultFacScale =
facscale_.facscale;
7673 int defaultNflav =
nflav_.nflav;
7674 string defaultPdflabel =
pdlabel_.pdlabel;
7678 double alphasVal, alphasmzVal;
7688 <<
"TUtil::BBHiggsMatEl: Set AlphaS:\n"
7689 <<
"\tBefore set, alphas scale: " << defaultRenScale <<
", PDF scale: " << defaultFacScale <<
'\n'
7691 <<
"\tAfter set, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7695 if (isUnknown[0] && isUnknown[1]){
7696 for (
unsigned int ix=0; ix<4; ix++) swap(p4[3][ix], p4[4][ix]);
7698 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.; }
7701 int GeVexponent_MEsq = 4-(1+nRequested_AssociatedJets)*2;
7702 double constant = pow(GeV, -GeVexponent_MEsq);
7703 for (
int ii=0; ii<
nmsq; ii++){
for (
int jj=0; jj<
nmsq; jj++) MatElsq[jj][ii] *= constant; }
7705 MELAout <<
"TUtil::BBHiggsMatEl: MEsq[ip][jp] = " << endl;
7706 for (
int iquark=-5; iquark<=5; iquark++){
7707 for (
int jquark=-5; jquark<=5; jquark++)
MELAout << MatElsq[jquark+5][iquark+5] <<
'\t';
7711 sum_msqjk =
SumMEPDF(MomStore[0], MomStore[1], MatElsq, RcdME, EBEAM, verbosity);
7715 <<
"TUtil::BBHiggsMatEl: Reset AlphaS:\n"
7716 <<
"\tBefore reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale << endl;
7718 SetAlphaS(defaultRenScale, defaultFacScale, 1., 1., defaultNloop, defaultNflav, defaultPdflabel);
7722 <<
"TUtil::BBHiggsMatEl: Reset AlphaS result:\n"
7723 <<
"\tAfter reset, alphas scale: " <<
scale_.scale <<
", PDF scale: " <<
facscale_.facscale <<
", alphas(Qren): " << alphasVal <<
", alphas(MZ): " << alphasmzVal << endl;
7731 MELAout <<
"TUtil::WipeMEArray: Initial MEsq:" << endl;
7732 for (
int iquark=-5; iquark<=5; iquark++){
7733 for (
int jquark=-5; jquark<=5; jquark++)
MELAout <<
msq[jquark+5][iquark+5] <<
'\t';
7745 for (
int iquark=-5; iquark<=5; iquark++){
7746 for (
int jquark=-5; jquark<=5; jquark++){
7752 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
7754 else msq[jquark+5][iquark+5]=0;
7764 if (
msq[5][5]>0.) nInstances=1;
7782 for (
int iquark=-5; iquark<=5; iquark++){
7783 for (
int jquark=-5; jquark<=5; jquark++){
7789 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
7791 else msq[jquark+5][iquark+5]=0;
7797 for (
int iquark=-5; iquark<=5; iquark++){
7798 for (
int jquark=-5; jquark<=5; jquark++){
7818 (TMath::Sign(1, iquark)==TMath::Sign(1, jquark))
7835 )
msq[jquark+5][iquark+5]=0;
7838 int order[2]={ -1, -1 };
7840 if (order[0]<0 || order[1]<0)
msq[jquark+5][iquark+5]=0;
7842 if (
msq[jquark+5][iquark+5]>0.) nInstances++;
7844 else msq[jquark+5][iquark+5]=0;
7858 xx[0]>1. || xx[0]<=
xmin_.xmin
7860 xx[1]>1. || xx[1]<=
xmin_.xmin
7865 if (xx[0]>1. || xx[1]>1.)
MELAerr <<
"TUtil::CheckPartonMomFraction: At least one of the parton momentum fractions is greater than 1." << endl;
7866 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;
7867 else MELAerr <<
"TUtil::CheckPartonMomFraction: EBEAM=" << EBEAM <<
"<=0." << endl;
7872 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::CheckPartonMomFraction: xx[0]: " << xx[0] <<
", xx[1] = " << xx[1] <<
", xmin = " <<
xmin_.xmin << endl;
7882 double fx1x2_jhu[2][13]={ { 0 } };
7883 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::ComputePDF: Calling setpdfs with xx[0]: " << xx[0] <<
", xx[1] = " << xx[1] << endl;
7886 for (
int ip=-6; ip<=6; ip++){
7889 if (ip!=0 && (abs(ip)%2==0)) fac=-1;
7890 else if (ip!=0) fac=+1;
7893 fx1[jp+5]=fx1x2_jhu[0][ip+6];
7894 fx2[jp+5]=fx1x2_jhu[1][ip+6];
7906 MELAout <<
"End TUtil::ComputePDF:"<< endl;
7907 for (
int ip=-
nf; ip<=
nf; ip++)
MELAout <<
"(fx1, fx2)[" << ip <<
"] = (" << fx1[ip+5] <<
" , " << fx2[ip+5] <<
")" << endl;
7913 double fx1[
nmsq]={ 0 };
7914 double fx2[
nmsq]={ 0 };
7917 ComputePDF(p0, p1, fx1, fx2, EBEAM, verbosity);
7932 const double GeV=1./100.;
7933 int isch=(
int)scheme;
7934 double shat_jhu = pow(
sqrts*GeV, 2);
7956 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: code=" << code << endl;
7958 TLorentzVector
const nullFourVector(0, 0, 0, 0);
7961 vector<int> idVstar; idVstar.reserve(2);
7965 idVstar.push_back(melaCand->
id);
7966 idVstar.push_back(-9000);
7973 for (
unsigned char iv=0; iv<2; iv++){
7976 int idtmp = Vdau->
id;
7982 idVstar.push_back(idtmp);
7985 else idVstar.push_back(-9000);
7990 for (
size_t iv=0; iv<nffs; iv++){
7995 daughters.at(2*iv+0).second = corrPair.first;
7996 daughters.at(2*iv+1).second = corrPair.second;
7999 TLorentzVector tmp = nullFourVector;
8002 lastDau.second, lastDau.first,
8005 lastDau.second = corrPair.first;
8015 int nsatisfied_jets=0;
8016 int nsatisfied_lnus=0;
8017 int nsatisfied_gammas=0;
8018 vector<MELAParticle*> candidateVs;
8021 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo!=0 case start" << endl;
8028 if ((abs(idV)==aVhypo || idV==0) && Vdau->getNDaughters()>0 && Vdau->passSelection){
8031 if (!Vdau_i){ doAdd=
false;
break; }
8037 !Vdau_i->passSelection
8043 if (doAdd) candidateVs.push_back(Vdau);
8047 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: candidateVs size = " << candidateVs.size() << endl;
8075 if (associated_tmp.size()>=2 || (associated_tmp.size()==1 &&
PDGHelpers::isAPhoton(associated_tmp.at(0).first))) idVstar.push_back(Vdau->id);
8077 if (associated_tmp.size()>=2){
8078 unsigned int nffs = associated_tmp.size()/2;
8079 for (
unsigned int iv=0; iv<nffs; iv++){
8081 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8082 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8084 associated_tmp.at(2*iv+0).second = corrPair.first;
8085 associated_tmp.at(2*iv+1).second = corrPair.second;
8087 if (2*nffs<associated_tmp.size()){
8088 TLorentzVector tmp = nullFourVector;
8090 associated_tmp.at(associated_tmp.size()-1).second, associated_tmp.at(associated_tmp.size()-1).first,
8093 associated_tmp.at(associated_tmp.size()-1).second = corrPair.first;
8099 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo!=0 case associated.size=" <<
associated.size() << endl;
8103 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 case begin" << endl;
8114 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 lep case associated_tmp.size=" << associated_tmp.size() << endl;
8117 if (associated_tmp.size()>=1){
8118 size_t nffs = associated_tmp.size()/2;
8119 for (
size_t iv=0; iv<nffs; iv++){
8120 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from lepton pair " << 2*iv+0 <<
'\t' << 2*iv+1 << endl;
8122 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8123 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8125 associated_tmp.at(2*iv+0).second = corrPair.first;
8126 associated_tmp.at(2*iv+1).second = corrPair.second;
8128 if (2*nffs<associated_tmp.size()){
8129 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from last lepton " << associated_tmp.size()-1 << endl;
8130 TLorentzVector tmp = nullFourVector;
8131 SimpleParticle_t& lastAssociated = associated_tmp.at(associated_tmp.size()-1);
8133 lastAssociated.second, lastAssociated.first,
8136 lastAssociated.second = corrPair.first;
8151 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: aVhypo==0 jet case associated_tmp.size=" << associated_tmp.size() << endl;
8154 if (associated_tmp.size()>=1){
8155 unsigned int nffs = associated_tmp.size()/2;
8156 for (
unsigned int iv=0; iv<nffs; iv++){
8157 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from jet pair " << 2*iv+0 <<
'\t' << 2*iv+1 << endl;
8159 associated_tmp.at(2*iv+0).second, associated_tmp.at(2*iv+0).first,
8160 associated_tmp.at(2*iv+1).second, associated_tmp.at(2*iv+1).first
8162 associated_tmp.at(2*iv+0).second = corrPair.first;
8163 associated_tmp.at(2*iv+1).second = corrPair.second;
8165 if (2*nffs<associated_tmp.size()){
8166 if (verbosity>=
TVar::DEBUG)
MELAout <<
"TUtil::GetBoostedParticleVectors: Removing mass from last jet " << associated_tmp.size()-1 << endl;
8167 TLorentzVector tmp = nullFourVector;
8168 SimpleParticle_t& lastAssociated = associated_tmp.at(associated_tmp.size()-1);
8170 lastAssociated.second, lastAssociated.first,
8173 lastAssociated.second = corrPair.first;
8184 nsatisfied_gammas++;
8194 int nsatisfied_tops=0;
8195 int nsatisfied_antitops=0;
8196 vector<SimpleParticleCollection_t> topDaughters;
8197 vector<SimpleParticleCollection_t> antitopDaughters;
8201 vector<MELATopCandidate_t*> tops;
8202 vector<MELATopCandidate_t*> topbars;
8203 vector<MELATopCandidate_t*> unknowntops;
8208 if (theTop && theTop->passSelection){
8209 vector<MELATopCandidate_t*>* particleArray;
8210 if (theTop->id==6) particleArray = &tops;
8211 else if (theTop->id==-6) particleArray = &topbars;
8212 else particleArray = &unknowntops;
8217 ) particleArray->push_back(theTop);
8220 if (verbosity>=
TVar::DEBUG){
MELAout <<
"TUtil::GetBoostedParticleVectors: tops.size=" << tops.size() <<
", topbars.size=" << topbars.size() <<
", unknowntops.size=" << unknowntops.size() << endl; }
8240 if (vdaughters.size()==3) topDaughters.push_back(vdaughters);
8245 nsatisfied_antitops++;
8249 nsatisfied_antitops++;
8260 if (vdaughters.size()==3) antitopDaughters.push_back(vdaughters);
8273 nsatisfied_antitops++;
8289 if (vdaughters.size()==3) topDaughters.push_back(vdaughters);
8292 nsatisfied_antitops++;
8303 if (vdaughters.size()==3) antitopDaughters.push_back(vdaughters);
8310 MELAout <<
"TUtil::GetBoostedParticleVectors: stableTops.size=" << stableTops.size() << endl;
8311 MELAout <<
"TUtil::GetBoostedParticleVectors: stableAntitops.size=" << stableAntitops.size() << endl;
8312 MELAout <<
"TUtil::GetBoostedParticleVectors: topDaughters.size=" << topDaughters.size() << endl;
8313 for (
unsigned int itop=0; itop<topDaughters.size(); itop++)
MELAout <<
"TUtil::GetBoostedParticleVectors: topDaughters.at(" << itop <<
").size=" << topDaughters.at(itop).size() << endl;
8314 MELAout <<
"TUtil::GetBoostedParticleVectors: antitopDaughters.size=" << antitopDaughters.size() << endl;
8315 for (
unsigned int itop=0; itop<antitopDaughters.size(); itop++)
MELAout <<
"TUtil::GetBoostedParticleVectors: antitopDaughters.at(" << itop <<
").size=" << antitopDaughters.at(itop).size() << endl;
8320 TLorentzVector pTotal(0, 0, 0, 0);
8324 for (
SimpleParticle_t const& sp:stableAntitops) pTotal = pTotal + sp.second;
8329 double qX = pTotal.X();
8330 double qY = pTotal.Y();
8331 double qE = pTotal.T();
8332 if ((qX*qX+qY*qY)>0.){
8333 TVector3 boostV(-qX/qE, -qY/qE, 0.);
8340 pTotal.Boost(boostV);
8348 double sysPz= pTotal.Z();
8349 double sysE = pTotal.T();
8350 double pz0 = (sysE+sysPz)/2.;
8351 double pz1 = -(sysE-sysPz)/2.;
8354 int motherId[2]={ 0, 0 };
8356 for (
int ip=0; ip<2; ip++) motherId[ip]=melaCand->
getMother(ip)->
id;
8359 if (TMath::Sign(1., melaCand->
getMother(0)->
z()-melaCand->
getMother(1)->
z())!=TMath::Sign(1., pz0-pz1)){ swap(pz0, pz1); swap(E0, E1); }
8361 if ((motherId[0]<0 && motherId[1]>=0) || (motherId[1]>0 && motherId[0]<=0)){
8364 swap(motherId[0], motherId[1]);
8367 TLorentzVector pM[2];
8368 pM[0].SetXYZT(0., 0., pz0, E0);
8369 pM[1].SetXYZT(0., 0., pz1, E1);
8396 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.intermediateVid.size=" << mela_event.
intermediateVid.size() << endl;
8397 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pMothers.size=" << mela_event.
pMothers.size() << endl;
8398 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pDaughters.size=" << mela_event.
pDaughters.size() << endl;
8399 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pAssociated.size=" << mela_event.
pAssociated.size() << endl;
8400 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pStableTops.size=" << mela_event.
pStableTops.size() << endl;
8401 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pStableAntitops.size=" << mela_event.
pStableAntitops.size() << endl;
8402 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pTopDaughters.size=" << mela_event.
pTopDaughters.size() << endl;
8403 MELAout <<
"TUtil::GetBoostedParticleVectors mela_event.pAntitopDaughters.size=" << mela_event.
pAntitopDaughters.size() << endl;
8404 MELAout <<
"End GetBoostedParticleVectors" << endl;
8417 std::vector<MELAParticle*>* particleList,
8418 std::vector<MELACandidate*>* candList
8422 if (!pDaughters){
MELAerr <<
"TUtil::ConvertVectorFormat: No daughters!" << endl;
return cand; }
8423 else if (pDaughters->size()==0){
MELAerr <<
"TUtil::ConvertVectorFormat: Daughter size==0!" << endl;
return cand; }
8424 else if (pDaughters->size()>4){
MELAerr <<
"TUtil::ConvertVectorFormat: Daughter size " << pDaughters->size() <<
">4 is not supported!" << endl;
return cand; }
8425 if (pMothers && pMothers->size()!=2){
MELAerr <<
"TUtil::ConvertVectorFormat: Mothers momentum size (" << pMothers->size() <<
") has to have had been 2! Continuing by omitting mothers." << endl; }
8429 std::vector<MELAParticle*> aparticles;
8430 std::vector<MELAParticle*>
mothers;
8431 for (
auto& spart:(*pDaughters)){
8434 if (particleList) particleList->push_back(onePart);
8438 for (
auto& spart:(*pAssociated)){
8441 if (particleList) particleList->push_back(onePart);
8442 aparticles.push_back(onePart);
8445 if (pMothers && pMothers->size()==2){
8446 for (
auto& spart:(*pMothers)){
8449 if (particleList) particleList->push_back(onePart);
8471 TLorentzVector pH = F1->
p4+F2->
p4;
8496 TLorentzVector pH = F1->
p4+F2->
p4+gamma->
p4;
8510 TLorentzVector pH(0, 0, 0, 0);
8512 for (
unsigned char ip=0; ip<4; ip++){ pH = pH + (
daughters.at(ip))->p4; charge += (
daughters.at(ip))->charge(); }
8523 MELAerr <<
"TUtil::ConvertVectorFormat: PDGHelpers::HDecayMode = " <<
PDGHelpers::HDecayMode <<
" is set incorrectly for Ndaughters=" <<
daughters.size() <<
". There is no automatic mechnism for this scenario. ";
8524 MELAerr <<
"Suggestion: Call PDGHelpers::setCandidateDecayMode(TVar::CandidateDecay_XY) before this call." << endl;
8537 if (!aparticles.empty()){
8538 for (
auto& part:aparticles){
8539 const int& partId = part->id;
8548 if (candList && cand) candList->push_back(cand);
8558 std::vector<MELAParticle*>* particleList,
8559 std::vector<MELAThreeBodyDecayCandidate*>* tbdCandList
8563 if (!tbdDaughters){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: No daughters!" << endl;
return cand; }
8564 else if (tbdDaughters->empty()){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: Daughter size==0!" << endl;
return cand; }
8565 else if (!(tbdDaughters->size()==1 || tbdDaughters->size()==3)){
MELAerr <<
"TUtil::ConvertThreeBodyDecayCandidate: Daughter size " << tbdDaughters->size() <<
"!=1 or 3 is not supported!" << endl;
return cand; }
8567 if (tbdDaughters->size()==1){
8568 if (abs((tbdDaughters->at(0)).first)==6 || (tbdDaughters->at(0)).first==0){
8570 tbdCandList->push_back(cand);
8573 else if (tbdDaughters->size()==3){
8578 if (Wf->
id<0 || Wfb->
id>0) std::swap(Wf, Wfb);
8580 particleList->push_back(partnerPart);
8581 particleList->push_back(Wf);
8582 particleList->push_back(Wfb);
8585 tbdCandList->push_back(cand);
8591 MELAout <<
"***** TUtil::PrintCandidateSummary *****" << endl;
8592 MELAout <<
"Candidate: " << cand << endl;
8597 MELAout <<
"***** TUtil::PrintCandidateSummary (Simple Event Record) *****" << endl;
8598 MELAout <<
"Candidate: " << cand << endl;
8608 for (
unsigned int ip=0; ip<cand->
pMothers.size(); ip++){
8611 <<
"\t\tV" << ip <<
" (" << part->first <<
") (X,Y,Z,T)=( "
8612 << part->second.X() <<
" , "
8613 << part->second.Y() <<
" , "
8614 << part->second.Z() <<
" , "
8615 << part->second.T() <<
" )" << endl;
8620 for (
unsigned int ip=0; ip<cand->
pDaughters.size(); ip++){
8623 <<
"\t\tDau[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8624 << part->second.X() <<
" , "
8625 << part->second.Y() <<
" , "
8626 << part->second.Z() <<
" , "
8627 << part->second.T() <<
" )" << endl;
8630 for (
unsigned int ip=0; ip<cand->
pAssociated.size(); ip++){
8633 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8634 << part->second.X() <<
" , "
8635 << part->second.Y() <<
" , "
8636 << part->second.Z() <<
" , "
8637 << part->second.T() <<
" )" << endl;
8640 for (
unsigned int ip=0; ip<cand->
pStableTops.size(); ip++){
8643 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8644 << part->second.X() <<
" , "
8645 << part->second.Y() <<
" , "
8646 << part->second.Z() <<
" , "
8647 << part->second.T() <<
" )" << endl;
8653 <<
"\t\tAPart[" << ip <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8654 << part->second.X() <<
" , "
8655 << part->second.Y() <<
" , "
8656 << part->second.Z() <<
" , "
8657 << part->second.T() <<
" )" << endl;
8661 for (
unsigned int ip=0; ip<cand->
pTopDaughters.size(); ip++){
8662 MELAout <<
"\t\tTop[" << ip <<
"] daughters:" << endl;
8663 for (
unsigned int jp=0; jp<cand->
pTopDaughters.at(ip).size(); jp++){
8666 <<
"\t\t- Top daughter[" << ip << jp <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8667 << part->second.X() <<
" , "
8668 << part->second.Y() <<
" , "
8669 << part->second.Z() <<
" , "
8670 << part->second.T() <<
" )" << endl;
8675 MELAout <<
"\t\tAntitop[" << ip <<
"] daughters:" << endl;
8679 <<
"\t\t- Antitop daughter[" << ip << jp <<
"] (" << part->first <<
") (X,Y,Z,T)=( "
8680 << part->second.X() <<
" , "
8681 << part->second.Y() <<
" , "
8682 << part->second.Z() <<
" , "
8683 << part->second.T() <<
" )" << endl;