include/TauSpinner/tau_reweight_lib.h

#ifndef _TAU_REWEIGHT_LIB_H_
#define _TAU_REWEIGHT_LIB_H_
 
// Debug mode
#ifdef DEBUG_MODE
#define DEBUG(arg) arg;
#else
#define DEBUG(arg)
#endif
 
// TAUOLA header
#include "Tauola/Tauola.h"
 
#include "TauSpinner/Tauola_wrapper.h"
#include "TauSpinner/SimpleParticle.h"
#include "TauSpinner/Particle.h"
 
// LHAPDF header
#include "LHAPDF/LHAPDF.h"
 
#include <vector>
#include <iostream>
using std::vector;
using std::cout;
using std::endl;
 
namespace TauSpinner {
 
/** Definition of REAL*8 FUNCTION DISTH(S,T,H1,H2) from disth.f calculating 
    SM glue glue-->Higgs --> tau tau*/
extern "C" double disth_(double *SVAR, double *COSTHE, int *TA, int *TB);
 
/** routine to modifiey electroweak FF */
 void  GSWadapt(int ID,double S, double cost, double GSWr[7],double GSWi[7]);
  /** routine to rotate R for mustraal frames */
  void Rotin(double phi, double R[4][4]);
/** Initialize TauSpinner
 
    Print info and set global variables */
void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE);
 
/**
  Initialize anomalous dipole moments 8=2*2*2 constants: 
  magnetic/electric * real/imaginary parts * for virtual photon/ for Z, 
  and its switches ifGSW, ifkorch, iqed
  ifGSW electroweak form-factors need to  be on for Korchin calculations?
*/
  void initialize_GSW( int _ifGSW, int _ifkorch, int _iqed,
                       double _ReAini,
                       double _ImAini,
                       double _ReBini,
                       double _ImBini,
                       double _ReXini,
                       double _ImXini,
                       double _ReYini,
                       double _ImYini);
 
/**
 Initialize anomalous gamma gamma couplings target  A and B, in sample A0i B0i
 if sample is of no anonmalous couplings set A0i B0i to zero.
 Also how much of gamma gamma process there should be, with respect to quarks.
 In sample GAMfraci. In target GAMfrac2i.
*/  
  void initialize_gamagama(
                       double _GAMfraci,
                       double _GAMfrac2i,
                       double _A0i,
                       double _B0i,
                       double _Ai,
                       double _Bi);
 
 
/**
User may want to modify  form-factors, this routine enables that.
*/
  void initialize_GSW_norm(int   _keyGSW,
                          double _ReGSW1,
                          double _ImGSW1,
                          double _ReGSW2,
                          double _ImGSW2,
                          double _ReGSW3,
                          double _ImGSW3,
                          double _ReGSW4,
                          double _ImGSW4,
                          double _ReGSW5,
                          double _ImGSW5,
                          double _ReGSW6,
                          double _ImGSW6);
   
/** Set flag for calculating relative (NONSM/SM) or absolute that is matrix element^2 (spin averaged), 
  weight for X-section calculated as by product in longitudinal polarization method. */
void setRelWTnonSM(int _relWTnonSM);
 
/** Set flag for type of effective Born kinematic approximation 
 1:  Mustraal frame (NLO like)   0:  old default, similar to Collins-Soper. */
void setFrameType(int _FrameType);
 
/** set flag for simple formula (just Breit Wigner) of Higgs cross section, 
    set its parameters: Higgs mass, width and normalization for Higgs born (BW) function */
 void setHiggsParameters(int jak, double mass, double width, double normalization);
 
/** Get Higgs mass, width and normalization of Higgs born function */
void getHiggsParameters(double *mass, double *width, double *normalization);
 
/**  Set flag defining what type of spin treatment was used in analyzed  sample */
void setSpinOfSample(bool _Ipol);
 
/**  Turn  non Standard Model calculation  on/off */
void setNonSMkey(int key);
 
