#include "TauSpinner/tau_reweight_lib.h" #include "TauSpinner/Particle.h" #include "TauSpinner/vbftests.h" #include #include using namespace Tauolapp; /*This is a group of methods prepared for the testing purpose. The tests include those where fixed kinematics is forced. Then external calculation can be used for matrix element calculation and/or PDFs etc. Gradually better documentation will be introduced. It will be documented with comment lines and in the draft of the paper. To use one has to modify the Makefile the appropriate line of objects has to be activated (# removed). Then invoking lines in code have to be commented out as well. At a certain moment such commented out lines will be removed from the code or will remain. We will have to make up the mind. -*/ namespace TauSpinner { extern double CMSENE; extern int relWTnonSM; extern double WTnonSM; extern double Polari; extern double WTamplit; extern double WTamplitP; extern double WTamplitM; extern double CMSENE; extern double f(double x, int ID, double SS, double cmsene); const bool DEBUG = 0; /** Wrapper to VBDISTR */ double vbfdistrWRAP(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY) { double P_copy[6][4]; memcpy(P_copy,P,sizeof(double)*6*4); return vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY); } //--------------------------------------------------------------------------------------------------------------------------------------------------- void calcXsect(int IDPROD, SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY) { /* // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation! // FSR photons may need to be treated explicitely or with interpolation procedures. Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(), p3.py()+p4.py()+tau1.py()+tau2.py(), p3.pz()+p4.pz()+tau1.pz()+tau2.pz(), p3.e() +p4.e() +tau1.e() +tau2.e(), 0 ); */ Particle P_QQ( p3.px()+p4.px()+sp_X.px(), p3.py()+p4.py()+sp_X.py(), p3.pz()+p4.pz()+sp_X.pz(), p3.e() +p4.e() +sp_X.e() , 0 ); double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass(); double x1x2 = SS/CMSENE/CMSENE; double x1Mx2 = P_QQ.pz()/CMSENE*2; double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; //--------------------------------------------------------------------------- // Construct the matrix for FORTRAN function // NOTE: different order of indices than in FORTRAN! double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 }, { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 }, { p3.e(), p3.px(), p3.py(), p3.pz() }, { p4.e(), p4.px(), p4.py(), p4.pz() }, { tau1.e(), tau1.px(), tau1.py(), tau1.pz() }, { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } }; W[0][0]=0.0; W[0][1]=0.0; W[1][0]=0.0; W[1][1]=0.0; // // Calculate 'f' function for all x1 and all ID1, ID2 // double f_x1_ARRAY[11] = { 0.0 }; double f_x2_ARRAY[11] = { 0.0 }; double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5 double SSfix = 91.188*91.188; double CMSENEfix = 91.188; for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,SSfix,CMSENEfix); f_x2[i] = f(x2,i,SSfix,CMSENEfix); // std::cout << "ERW: x, Q2, f(x,Q2), f(x,Q)=" << x1 << " " << SS << " " << f(x1,i,SS,CMSENE) << " " << f(x1,i,sqrt(SS),CMSENE) << std::endl; } W[0][0]=0.0; W[0][1]=0.0; W[1][0]=0.0; W[1][1]=0.0; // these loops need to be cleaned from zero contributions! for(int I1=-4;I1<=4;I1++){ // for test of single production process fix flavour for(int I2=-4;I2<=4;I2++){ // for test of single production process fix flavour for(int I3=-4;I3<=4;I3++){ for(int I4=-4;I4<=4;I4++){ int ID1 = I1; if( ID1 == 0 ) ID1 = 21; int ID2 = I2; if( ID2 == 0 ) ID2 = 21; int ID3 = I3; if( ID3 == 0 ) ID3 = 21; int ID4 = I4; if( ID4 == 0 ) ID4 = 21; //preselect production group bool accept = false; // any process if( IDPROD == 0) accept = true; // only GG processes if( IDPROD == 1 && ID1 == 21 && ID2 == 21 && abs (ID3) < 5 && abs (ID4) < 5 ) accept = true; // only GG processes ... details if( IDPROD == 110 && ID1 == 21 && ID2 == 21 ) { accept = true; if( abs (ID3) != 2 || abs (ID4) != 2 ) accept = false; } if( IDPROD == 111 && ID1 == 21 && ID2 == 21 ) { accept = true; if( abs (ID3) != 1 || abs (ID4) != 1 ) accept = false; } if( IDPROD == 112 && ID1 == 21 && ID2 == 21 ) { accept = true; if( abs (ID3) != 3 || abs (ID4) != 3 ) accept = false; } if( IDPROD == 113 && ID1 == 21 && ID2 == 21 ) { accept = true; if( abs (ID3) != 4 || abs (ID4) != 4 ) accept = false; } // only GQ processes if( IDPROD == 2 && ID1 == 21 && abs(ID2) < 5 ) accept = true ; if( IDPROD == 2 && ID2 == 21 && abs(ID1) < 5 ) accept = true ; // only GQ processes ... details if( (IDPROD == 210) && ( ID1 == 21) && ( ID2 < 5 ) && (ID2 > 0) ) accept = true; if( (IDPROD == 211) && ( ID1 == 21) && ( ID2 < 5 ) && (ID2 < 0) ) accept = true; if( (IDPROD == 212) && ( ID2 == 21) && ( ID1 < 5 ) && (ID1 > 0) ) accept = true; if( (IDPROD == 213) && ( ID2 == 21) && ( ID1 < 5 ) && (ID1 < 0) ) accept = true; // only QQ, QXQX processes if( IDPROD == 3 && (ID1 * ID2 > 0 ) && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true; // only QQ processes .. details if( IDPROD == 31 && ID1 > 0 && ID2 > 0 && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true; if( IDPROD == 32 && ID1 < 0 && ID2 < 0 && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true; // only QQ processes .. details if( IDPROD == 310 && ID1 > 0 && ID2 > 0 && ID1 < 5 && ID2 < 5 && ID1 == ID2 ) accept = true; if( IDPROD == 311 && ID1 > 0 && ID2 > 0 && ID1 < 5 && ID2 < 5 && ID1 != ID2 ) accept = true; if( IDPROD == 312 && ID1 < 0 && ID2 < 0 && ID1 == ID2 ) accept = true; if( IDPROD == 313 && ID1 < 0 && ID2 < 0 && ID1 != ID2 ) accept = true; // only QQ processes .. details if( IDPROD == 320 && ID1 * ID2 > 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) == abs(ID2) ) accept = true; if( IDPROD == 321 && ID1 * ID2 > 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) != abs(ID2) ) accept = true; // only QQ processes .. details if( IDPROD == 3101 && ID1 == 2 && ID2 == 2 && ID3 == 2 && ID4 == 2 ) accept = true; if( IDPROD == 3102 && ID1 == -2 && ID2 == -2 && ID3 == -2 && ID4 == -2 ) accept = true; if( IDPROD == 3103 && ID1 == 2 && ID2 == 2 ) accept = true; if( IDPROD == 3104 && ID1 == -2 && ID2 == -2 ) accept = true; // only QQX processes if( IDPROD == 4 && (ID1 * ID2 < 0 ) && abs (ID1) < 5 && abs (ID2) < 5 ) accept = true; // only QQX processes .. details if( IDPROD == 420 && ID1 * ID2 < 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) == abs(ID2) ) accept = true; if( IDPROD == 421 && ID1 * ID2 < 0 && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID1) != abs(ID2) ) accept = true; // only QQX processes, first family in initial/final if( IDPROD == 41 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) { accept = true; if( abs (ID1) > 2 || abs (ID2) > 2 || abs (ID3) >2 || abs (ID4) > 2 ) accept = false; } // only QQX processes, second family in initial/final state if( IDPROD == 42 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) { accept = true; if( abs (ID1) < 3 || abs (ID2) < 3 || abs (ID3) < 3 || abs (ID4) < 3 ) accept = false; } // only QQX processes, mixed families in initial/final state if( IDPROD == 43 && ( ID1 * ID2 < 0) && abs(ID1) < 5 && abs(ID2) < 5 && abs(ID3) < 5 && abs(ID4) < 5 ) { accept = true; if( abs (ID1) < 3 && abs (ID2) < 3 && abs (ID3) < 3 && abs (ID4) < 3 ) accept = false; if( abs (ID1) > 2 && abs (ID2) > 2 && abs (ID3) > 2 && abs (ID4) > 2 ) accept = false; } // only QQX processes, first family in initial state if( IDPROD == 44 && (ID1 * ID2 < 0 ) && abs(ID1) < 5 && abs(ID2) < 5 ) { accept = true; if( abs (ID1) > 2 || abs (ID2) > 2 ) accept = false; } // only QQX processes ... details if( IDPROD == 4101 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> U UX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false; if( ID1 != 2 ) accept = false; } if( IDPROD == 4102 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> U UX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false; if( ID1 != -2 ) accept = false; } if( IDPROD == 4103&& (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> Q QX if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false; if( ID1 != 2 ) accept = false; } if( IDPROD == 4104&& (ID1 * ID2 < 0 ) ) { accept = true; // only UX U --> Q QX if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false; if( ID1 != -2 ) accept = false; } if( IDPROD == 4105&& (ID1 * ID2 < 0 ) ) { accept = true; // only D DX --> Q QX if( abs (ID1) != 1 || abs (ID2) != 1 ) accept = false; if( ID1 != 1 ) accept = false; } if( IDPROD == 4106&& (ID1 * ID2 < 0 ) ) { accept = true; // only DX D --> Q QX if( abs (ID1) != 1 || abs (ID2) != 1 ) accept = false; if( ID1 != -1 ) accept = false; } if( IDPROD == 4107&& (ID1 * ID2 < 0 ) ) { accept = true; // only Q QX --> Q QX if( ID1 < 0 || ID2 > 0) accept = false; } if( IDPROD == 4108&& (ID1 * ID2 < 0 ) ) { accept = true; // only QX Q --> Q QX if( ID1 > 0 || ID2 < 0) accept = false; } // only QQX processes ... details if( IDPROD == 410 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> U UX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 2 || abs (ID4) != 2 ) accept = false; } if( IDPROD == 411 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> C CX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 4 || abs (ID4) != 4 ) accept = false; } if( IDPROD == 412 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> S DX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 3 || abs (ID4) != 1 ) accept = false; } if( IDPROD == 413 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> D SX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 1 || abs (ID4) != 3 ) accept = false; } if( IDPROD == 414 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> D DX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 1 || abs (ID4) != 1 ) accept = false; } if( IDPROD == 415 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> S SX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 3 || abs (ID4) != 3 ) accept = false; } if( IDPROD == 416 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> Q QX if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 417 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> G G if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) == 21 || abs (ID4) == 21 ) accept = false; } if( IDPROD == 418 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> J J if( abs (ID1) != 2 || abs (ID2) != 2 ) accept = false; } if( IDPROD == 419 && (ID1 * ID2 < 0 ) ) { accept = true; // only D DX --> D DX if( abs (ID1) != 1 || abs (ID2) != 1 || abs (ID3) != 1 || abs (ID4) != 1 ) accept = false; } if( IDPROD == 511 && (ID1 * ID2 < 0 ) ) { accept = true; // only D DX, DX D --> G G if( abs (ID1) != 1 || abs (ID2) != 1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 512 && (ID1 * ID2 < 0 ) ) { accept = true; // only D DX --> G G if( ID1 != 1 || ID2 != -1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 513 && (ID1 * ID2 < 0 ) ) { accept = true; // only DX D --> G G if( ID1 != -1 || ID2 != 1 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 521 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX, UX U --> G G if( abs (ID1) != 2 || abs (ID2) != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 522 && (ID1 * ID2 < 0 ) ) { accept = true; // only U UX --> G G if( ID1 != 2 || ID2 != -2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 523 && (ID1 * ID2 < 0 ) ) { accept = true; // only UX U --> G G if( ID1 != -2 || ID2 != 2 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 531 && (ID1 * ID2 < 0 ) ) { accept = true; // only Q QX, QX Q --> G G if( abs(ID1) > 4 || abs(ID2) > 4 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 532 && (ID1 * ID2 < 0 ) ) { accept = true; // only Q QX --> G G if( abs(ID1) > 4 || abs(ID2) > 4 || ID1 < 0 || abs(ID3) != 21 || abs(ID4) != 21 ) accept = false; } if( IDPROD == 533 && (ID1 * ID2 < 0 ) ) { accept = true; // only QX Q --> G G if( abs(ID1) > 4 || abs(ID2) > 4 || ID2 < 0 || abs (ID3) != 21 || abs (ID4) != 21 ) accept = false; } if( IDPROD == 541 ) { accept = false; // only U SX --> D CX if( ID1 == 2 && ID2 == -3 && ID3 == 1 && ID4 == -4 ) accept = true; if( ID1 == -3 && ID2 == 2 && ID3 == 1 && ID4 == -4 ) accept = true; if( ID1 == 2 && ID2 == -3 && ID3 ==-4 && ID4 == 1 ) accept = true; if( ID1 == -3 && ID2 == 2 && ID3 ==-4 && ID4 == 1 ) accept = true; } if(accept ){ W[0][0] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY); // as in case of nonSM_adopt we change the sign for W[0][1] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY); // second tau helicity assuming without checking that W[1][0] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY); // for VBF quantization frame conventions are the same. W[1][1] += f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY); } } } } } } void makeSimpleTestME2() { //--------------------------------------------------------------------------- // Construct the matrix for FORTRAN function // NOTE: different order of indices than in FORTRAN! std::cout << "ERW: ----------------------------- " << std::endl; double P1[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 }, { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 }, { 0.8855009E+02, -0.2210038E+02, 0.4007979E+02, -0.7580437E+02 }, { 0.3283248E+03, -0.1038482E+03, -0.3019295E+03, 0.7649385E+02 }, { 0.1523663E+03, -0.1058795E+03, -0.9770827E+02, 0.4954769E+02 }, { 0.4307588E+03, 0.2318280E+03, 0.3595580E+03, -0.5023717E+02 } }; calcTestME2(1, P1); Particle p1 ( 0.0, 0.0, 0.5000000E+03, 0.5000000E+03, 1 ); Particle p2 ( 0.0, 0.0, -0.5000000E+03, 0.5000000E+03, -1 ); Particle p3 ( -0.2210038E+02, 0.4007979E+02, -0.7580437E+02, 0.8855009E+02, 2 ); Particle p4 ( -0.1038482E+03, -0.3019295E+03, 0.7649385E+02, 0.3283248E+03, -2 ); Particle tau1 ( -0.1058795E+03, -0.9770827E+02, 0.4954769E+02, 0.1523663E+03, 15 ); Particle tau2 ( 0.2318280E+03, 0.3595580E+03, -0.5023717E+02, 0.4307588E+03, -15 ); std::cout << "ERW: 4-momenta before boost ----------------------------- " << std::endl; p1.print(); p2.print(); p3.print(); p4.print(); tau1.print(); tau2.print(); double x1=0.03; // x1 for the first parton double xm = 1000.0; double cmsene = 13000.0; double x2=xm/cmsene*xm/cmsene/x1; // x2 adjusted from cms energy (13000), and virtuality of the whole system 1000 std::cout << "x1, x2= " << x1 << " " << x2 << std::endl; if (x2>=1.0) std::cout << "ERW: wrong x1, x2 cannot be defined (is above 1); x2=" << x2 << std::endl; Particle P_Z0( tau1.px()+tau2.px(), tau1.py()+tau2.py(), tau1.pz()+tau2.pz(), tau1.e()+tau2.e(), 23 ); P_Z0.print(); double Q2 = P_Z0.recalculated_mass() * P_Z0.recalculated_mass(); std::cout << " pdfs = " << f(x1,0,Q2,cmsene) * f(x2,1,Q2,cmsene) << std::endl; Particle P_QQ( 0.0, 1.0e-10, (x1-x2)*cmsene/2.0, (x1+x2)*cmsene/2.0, 10 ); P_QQ.print(); p1.boostFromRestFrame(P_QQ); p2.boostFromRestFrame(P_QQ); p3.boostFromRestFrame(P_QQ); p4.boostFromRestFrame(P_QQ); P_Z0.boostFromRestFrame(P_QQ); tau1.boostFromRestFrame(P_QQ); tau2.boostFromRestFrame(P_QQ); std::cout << "ERW: 4-momenta after boost ----------------------------- " << std::endl; p1.print(); p2.print(); p3.print(); p4.print(); P_Z0.print(); tau1.print(); tau2.print(); std::cout << "ERW: ----------------------------- " << std::endl; double P2[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 }, { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 }, { 0.1177462E+03, -0.6070512E+02, 0.7123011E+02, 0.7145150E+02 }, { 0.3509495E+03, -0.3178775E+02, 0.8393832E+02, 0.3392779E+03 }, { 0.3493321E+03, 0.1840069E+03, -0.5152712E+02, -0.2924315E+03 }, { 0.1819722E+03, -0.9151401E+02, -0.1036413E+03, -0.1182978E+03 } }; calcTestME2(2, P2); std::cout << "ERW: ----------------------------- " << std::endl; double P3[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 }, { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 }, { 0.2586900E+03, 0.1324670E+03, -0.1696171E+03, -0.1435378E+03 }, { 0.1084567E+03, -0.5735712E+02, -0.2162482E+02, -0.8947281E+02 }, { 0.4005742E+03, -0.1580760E+03, 0.3563160E+03, 0.9223569E+02 }, { 0.2322791E+03, 0.8296613E+02, -0.1650741E+03, 0.1407749E+03 } }; calcTestME2(3, P3); std::cout << "ERW: ----------------------------- " << std::endl; double P4[6][4] = {{ 0.5000000E+03, 0.0, 0.0, 0.5000000E+03 }, { 0.5000000E+03, 0.0, 0.0, -0.5000000E+03 }, { 0.1595700E+03, -0.6917808E+02, -0.1395175E+03, -0.3481123E+02 }, { 0.2247758E+03, -0.1360140E+03, 0.1650340E+03, -0.6919641E+02 }, { 0.2508802E+03, 0.1447863E+01, 0.2499830E+03, -0.2107335E+02 }, { 0.3647740E+03, 0.2037442E+03, -0.2754995E+03, 0.1250810E+03 } }; calcTestME2(4, P4); } /*******************************************************************************/ void calcTestME2(int iter, double P1[6][4]) { double ME2ref, ME2; int ID1, ID2, ID3, ID4; int KEY=0; // case GD -> GD ID1 = 21; ID2 = 1; ID3 = 21; ID4 = 1; if(iter == 1) ME2ref = 1.5961579933867344E-010; if(iter == 2) ME2ref = 3.8749742050329834E-010; if(iter == 3) ME2ref = 5.0434636937545397E-012; if(iter == 4) ME2ref = 7.9077184257060427E-012; ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY); std::cout << "ERW: ME2 GD->GD (ref) = " << ME2ref << std::endl; std::cout << "ERW: ME2 GD->GD = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl; // case GU -> GU ID1 = 21; ID2 = 2; ID3 = 21; ID4 = 2; if(iter == 1) ME2ref = 2.9195503763051040E-010; ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY); if(iter == 1)std::cout << "ERW: ME2 GU->GU (ref) = " << ME2ref << std::endl; if(iter == 1)std::cout << "ERW: ME2 GU->GU = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl; // case DD -> DD ID1 = 1; ID2 = 1; ID3 = 1; ID4 = 1; if(iter == 1) ME2ref = 3.3953129762581284E-017; if(iter == 2) ME2ref = 6.0054072781075002E-017; if(iter == 3) ME2ref = 3.9707833415682912E-018; if(iter == 4) ME2ref = 1.5342177192807347E-018; ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY); std::cout << "ERW: ME2 DD->DD (ref) = " << ME2ref << std::endl; std::cout << "ERW: ME2 DD->DD = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl; // case UU -> UU ID1 = 2; ID2 = 2; ID3 = 2; ID4 = 2; if(iter == 1) ME2ref = 1.9412924824120248E-017; if(iter == 2) ME2ref = 4.0830132078559964E-017; if(iter == 3) ME2ref = 2.1297931556738857E-018; if(iter == 4) ME2ref = 7.9215386777281075E-019; ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY); std::cout << "ERW: ME2 UU->UU (ref) = " << ME2ref << std::endl; std::cout << "ERW: ME2 UU->UU = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl; // case DDX -> CCX ID1 = 1; ID2 = -1; ID3 = 4; ID4 = -4; if(iter == 1) ME2ref = 3.6780119739685137E-020; if(iter == 2) ME2ref = 5.2280923900274694E-018; if(iter == 3) ME2ref = 2.9001589049209953E-019; if(iter == 4) ME2ref = 5.8509026904432882E-020; ME2 = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, 1, P1, KEY) + vbfdistrWRAP(ID1,ID2,ID3,ID4, -1, -1, P1, KEY); std::cout << "ERW: ME2 DDX->CCX (ref) = " << ME2ref << std::endl; std::cout << "ERW: ME2 DDX->CCX = " << ME2 << " ratio to ref = " << ME2/ME2ref << std::endl; } // // ------------------------------------------------------------------------------------------------------- // ------------------------------------------------------------------------------------------------------- // ------------------------------------------------------------------------------------------------------- // void makeSimpleTestPDF() { std::cout << "------------------------------------" << std::endl; std::cout << "ERW: TEST PDF only = " << std::endl; std::cout << "ERW: PDF:LHAPDFset = cteq6ll.LHpdf " << std::endl; std::cout << "ERW: reference from PYTHIA printout " << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 =2 -2 " << std::endl; std::cout << "ERW: Q,Q2 = 97.868, 9578.118 " << std::endl; double pdf1= f(2.51297e-01,2,9578.118,97.868); std::cout << "ERW: x1, xpdf1= " << 2.51297e-01 << " " << 2.51297e-01 *pdf1 << " ref pdf = " << 0.411 << std::endl; double pdf2= f(1.94463e-04,-2,9578.118,97.868); std::cout << "ERW: x2, xpdf2= " << 1.94463e-04 << " " << 1.94463e-04 *pdf2 << " ref pdf = " << 2.989 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 =2 -2 " << std::endl; std::cout << "ERW: Q,Q2 = 38.747, 1501.314 " << std::endl; pdf1= f(4.17758e-04,-2,1501.314,38.747); std::cout << "ERW: x1, xpdf1 = " << 4.17758e-04 << " " << 4.17758e-04 *pdf1 << " ref pdf = " << 1.768 << std::endl; pdf2= f(1.83354e-02, 2,1501.314,38.747); std::cout << "ERW: x2, xpdf2 = " << 1.83354e-02 << " " << 1.83354e-02 *pdf2 << " ref pdf = " << 0.610 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 =1 -1 " << std::endl; std::cout << "ERW: Q,Q2 = 91.585, 8387.831 " << std::endl; pdf1= f(6.29049e-02, 1,8387.831,91.585); std::cout << "ERW: x1, xpdf1 = " << 6.29049e-02 << " " << 6.29049e-02 *pdf1 << " ref pdf = " << 0.386 << std::endl; pdf2= f(6.80314e-04, -1,8387.831,91.585); std::cout << "ERW: x2, xpdf2 = " << 6.80314e-04 << " " << 6.80314e-04 *pdf2 << " ref pdf = " << 1.751 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 =-2 2 " << std::endl; std::cout << "ERW: Q,Q2 = 12.979, 168.453 " << std::endl; pdf1= f( 1.29739e-01, -2,168.453,12.979); std::cout << "ERW: x1, xpdf1 = " << 1.29739e-01 << " " << 1.29739e-01*pdf1 << " ref pdf = " << 0.066 << std::endl; pdf2= f( 6.62449e-06, 2,168.453,12.979); std::cout << "ERW: x2, xpdf2 = " << 6.62449e-06 << " " << 6.62449e-06 *pdf2 << " ref pdf = " << 4.559 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 = 2 -2 " << std::endl; std::cout << "ERW: Q,Q2 = 91.707, 8410.127 " << std::endl; pdf1= f(1.73131e-02 , 2,8410.127,91.707); std::cout << "ERW: x1, xpdf1 = " << 1.73131e-02 << " " << 1.73131e-02 *pdf1 << " ref pdf = " << 0.643 << std::endl; pdf2= f( 2.47840e-03, -2,8410.127,91.707); std::cout << "ERW: x2, xpdf2 = " << 2.47840e-03 << " " << 2.47840e-03 *pdf2 << " ref pdf = " << 0.986 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 = 1 -1 " << std::endl; std::cout << "ERW: Q,Q2 = 91.203 8317.936 " << std::endl; pdf1= f( 3.87999e-04 , 1, 8317.936, 91.203); std::cout << "ERW: x1, xpdf1 = " << 3.87999e-04 << " " << 3.87999e-04 *pdf1 << " ref pdf = " << 2.256 << std::endl; pdf2= f( 1.09378e-01, -1,8317.936, 91.203); std::cout << "ERW: x2, xpdf2 = " << 1.09378e-01 << " " << 1.09378e-01 *pdf2 << " ref pdf = " << 0.106 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 = 1 0 " << std::endl; std::cout << "ERW: Q,Q2 = 51.963, 2700.201 " << std::endl; pdf1= f( 9.22157e-02 , 1, 2700.201, 51.963); std::cout << "ERW: x1, xpdf1 = " << 9.22157e-02 << " " << 9.22157e-02 *pdf1 << " ref pdf = " << 0.353 << std::endl; pdf2= f( 2.50786e-03 , 0,2700.201, 51.963); std::cout << "ERW: x2, xpdf2 = " << 2.50786e-03 << " " << 2.50786e-03 *pdf2 << " ref pdf = " << 20.362 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 = 1 0 " << std::endl; std::cout << "ERW: Q,Q2 = 55.719, 3104.648 " << std::endl; pdf1= f( 1.28045e-01 , 1, 3104.648, 55.719); std::cout << "ERW: x1, xpdf1 = " << 1.28045e-01 << " " << 1.28045e-01 *pdf1 << " ref pdf = " << 0.312 << std::endl; pdf2= f( 3.12973e-03 , 0, 3104.648, 55.719); std::cout << "ERW: x2, xpdf2 = " << 3.12973e-03 << " " << 3.12973e-03 *pdf2 << " ref pdf = " << 17.938 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: ID1,ID2 = 3 0 " << std::endl; std::cout << "ERW: Q,Q2 = 106.126, 11262.651 " << std::endl; pdf1= f( 8.02832e-02 , 3, 11262.651, 106.126); std::cout << "ERW: x1, xpdf1 = " << 8.02832e-02 << " " << 8.02832e-02 *pdf1 << " ref pdf = " << 0.077 << std::endl; pdf2= f( 3.85011e-03 , 0, 11262.651, 106.126); std::cout << "ERW: x2, xpdf2 = " << 3.85011e-03 << " " << 3.85011e-03 *pdf2 << " ref pdf = " << 16.542 << std::endl; std::cout << "------------------------------------" << std::endl; std::cout << "ERW: Q= 98.069051344553060, x =1.36851248162658170E-002 " << std::endl; double x = 1.36851248162658170E-002; std::cout << "ERW: f(-5) = " << f( 1.36851248162658170E-002 , -5, 9617.53, 98.069 ) << " ref=12.649157008859417" << std::endl; std::cout << "ERW: f(-4) = " << f( 1.36851248162658170E-002 , -4, 9617.53, 98.069 ) << " ref=19.196982264477899" << std::endl; std::cout << "ERW: f(-3) = " << f( 1.36851248162658170E-002 , -3, 9617.53, 98.069 ) << " ref=23.971647011421041" << std::endl; std::cout << "ERW: f(-2) = " << f( 1.36851248162658170E-002 , -2, 9617.53, 98.069 ) << " ref=30.602439987620933" << std::endl; std::cout << "ERW: f(-1) = " << f( 1.36851248162658170E-002 , -1, 9617.53, 98.069 ) << " ref=31.664848276050574" << std::endl; std::cout << "ERW: f(0) = " << f( 1.36851248162658170E-002 , 0, 9617.53, 98.069 ) << " ref=483.13004227606962" << std::endl; std::cout << "ERW: f( 1) = " << f( 1.36851248162658170E-002 , 1, 9617.53, 98.069 ) << " ref=41.868651977841559" << std::endl; std::cout << "ERW: f( 2) = " << f( 1.36851248162658170E-002 , 2, 9617.53, 98.069 ) << " ref=49.256024133250129" << std::endl; std::cout << "ERW: f( 3) = " << f( 1.36851248162658170E-002 , 3, 9617.53, 98.069 ) << " ref=23.971647011421041" << std::endl; std::cout << "ERW: f( 4) = " << f( 1.36851248162658170E-002 , 4, 9617.53, 98.069 ) << " ref=19.196982264477899" << std::endl; std::cout << "ERW: f( 5) = " << f( 1.36851248162658170E-002 , 5, 9617.53, 98.069 ) << " ref=12.649157008859417" << std::endl; std::cout << "ERW: f(-5) = " << f( 1.36851248162658170E-002 , -5, 98.069, 98.069 ) << " ref=12.649157008859417" << std::endl; std::cout << "ERW: f(-4) = " << f( 1.36851248162658170E-002 , -4, 98.069, 98.069 ) << " ref=19.196982264477899" << std::endl; std::cout << "ERW: f(-3) = " << f( 1.36851248162658170E-002 , -3, 98.069, 98.069 ) << " ref=23.971647011421041" << std::endl; std::cout << "ERW: f(-2) = " << f( 1.36851248162658170E-002 , -2, 98.069, 98.069 ) << " ref=30.602439987620933" << std::endl; std::cout << "ERW: f(-1) = " << f( 1.36851248162658170E-002 , -1, 98.069, 98.069 ) << " ref=31.664848276050574" << std::endl; std::cout << "ERW: f(0) = " << f( 1.36851248162658170E-002 , 0, 98.069, 98.069 ) << " ref=483.13004227606962" << std::endl; std::cout << "ERW: f( 1) = " << f( 1.36851248162658170E-002 , 1, 98.069, 98.069 ) << " ref=41.868651977841559" << std::endl; std::cout << "ERW: f( 2) = " << f( 1.36851248162658170E-002 , 2, 98.069, 98.069 ) << " ref=49.256024133250129" << std::endl; std::cout << "ERW: f( 3) = " << f( 1.36851248162658170E-002 , 3, 98.069, 98.069 ) << " ref=23.971647011421041" << std::endl; std::cout << "ERW: f( 4) = " << f( 1.36851248162658170E-002 , 4, 98.069, 98.069 ) << " ref=19.196982264477899" << std::endl; std::cout << "ERW: f( 5) = " << f( 1.36851248162658170E-002 , 5, 98.069, 98.069 ) << " ref=12.649157008859417" << std::endl; } //--------------------------------------------------------------------------------------------------------------------------------------------------- void calcSumME2(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4) { /* // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation! // FSR photons may need to be treated explicitely or with interpolation procedures. Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(), p3.py()+p4.py()+tau1.py()+tau2.py(), p3.pz()+p4.pz()+tau1.pz()+tau2.pz(), p3.e() +p4.e() +tau1.e() +tau2.e(), 0 ); */ Particle P_QQ( p3.px()+p4.px()+sp_X.px(), p3.py()+p4.py()+sp_X.py(), p3.pz()+p4.pz()+sp_X.pz(), p3.e() +p4.e() +sp_X.e() , 0 ); double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass(); double x1x2 = SS/CMSENE/CMSENE; double x1Mx2 = P_QQ.pz()/CMSENE*2; double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; //--------------------------------------------------------------------------- // Construct the matrix for FORTRAN function // NOTE: different order of indices than in FORTRAN! double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 }, { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 }, { p3.e(), p3.px(), p3.py(), p3.pz() }, { p4.e(), p4.px(), p4.py(), p4.pz() }, { tau1.e(), tau1.px(), tau1.py(), tau1.pz() }, { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } }; W[0][0] = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY); W[0][1] = vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY); W[1][0] = vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY); W[1][1] = vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY); double sumW =(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); if( sumW < 0 ){ std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: W[] = " << W[0][0] << " " << W[0][1] << " " << W[1][0] << " " << W[1][1] << std::endl; std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: p1 = " << CMSENE/2*x1 << " " << 0.0 << " " << 0.0 << " " << CMSENE/2*x1 << std::endl; std::cout << "ERW: p2 = " << CMSENE/2*x2 << " " << 0.0 << " " << 0.0 << " " << -CMSENE/2*x2 << std::endl; std::cout << "ERW: X = " << sp_X.e() << " " << sp_X.px() << " " << sp_X.py() << " " << sp_X.pz() << std::endl; std::cout << "ERW: p3 = " << p3.e() << " " << p3.px() << " " << p3.py() << " " << p3.pz() << std::endl; std::cout << "ERW: p4 = " << p4.e() << " " << p4.px() << " " << p4.py() << " " << p4.pz() << std::endl; std::cout << "ERW: tau1 = " << tau1.e() << " " << tau1.px() << " " << tau1.py() << " " << tau1.pz() << std::endl; std::cout << "ERW: tau2 = " << tau2.e() << " " << tau2.px() << " " << tau2.py() << " " << tau2.pz() << std::endl; std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: p1+p2= " << P[0][0]+P[1][0] << " " << P[0][1]+P[1][1] << " " << P[0][2]+P[1][2] << " " << P[0][3]+P[1][3] << std::endl; std::cout << "ERW: X+p3+p4= " << P[2][0]+P[3][0]+P[4][0]+P[5][0] << " " << P[2][1]+P[3][1]+P[4][1]+P[5][1] << " " << P[2][2]+P[3][2]+P[4][2]+P[5][2] << " " << P[2][3]+P[3][3]+P[4][3]+P[5][3] << std::endl; double p3sqr = p3.e()*p3.e()-p3.px()*p3.px()-p3.py()*p3.py()-p3.pz()*p3.pz(); double p4sqr = p4.e()*p4.e()-p4.px()*p4.px()-p4.py()*p4.py()-p4.pz()*p4.pz(); std::cout << "ERW: p3^2, p4^2 =" << p3sqr << " " << p4sqr << std::endl; } } //--------------------------------------------------------------------------------------------------------------------------------------------------- void calcPDFs(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4, int pdfOpt) { /* // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation! // FSR photons may need to be treated explicitely or with interpolation procedures. Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(), p3.py()+p4.py()+tau1.py()+tau2.py(), p3.pz()+p4.pz()+tau1.pz()+tau2.pz(), p3.e() +p4.e() +tau1.e() +tau2.e(), 0 ); */ Particle P_QQ( p3.px()+p4.px()+sp_X.px(), p3.py()+p4.py()+sp_X.py(), p3.pz()+p4.pz()+sp_X.pz(), p3.e() +p4.e() +sp_X.e() , 0 ); double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass(); double x1x2 = SS/CMSENE/CMSENE; double x1Mx2 = P_QQ.pz()/CMSENE*2; double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; //--------------------------------------------------------------------------- // Construct the matrix for FORTRAN function // NOTE: different order of indices than in FORTRAN! double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 }, { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 }, { p3.e(), p3.px(), p3.py(), p3.pz() }, { p4.e(), p4.px(), p4.py(), p4.pz() }, { tau1.e(), tau1.px(), tau1.py(), tau1.pz() }, { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } }; // // Calculate 'f' function for all x1 and all ID1, ID2 // double f_x1_ARRAY[11] = { 0.0 }; double f_x2_ARRAY[11] = { 0.0 }; double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5 Particle P_X( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e() , 0 ); double Q2 = P_X.recalculated_mass()*P_X.recalculated_mass(); double PT2 = sp_X.px()* sp_X.px() + sp_X.py()* sp_X.py(); if(pdfOpt == 0){ for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,Q2,CMSENE); f_x2[i] = f(x2,i,Q2,CMSENE); } } else if (pdfOpt == 1) { for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,PT2,CMSENE); f_x2[i] = f(x2,i,PT2,CMSENE); } } else if (pdfOpt == 2) { double PT24 = PT2/4.; for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,PT24,CMSENE); f_x2[i] = f(x2,i,PT24,CMSENE); } } int I1=ID1; int I2=ID2; if( I1 == 21) I1=0; if( I2 == 21) I2=0; W[0][0] = f_x1[I1]*f_x2[I2]; W[0][1] = f_x1[I1]*f_x2[I2]; W[1][0] = f_x1[I1]*f_x2[I2]; W[1][1] = f_x1[I1]*f_x2[I2]; } //--------------------------------------------------------------------------------------------------------------------------------------------------- void calcProdMatrix(SimpleParticle &p3, SimpleParticle &p4, SimpleParticle &sp_X,SimpleParticle &tau1, SimpleParticle &tau2, double (&W)[2][2], int KEY, int ID1, int ID2, int ID3, int ID4, int pdfOpt) { /* // this may be necessary because of matrix element calculation may require absolute energy-momentum conservation! // FSR photons may need to be treated explicitely or with interpolation procedures. Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px(), p3.py()+p4.py()+tau1.py()+tau2.py(), p3.pz()+p4.pz()+tau1.pz()+tau2.pz(), p3.e() +p4.e() +tau1.e() +tau2.e(), 0 ); */ Particle P_QQ( p3.px()+p4.px()+sp_X.px(), p3.py()+p4.py()+sp_X.py(), p3.pz()+p4.pz()+sp_X.pz(), p3.e() +p4.e() +sp_X.e() , 0 ); double SS = P_QQ.recalculated_mass()*P_QQ.recalculated_mass(); double x1x2 = SS/CMSENE/CMSENE; double x1Mx2 = P_QQ.pz()/CMSENE*2; double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2; //--------------------------------------------------------------------------- // Construct the matrix for FORTRAN function // NOTE: different order of indices than in FORTRAN! double P[6][4] = { { CMSENE/2*x1, 0.0, 0.0, CMSENE/2*x1 }, { CMSENE/2*x2, 0.0, 0.0, -CMSENE/2*x2 }, { p3.e(), p3.px(), p3.py(), p3.pz() }, { p4.e(), p4.px(), p4.py(), p4.pz() }, { tau1.e(), tau1.px(), tau1.py(), tau1.pz() }, { tau2.e(), tau2.px(), tau2.py(), tau2.pz() } }; // // Calculate 'f' function for all x1 and all ID1, ID2 // double f_x1_ARRAY[11] = { 0.0 }; double f_x2_ARRAY[11] = { 0.0 }; double *f_x1 = f_x1_ARRAY+5; // This way f_x1[i],f_x2[i] can be used with double *f_x2 = f_x2_ARRAY+5; // i going from -5 to 5 Particle P_X( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e() , 0 ); double Q2 = P_X.recalculated_mass()*P_X.recalculated_mass(); double PT2 = sp_X.px()* sp_X.px() + sp_X.py()* sp_X.py(); if(pdfOpt == 0){ for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,Q2,CMSENE); f_x2[i] = f(x2,i,Q2,CMSENE); } } else if (pdfOpt == 1) { for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,PT2,CMSENE); f_x2[i] = f(x2,i,PT2,CMSENE); } } else if (pdfOpt == 2) { double PT24 = PT2/4.; for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,PT24,CMSENE); f_x2[i] = f(x2,i,PT24,CMSENE); } } int I1=ID1; int I2=ID2; if( I1 == 21) I1=0; if( I2 == 21) I2=0; W[0][0] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, -1, P, KEY); W[0][1] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4, 1, 1, P, KEY); W[1][0] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, -1, P, KEY); W[1][1] = f_x1[I1]*f_x2[I2]*vbfdistrWRAP(ID1,ID2,ID3,ID4,-1, 1, P, KEY); double sumW =(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); if( sumW < 0 ){ std::cout << "ERW: f1*f2= " << f_x1[I1]*f_x2[I2] << " x1=" << x1 << " x2=" << x2 << " f1=" << f_x1[I1] << " f2=" << f_x2[I1] << std::endl; std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: W[] = " << W[0][0] << " " << W[0][1] << " " << W[1][0] << " " << W[1][1] << std::endl; std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: p1 = " << CMSENE/2*x1 << " " << 0.0 << " " << 0.0 << " " << CMSENE/2*x1 << std::endl; std::cout << "ERW: p2 = " << CMSENE/2*x2 << " " << 0.0 << " " << 0.0 << " " << -CMSENE/2*x2 << std::endl; std::cout << "ERW: X = " << sp_X.e() << " " << sp_X.px() << " " << sp_X.py() << " " << sp_X.pz() << std::endl; std::cout << "ERW: p3 = " << p3.e() << " " << p3.px() << " " << p3.py() << " " << p3.pz() << std::endl; std::cout << "ERW: p4 = " << p4.e() << " " << p4.px() << " " << p4.py() << " " << p4.pz() << std::endl; std::cout << "ERW: tau1 = " << tau1.e() << " " << tau1.px() << " " << tau1.py() << " " << tau1.pz() << std::endl; std::cout << "ERW: tau2 = " << tau2.e() << " " << tau2.px() << " " << tau2.py() << " " << tau2.pz() << std::endl; std::cout << "ERW: ==== " << std::endl; std::cout << "ERW: p1+p2= " << P[0][0]+P[1][0] << " " << P[0][1]+P[1][1] << " " << P[0][2]+P[1][2] << " " << P[0][3]+P[1][3] << std::endl; std::cout << "ERW: X+p3+p4= " << P[2][0]+P[3][0]+P[4][0]+P[5][0] << " " << P[2][1]+P[3][1]+P[4][1]+P[5][1] << " " << P[2][2]+P[3][2]+P[4][2]+P[5][2] << " " << P[2][3]+P[3][3]+P[4][3]+P[5][3] << std::endl; double p3sqr = p3.e()*p3.e()-p3.px()*p3.px()-p3.py()*p3.py()-p3.pz()*p3.pz(); double p4sqr = p4.e()*p4.e()-p4.px()*p4.px()-p4.py()*p4.py()-p4.pz()*p4.pz(); std::cout << "ERW: p3^2, p4^2 =" << p3sqr << " " << p4sqr << std::endl; } } /******************************************************************************* Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay plus 2 jets. Determine taus decay channel, calculates all necessary components from decay HHp HHm and production W[][] for calculation of all weights, later calculates weights. Input: X four momentum of Z/Higgs may be larger than sum of tau1 tau2, missing component is assumed to be QED brem four momenta of outgoing partons and taus; vectors of decay products for first and second tau Hidden input: relWTnonSM, switch of what is WTnonSM Hidden output: WTamplitM or WTamplitP matrix elements of tau1 tau2 decays Polari - helicity attributed to taus, 100% correlations between tau+ and tau- accesible with getTauSpin() WTnonSM weight for introduction of matrix element due to nonSM in cross section, or just ME^2*PDF's for relWTnonS==false Explicit output: WT spin correlation weight *******************************************************************************/ double calculateWeightFromParticlesVBFPROD(int IDPROD, SimpleParticle &p3, SimpleParticle &p4,SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector &sp_tau1_daughters, vector &sp_tau2_daughters, int KEY) { SimpleParticle sp_tau; SimpleParticle sp_nu_tau; vector sp_tau_daughters; // First iteration is for tau plus, so the 'nu_tau' is tau minus if (sp_tau1.pdgid() == -15 ) { sp_tau = sp_tau1; sp_nu_tau = sp_tau2; sp_tau_daughters = sp_tau1_daughters; } else { sp_tau = sp_tau2; sp_nu_tau = sp_tau1; sp_tau_daughters = sp_tau2_daughters; } double *HHp, *HHm; // We use this to separate namespace for tau+ and tau- if(true) { // Create Particles from SimpleParticles Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() ); Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() ); Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() ); vector tau_daughters; // tau pdgid int tau_pdgid = sp_tau.pdgid(); // Create list of tau daughters for(unsigned int i=0; i "; for(unsigned int i=0;i tau_daughters; // tau pdgid int tau_pdgid = sp_tau.pdgid(); // Create list of tau daughters for(unsigned int i=0; i "; for(unsigned int i=0;i=2) { // now we re-calculate W using anomalous variant of matrix elements. calcXsect(IDPROD, p3, p4, sp_X, sp_tau1, sp_tau2, W, KEY); double sum2=(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); WTnonSM=sum2/sum; if(relWTnonSM==0) WTnonSM=sum2; } double WT = W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2])+ W[1][0]*(1-HHp[2])*(1+HHm[2])+W[1][1]*(1-HHp[2])*(1-HHm[2]); WT = WT/(W[0][0]+W[0][1]+ W[1][0]+W[1][1]); // we separate cross section into first tau helicity + and - parts. Second tau follow. // This must be tested especially for KEY >=2 // case for scalar (H) has flipped sign of second tau spin? Tests needed double RRR = Tauola::randomDouble(); Polari=1.0; if (RRR<(W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2]))/(W[0][0]*(1+HHp[2])*(1+HHm[2])+W[0][1]*(1+HHp[2])*(1-HHm[2])+W[1][0]*(1-HHp[2])*(1+HHm[2])+W[1][1]*(1-HHp[2])*(1-HHm[2]))) Polari=-1.0; // Print out some info about the channel DEBUG( cout<<" WT: "<4.0) { printf("WT is: %13.10f. Setting WT = 4.0\n",WT); WT = 4.0; } delete[] HHp; delete[] HHm; return WT; } } // namespace TauSpinner