#include "TauSpinner/tau_reweight_lib.h" #include "TauSpinner/Particle.h" #include "TauSpinner/vbfdistr.h" #include #include using namespace Tauolapp; namespace TauSpinner { extern double CMSENE; extern int nonSM2; 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; int _QCDdefault=1; int _QCDvariant=1; /** @brief Pointer to a function that modify rezult of alpha_ calc. in vbfdist * * By default this function is not used (pointer is NULL). * It can be changed by the user through TauSpinner::set_vbfdistrModif */ void (*alphasModif)(double, int, int) = NULL; /** @brief Set vbfdistrModif function */ void set_alphasModif(void (*function)(double, int, int) ) { alphasModif = function; } /** @brief Pointer to a function that modify rezult of Matrix Element calc. in vbfdist * * By default this function is not used (pointer is NULL). * It can be changed by the user through TauSpinner::set_vbfdistrModif */ double (*vbfdistrModif)(int, int, int, int, int, int, double[6][4], int, double) = NULL; /** @brief Set vbfdistrModif function */ void set_vbfdistrModif(double (*function)(int, int, int, int, int, int, double[6][4], int, double) ) { vbfdistrModif = function; } void setPDFOpt(int QCDdefault,int QCDvariant){ _QCDdefault=QCDdefault; _QCDvariant=QCDvariant; } int getPDFOpt(int KEY){ if (KEY==0 || KEY==1) return _QCDdefault; else return _QCDvariant; } void alphas(double Q2,int scalePDFOpt, int KEY){ if (alphasModif) { alphasModif(Q2,scalePDFOpt,KEY); } else { // any textbook LL calculation const double PI=3.14159265358979324; // number of flavors const double Nf = 5; // alphas_s at mZ const double AlphasMZ = 0.118; const double MZ =91.1876; // alphas_s at scale Q2 (GeV^2) double alfas=AlphasMZ / ( 1 + AlphasMZ/(4*PI) * (11 - 2./3 * Nf) * log(Q2/(MZ*MZ))); // test Alphas ( Q2 = 1000^2 ) = 0.0877445 if(scalePDFOpt==0) alfas = 0.118; if(params_r_.as != alfas){// we pass alpha_s to calculation of amplitudes params_r_.as = alfas; // reinitialize constants couplings for amplitude calculations vbf_reinit_(&KEY); } } } /** Wrapper to VBDISTR and place for interface to user provided modification*/ double vbfdistr(int I1, int I2, int I3, int I4, int H1, int H2, double P[6][4], int KEY) { double P_copy[6][4]; double original_result = 0.; memcpy(P_copy,P,sizeof(double)*6*4); if( KEY<2) { // SM mode return vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY); } else if( !vbfdistrModif) { printf("TauSpinner::vbfdistr: User function vbfdistrModif not declared. Setting WT_contrib = 0.0 Failed attempt with KEY = %i5. \n",KEY); return original_result; } else if( KEY<4) { // modification mode int KEY_BUF=KEY-2; original_result= vbfdistr_(&I1, &I2, &I3, &I4, &H1, &H2, P_copy, &KEY_BUF); return vbfdistrModif(I1,I2,I3,I4,H1,H2,P_copy,KEY,original_result); } else { // replacement mode return vbfdistrModif(I1,I2,I3,I4,H1,H2,P_copy,KEY,original_result); } } /** Get VBF ME2 Returns array W[2][2] */ //--------------------------------------------------------------------------------------------------------------------------------------------------- void getME2VBF(SimpleParticle &p3i, SimpleParticle &p4i, SimpleParticle &sp_X,SimpleParticle &tau1i, SimpleParticle &tau2i, 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 p3(p3i.px(),p3i.py(),p3i.pz(),p3i.e(),0); Particle p4(p4i.px(),p4i.py(),p4i.pz(),p4i.e(),0); double mtaul=sqrt(tau1i.e()*tau1i.e()-tau1i.px()*tau1i.px()-tau1i.py()*tau1i.py()-tau1i.pz()*tau1i.pz()); Particle tau1(tau1i.px(),tau1i.py(),tau1i.pz(),tau1i.e(),0); Particle tau2(tau2i.px(),tau2i.py(),tau2i.pz(),tau2i.e(),0); // TEST begin /* if(KEY>1){ // this is special arrangement to used alternative calculation to introduce effect of smearing on tau lepton pair (hidden in sp_X) due to showering // in this way we get weights for corelated configs; with and without the shower-kick. Particle P_tautau(tau1.px()+tau2.px(),tau1.py()+tau2.py(),tau1.pz()+tau2.pz(),tau1.e()+tau2.e(),0); tau1.boostToRestFrame(P_tautau); tau2.boostToRestFrame(P_tautau); p3.boostToRestFrame(P_tautau); p4.boostToRestFrame(P_tautau); //v1.boostToRestFrame(P_tautau); //v2.boostToRestFrame(P_tautau); Particle P_X(sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(),0); tau1.boostFromRestFrame(P_X); tau2.boostFromRestFrame(P_X); p3.boostFromRestFrame(P_X); p4.boostFromRestFrame(P_X); //v1.boostFromRestFrame(P_X); //v2.boostFromRestFrame(P_X); } */ // TEST end // we may want to force p3,p4 to be masless. Particle v1(0.,0., 1.,1.,0.); Particle v2(0.,0.,-1.,1.,0.); Particle P_QQ( p3.px()+p4.px()+tau1.px()+tau2.px()+1e-8*p3.pz(), p3.py()+p4.py()+tau1.py()+tau2.py()-2e-8*p3.pz(), p3.pz()+p4.pz()+tau1.pz()+tau2.pz(), p3.e() +p4.e() +tau1.e() +tau2.e(), 0 ); tau1.boostToRestFrame(P_QQ); tau2.boostToRestFrame(P_QQ); p3.boostToRestFrame(P_QQ); p4.boostToRestFrame(P_QQ); v1.boostToRestFrame(P_QQ); v2.boostToRestFrame(P_QQ); // Now we can define vers1 vers2 double xn1=sqrt(v1.px()*v1.px()+v1.py()*v1.py()+v1.pz()*v1.pz()); double xn2=sqrt(v2.px()*v2.px()+v2.py()*v2.py()+v2.pz()*v2.pz()); double xn12=sqrt( (v1.px()-v2.px())*(v1.px()-v2.px()) +(v1.py()-v2.py())*(v1.py()-v2.py())+(v1.pz()-v2.pz())*(v1.pz()-v2.pz())); 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! // four options for partonic beams. double P[6][4] = { { sqrt(SS)/2, sqrt(SS)/2/xn12*(v1.px()-v2.px()), sqrt(SS)/2/xn12*(v1.py()-v2.py()), sqrt(SS)/2/xn12*(v1.pz()-v2.pz()) }, { sqrt(SS)/2, -sqrt(SS)/2/xn12*(v1.px()-v2.px()), -sqrt(SS)/2/xn12*(v1.py()-v2.py()), -sqrt(SS)/2/xn12*(v1.pz()-v2.pz()) }, // double P[6][4] = { { sqrt(SS)/2, -sqrt(SS)/2/xn1*v2.px(), -sqrt(SS)/2/xn1*v2.py(), -sqrt(SS)/2/xn1*v2.pz() }, // { sqrt(SS)/2, sqrt(SS)/2/xn1*v2.px(), sqrt(SS)/2/xn1*v2.py(), sqrt(SS)/2/xn1*v2.pz() }, // double P[6][4] = { { sqrt(SS)/2, sqrt(SS)/2/xn1*v1.px(), sqrt(SS)/2/xn1*v1.py(), sqrt(SS)/2/xn1*v1.pz() }, // { sqrt(SS)/2, -sqrt(SS)/2/xn1*v1.px(), -sqrt(SS)/2/xn1*v1.py(), -sqrt(SS)/2/xn1*v1.pz() }, // double P[6][4] = { { sqrt(SS)/2, 0.0, 0.0, sqrt(SS)/2 }, // { sqrt(SS)/2, 0.0, 0.0, -sqrt(SS)/2 }, { 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_tautau( tau1.px()+tau2.px(), tau1.py()+tau2.py(), tau1.pz()+tau2.pz(), tau1.e()+tau2.e(),0 ); double Q2 = P_tautau.recalculated_mass()*P_tautau.recalculated_mass(); double QQ2 = P_QQ.recalculated_mass()*P_QQ.recalculated_mass(); double PT2 = P_QQ.px() * P_QQ.px() + P_QQ.py()* P_QQ.py(); double sumET = p3.e()*sqrt( p3.px()*p3.px()+p3.py()*p3.py())/sqrt( p3.px()*p3.px()+p3.py()*p3.py()+p3.pz()*p3.pz()) + p4.e()*sqrt( p4.px()*p4.px()+p4.py()*p4.py())/sqrt( p4.px()*p4.px()+p4.py()*p4.py()+p4.pz()*p4.pz()) + tau1.e()*sqrt( tau1.px()*tau1.px()+tau1.py()*tau1.py())/sqrt( tau1.px()*tau1.px()+tau1.py()*tau1.py()+tau1.pz()*tau1.pz() ) + tau2.e()*sqrt( tau2.px()*tau2.px()+tau2.py()*tau2.py())/sqrt( tau2.px()*tau2.px()+tau2.py()*tau2.py()+tau2.pz()*tau2.pz() ) ; double sumMT = sqrt( p3.recalculated_mass()*p3.recalculated_mass() + p3.px()*p3.px()+p3.py()*p3.py() ) + sqrt( p4.recalculated_mass()*p4.recalculated_mass() + p4.px()*p4.px()+p4.py()*p4.py() ) + sqrt( tau1.recalculated_mass()*tau1.recalculated_mass() + tau1.px()*tau1.px()+tau1.py()*tau1.py() ) + sqrt( tau2.recalculated_mass()*tau2.recalculated_mass() + tau2.px()*tau2.px()+tau2.py()*tau2.py() ); double fixed_scale = 91.17; // 200.; int scalePDFOpt = 1; // it has to be connected to initialization and flipper scalePDFOpt=getPDFOpt(KEY); double Q2pdf = fixed_scale*fixed_scale; if(scalePDFOpt == 0){ // alphas(Q2pdf,scalePDFOpt,KEY); } else if (scalePDFOpt == 1) { Q2pdf = QQ2; alphas(Q2pdf,scalePDFOpt,KEY); } else if (scalePDFOpt == 2) { Q2pdf = sumMT*sumMT; alphas(Q2pdf,scalePDFOpt,KEY); } else if (scalePDFOpt == 3) { Q2pdf = sumET*sumET; alphas(Q2pdf,scalePDFOpt,KEY); } else if (scalePDFOpt == 4) { Q2pdf = Q2; alphas(Q2pdf,scalePDFOpt,KEY); } for(int i=-5;i<=5;++i) { f_x1[i] = f(x1,i,Q2pdf,CMSENE); f_x2[i] = f(x2,i,Q2pdf,CMSENE); } // reset to zero W[0][0]=0.0; W[0][1]=0.0; W[1][0]=0.0; W[1][1]=0.0; if( DEBUG ){ // here you can overwrite kinematical configurations for tests std::cout << " " << std::endl; std::cout << "---podstawiamy 4-pedy na jakies ---" << std::endl; 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 } }; 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 } }; 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 } }; 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 } }; for (int I1=0;I1<=5;I1++){for (int I2=0;I2<=3;I2++){ P[I1][I2]=P1[I1][I2]; }} printf(" our event : P[0,i]= %16.10f %16.10f %16.10f %16.10f \n",P[0][0],P[0][1],P[0][2],P[0][3]); int II=1; printf(" P[1,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]); II=2; printf(" P[2,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]); II=3; printf(" P[3,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]); II=4; printf(" P[4,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]); II=5; printf(" P[4,i]= %16.10f %16.10f %16.10f %16.10f \n",P[II][0],P[II][1],P[II][2],P[II][3]); printf(" ===========================\n"); //################################################################# KEY=0; //################################################################# //################################################################# printf(" ============== KEY= %2i ==\n",KEY); printf(" ===========================\n"); printf(" \n"); printf(" ===========================================================\n"); printf(" ============== non-zero contributions to <|ME|^2>_spin ==\n"); printf(" ===========================================================\n"); printf(" \n"); } // these loops need to be cleaned from zero contributions! for(int I1=-5;I1<=5;I1++){ // for test of single production process fix flavour for(int I2=-5;I2<=5;I2++){ // for test of single production process fix flavour for(int I3=-5;I3<=5;I3++){ for(int I4=-5;I4<=5;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; W[0][0] += f_x1[I1]*f_x2[I2]*vbfdistr(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]*vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY); // second tau helicity assuming without checking that W[1][0] += f_x1[I1]*f_x2[I2]*vbfdistr(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]*vbfdistr(ID1,ID2,ID3,ID4,-1, 1, P, KEY); if( DEBUG ) { if (ID1==51 && ID2==-1 && ID3==4 && ID4==-4 ) /// THIS IS BLOCKED NOW ##################################################### { KEY=1; double cosik=vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4,-1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4,-1, 1, P, KEY); std::cout << " ID-s = " < 10) || (W[0][1] > 10) || (W[1][0] > 10) || ( W[1][1] > 10)) { std::cout << "ERWxxx: ID= " < 0) { float check = vbfdistr(ID1,ID2,ID4,ID3, 1, -1, P, KEY)+vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY); if( check != 0 || check == 0){ int CPSTATUSICP=0; printf(" ID-s= %2i %2i %2i %2i CP used= %1i ## VALUE: %16.10e ## Spin contr.: (+-)= %9.3e (-+)= %9.3e (--)= %9.3e (++)= %9.3e \n", ID1, ID2, ID3,ID4, cpstatus_.icp, vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, -1, -1, P, KEY)+ vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY), vbfdistr(ID1,ID2,ID3,ID4, 1, -1, P, KEY), vbfdistr(ID1,ID2,ID3,ID4, -1, 1, P, KEY), vbfdistr(ID1,ID2,ID3,ID4, -1, -1, P, KEY), vbfdistr(ID1,ID2,ID3,ID4, 1, 1, P, KEY) ); } } } } } } } // std::cout << "That is it"<< std::endl; //exit(-1); } /******************************************************************************* 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 calculateWeightFromParticlesVBF(SimpleParticle &p3, SimpleParticle &p4,SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector &sp_tau1_daughters, vector &sp_tau2_daughters) { SimpleParticle sp_tau; SimpleParticle sp_nu_tau; vector sp_tau_daughters; int KEY=0; // here we impose that it is for SM Higgs only if(sp_X.pdgid()==25) KEY=1; // 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 // 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