#include "TauSpinner/tau_reweight_lib.h" #include "TauSpinner/nonSM.h" #include #include "Tauola/TauolaParticlePair.h" using namespace Tauolapp; namespace TauSpinner { /***** GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/ double CMSENE = 7000.0;// Center of mass system energy (used by PDF's) bool Ipp = true; // pp collisions int Ipol = 1; // Is the sample in use polarized? (relevant for Z/gamma case only) int nonSM2 = 0; // Turn on/off nonSM calculations int nonSMN = 0; // Turn on, if calculations of nonSM weight, is to be used to modify shapes only int relWTnonSM = 1; // 1: relWTnonSM is relative to SM; 0: absolute double WTnonSM=1.0; // nonSM weight double Polari =0.0; // Helicity, attributed to the tau pair. If not attributed then 0. // Polari is attributed for the case when spin effects are taken into account. // Polari is available for the user program with getTauSpin() bool IfHiggs=false; // flag used in sigborn() double WTamplit=1; // AMPLIT weight for the decay last invoked double WTamplitP=1; // AMPLIT weight for tau+ decay double WTamplitM=1; // AMPLIT weight for tau- decay // Higgs parameters int IfHsimple=0; // switch for simple Higgs (just Breit-Wigner) double XMH = 125.0; // mass (GeV) double XGH = 1.0; // width, should be detector width, analysis dependent. double Xnorm = 0.15; // normalization of Higgs Born cross-section, at hoc // Transverse components of Higgs spin density matrix double RXX = 0.0; //-1.0; double RYY = 0.0; // 1.0; double RXY = 0.0; double RYX = 0.0; // Coefficients for transverse components of Z/gamma spin density matrix double RzXX = 0.0; //-1.0; double RzYY = 0.0; // 1.0; double RzXY = 0.0; double RzYX = 0.0; // Values of transverse components of Z/gamma spin density matrix calculated inside getLongitudinalPolarization double R11 = 0.0; double R22 = 0.0; double R12 = 0.0; // for future use double R21 = 0.0; // for future use /***** END: GLOBAL VARIABLES AND THEIR PRE-INITIALISATION *****/ double f(double x, int ID, double SS, double cmsene) // PDF's parton density function divided by x; // x - fraction of parton momenta // ID flavour of incoming quark // SS scale of hard process // cmsene center of mass for pp collision. { // LHAPDF manual: http://projects.hepforge.org/lhapdf/manual // double xfx(const double &x;, const double &Q;, int fl); // returns xf(x, Q) for flavour fl - this time the flavour encoding // is as in the LHAPDF manual... // -6..-1 = tbar,...,ubar, dbar // 1..6 = duscbt // 0 = g // xfx is the C++ wrapper for fortran evolvePDF(x,Q,f) // note that SS=Q^2, make the proper choice of PDFs arguments. return LHAPDF::xfx(x, sqrt(SS), ID)/x; //return x*(1-x); } // Calculates Born cross-section summed over final taus spins. // Input parameters: // incoming flavour ID // invariant mass^2 SS // scattering angle costhe // Hidden input (type of the process): IfHiggs, IfHsimple double sigborn(int ID, double SS, double costhe) { // cout << "ID : " << ID << " HgsWTnonSM = " << HgsWTnonSM << " IfHsimple = " << IfHsimple << endl; // BORN x-section. // WARNING: overall sign of costheta must be fixed int tauID = 15; // case of the Higgs boson if (IfHiggs) { double SIGggHiggs=0.; // for the time being only for gluon it is non-zero. if(ID==0){ int IPOne = 1; int IMOne = -1; SIGggHiggs=disth_(&SS, &costhe, &IPOne, &IPOne)+disth_(&SS, &costhe, &IPOne, &IMOne)+ disth_(&SS, &costhe, &IMOne, &IPOne)+disth_(&SS, &costhe, &IMOne, &IMOne); double PI=3.14159265358979324; SIGggHiggs *= XMH * XMH * XMH *XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) + XGH*XGH*XMH*XMH); // cout << "JK: SIGggHiggs = " << SS << " " << XMH << " " << XGH << " " << XMH * XMH * XMH * XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) + XGH*XGH*XMH*XMH) << " " << SIGggHiggs << endl; if(IfHsimple==1) SIGggHiggs = Xnorm * XMH * XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) +XGH*XGH*XMH*XMH); // cout << "ZW: SIGggHiggs = " << SS << " " << costhe << " " << SIGggHiggs << endl; } return SIGggHiggs; } // case of Drell-Yan if (ID==0) return 0.0 ; // for the time being for gluon it is zero. if (ID>0) initwk_( &ID, &tauID, &SS); else { ID = -ID; initwk_( &ID, &tauID, &SS); } int iZero = 0; double dOne = 1.0; double dMOne = -1.0; // sum DY Born over all tau helicity configurations: return ( t_born_(&iZero, &SS, &costhe, &dOne , &dOne) + t_born_(&iZero, &SS, &costhe, &dOne , &dMOne) + t_born_(&iZero, &SS, &costhe, &dMOne, &dOne) + t_born_(&iZero, &SS, &costhe, &dMOne, &dMOne))/SS/128.86674175/128.86674175/3.141559265358979324; //123231.; // overall norm. factor .../SS/123231 most probably it is alpha_QED**2/pi/2/SS is from comparison between Born we use and Born used in Warsaw group. } /******************************************************************************* Initialize TauSpinner Print info and set global variables *******************************************************************************/ void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE) { Ipp = _Ipp; Ipol = _Ipol; nonSM2 = _nonSM2; nonSMN = _nonSMN; CMSENE = _CMSENE; cout<<" ------------------------------------------------------"<=0 "<"< tau nu_tau decay. Function for W+/- and H+/- Determines decay channel, calculates all necessary components for calculation of all weights, calculates weights. Input: X four momentum may be larger than sum of tau nu_tau, missing component is assumed to be QED brem Hidden input: none Hidden output: WTamplitM or WTamplitP of tau decay (depending on tau charge) NOTE: weight for sp_X production matrix elements is not calculated for decays of charged intermediate W+-/H+-! Explicit output: WT spin correlation weight *******************************************************************************/ double calculateWeightFromParticlesWorHpn(SimpleParticle &sp_X, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau, vector &sp_tau_daughters) { // 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 vector of tau daughters for(unsigned int i=0; i "; for(unsigned int i=0;i2.0) { printf("Tauspinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 2.0\n",WT); WT = 2.0; } delete HH; return WT; } /******************************************************************************* Calculate weights, case of event record vertex like Z/gamma/H ... -> tau tau decay. Determine decay channel, calculates all necessary components for calculation of all weights, calculates weights. Input: X four momentum may be larger than sum of tau1 tau2, missing component is assumed to be QED brem Hidden input: relWTnonS, nonSM2 (used only in getLongitudinalPolarization(...) ) Hidden output: WTamplitM or WTamplitP of tau1 tau2 decays weight for sp_X production matrix element^2 is calculated inside plzap2 Polari - helicity attributed to taus, 100% correlations between tau+ and tau- WTnonSM Explicit output: WT spin correlation weight *******************************************************************************/ double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector &sp_tau1_daughters, vector &sp_tau2_daughters) { // cout << "sp_tau1_daughters = " << sp_tau1_daughters.size() << endl; // cout << "sp_tau2_daughters = " << sp_tau2_daughters.size() << endl; SimpleParticle sp_tau; SimpleParticle sp_nu_tau; vector sp_tau_daughters; // First we calculate HH for tau+ // We enforce that sp_tau is tau+ so the 'nu_tau' is tau- 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 artificial if(true){... } construction 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