/** set transverse components for Higgs spin density matrix */
void setHiggsParametersTR(double Rxx, double Ryy, double Rxy, double Ryx);
 
/** set multipliers for transverse components  for Drell-Yan spin density matrix */
void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx);
 
/** get transverse components  for Drell-Yan spin density matrix */
void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx);
 
/** get longitudinal (or time) - transverse  components  for Drell-Yan spin density matrix */
  void getZgamParametersL(double &Rzx, double &Rzy, double &Rzz, double &Rtx, double &Rty,  double &Rtz);
 
/** get transverse components  for Drell-Yan spin density matrix multipliers */
void getZgamMultipliersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx);
 
/** get polarimetric vectors*/
void getHH(double &HHmx, double &HHmy, double &HHmz, double &HHpx, double &HHpy, double &HHpz);
 
/** Get nonSM weight of production matrix element */
double getWtNonSM();
 
/** Get tau+ weight of its decay matrix element */
double getWtamplitP();
 
/** Get tau- weight of its decay matrix element */
double getWtamplitM();
 
/** Get tau Helicity
    Used after event  is analyzed to obtain  attributed tau helicity.
    Returns -1 or 1. Note that 0 is returned if attribution fails. Method
    assumes that sample has spin effects taken into account. Effects can 
    be taken into account with the help of spin weights. */
double getTauSpin();
 
/** Calculate weights.
 
  Determines decay channel, calculates all necessary components for 
  calculation of all weights, calculates weights.
  Function for W+/- and H+/-  */
double calculateWeightFromParticlesWorHpn(SimpleParticle &W, SimpleParticle &tau, SimpleParticle &nu_tau, vector<SimpleParticle> &tau_daughters);
 
/** Calculate weights.
 
  Determines decay channel, calculates all necessary components for 
  calculation of all weights, calculates weights.
  Function for H/Z */
double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector<SimpleParticle> &sp_tau1_daughters, vector<SimpleParticle> &sp_tau2_daughters);
 
/** Prepare kinematics for HH calculation
  
  Boost particles to effective pair rest frame, and rotate them so that tau is on Z axis.
  Then rotate again with theta2 phi2 so neutrino (first tau decay product) from tau decay is along Z. */
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector<Particle> &tau_daughters, double *phi2, double *theta2);
 
/** Calculate polarization vector.
 
  We use FORTRAN metdods to calculate HH.
  First decide what is the channel. After that, 4-vectors
  are moved to tau rest frame of tau.
  Polarimetric vector HH is then rotated using angles phi and theta.
  
  Order of the particles does not matter. */
double* calculateHH(int tau_pdgid, vector<Particle> &tau_daughters, double phi, double theta);
 
/**
 Get Longitudinal polarization
 
 Returns longitudinal polarization in Z/gamma* -> tau+ tau-
 S: invariant mass^2 of the bozon */
double getLongitudinalPolarization(double S, SimpleParticle &sp_tau, SimpleParticle &sp_nu_tau);
 
/** Check if the vector of particles match the listed order of pdgid's:
 
  1) Returns true if 'particles' contain all of the listed pdgid-s.
  2) If it does - 'particles' will be reordered so that they have
     the same order as listed pdgid-s.
 
  It is done so the order of particles is the same as the order mandatory
  for TAUOLA Fortran routines. */
bool channelMatch(vector<Particle> &particles, int p1, int p2=0, int p3=0, int p4=0, int p5=0, int p6=0);
 
//------------------------------------------------------------------------------
//-- Useful debug methods ------------------------------------------------------
//------------------------------------------------------------------------------
 
/** Prints out two vertices:
      W   -> tau, nu_tau
      tau -> tau_daughters */
void print(Particle &W, Particle &nu_tau, Particle &tau, vector<Particle> &tau_daughters);
 
/** Sums all 4-vectors of the particles on the list */
Particle *vector_sum(vector<Particle> &x);
 
} // namespace TauSpinner
#endif