TEvtGen/Tauola/Plots.cxx

   1 #include "Plots.h"
   2 #include <iostream>
   3 #include <fstream>
   4 #include <math.h>
   5 using namespace std;
   6
   7 namespace Tauolapp
   8 {
   9
  10 Plots::Plots():
  11     m_incoming_pdg_id(1),
  12     m_cosTheta       (-0.2),
  13     m_n_plot_points  (1000)
  14 {
  15 }
  16
  17 void Plots::SANCtest1(){
  18
  19   cout<<"SANC plot 1 (short)..."<<endl;
  20
  21   double smin = log(6.*6.)+0.0001;
  22   double smax = log(17000.*17000.);
  23   double step = (smax-smin)/(m_n_plot_points-1);
  24
  25   ofstream f1,f2,f3;
  26   f1.open("f-sanc.txt");
  27   f2.open("f-born.txt");
  28   f3.open("f-plzap0.txt");
  29   for(int i=0; i<m_n_plot_points; i++)
  30   {
  31     double s = exp(smin+i*step);
  32
  33     // Write SANC value
  34     t_pair.recalculateRij(m_incoming_pdg_id,15,s,m_cosTheta);
  35     f1<<sqrt(s)<<" "<<t_pair.m_R[0][3]<<endl;
  36
  37     // Write Born-level value
  38     // (assuming table 11 is filled with born-level data)
  39     t_pair.recalculateRij(11,15,s,m_cosTheta);
  40     f2<<sqrt(s)<<" "<<t_pair.m_R[0][3]<<endl;
  41
  42     int outgoing_pdg_id = 15;
  43
  44     // Write PLZAP0 value
  45     double pz = 1-plzap0_(&m_incoming_pdg_id,&outgoing_pdg_id,&s, &m_cosTheta);
  46     t_pair.m_R[0][3]=2*pz-1;
  47     f3<<sqrt(s)<<" "<<t_pair.m_R[0][3]<<endl;
  48   }
  49   f1.close();
  50   f2.close();
  51   f3.close();
  52 }
  53
  54 void Plots::SANCtest2(){
  55
  56   cout<<"SANC plot 2 (short)..."<<endl;
  57
  58   double smin = log(6.*6.)+0.0001;
  59   double smax = log(17000.*17000.);
  60   double step = (smax-smin)/(m_n_plot_points-1);
  61
  62   ofstream f1,f2,f3;
  63   f1.open("f-w-single-point.txt");
  64   f2.open("f-w0-single-point.txt");
  65   f3.open("f-ww0-single-point.txt");
  66
  67   for(int i=0; i<m_n_plot_points; i++){
  68
  69     double s=exp(smin+i*step);
  70     t_pair.recalculateRij(m_incoming_pdg_id,15,s,m_cosTheta);
  71
  72     // Write w, w0 and w/w0
  73     f1<<sqrt(s)<<" "<<Tauola::getEWwt()<<endl;
  74     f2<<sqrt(s)<<" "<<Tauola::getEWwt0()<<endl;
  75     f3<<sqrt(s)<<" "<<Tauola::getEWwt()/Tauola::getEWwt0()<<endl;
  76   }
  77   f1.close();
  78   f2.close();
  79   f3.close();
  80 }
  81
  82 void Plots::SANCtest3(){
  83
  84   cout<<"SANC plot 3 (long)..."<<endl;
  85
  86   double smin = log(6.*6.)+0.0001;
  87   double smax = log(17000.*17000.);
  88   double step = (smax-smin)/(m_n_plot_points-1);
  89
  90   ofstream f1;
  91   f1.open("f-err.txt");
  92   double costheta=-1.;
  93
  94   for(int i=0; i<m_n_plot_points; i++){
  95
  96     double buf=0.,err=0.;
  97
  98     for(int j=0; j<m_n_plot_points; j++){
  99
 100       double s = exp(smin+j*step);
 101
 102       // Get value from SANC table
 103       t_pair.recalculateRij(m_incoming_pdg_id,15,s,costheta);
 104       buf = t_pair.m_R[0][3];
 105       t_pair.recalculateRij(11,15,s,costheta);
 106
 107       // Calculate error between SANC and Born-level
 108       err += (buf-t_pair.m_R[0][3])*(buf-t_pair.m_R[0][3]);
 109     }
 110
 111     f1<<costheta<<" "<<err/m_n_plot_points<<endl;
 112     err=0.;
 113     costheta+=2./m_n_plot_points;
 114   }
 115
 116   f1.close();
 117 }
 118
 119 void Plots::SANCtest4(){
 120
 121   cout<<"SANC plot 4 (medium)..."<<endl;
 122
 123   double smin = log(6.*6.);
 124   double smax = log(17000.*17000.);
 125   double step = (smax-smin)/(m_n_plot_points-1);
 126
 127   ofstream f1,f2,f3;
 128   f1.open("f-cross.txt");
 129   f2.open("f-w.txt");
 130   f3.open("f-w0.txt");
 131
 132   for(int i=0; i<m_n_plot_points; i++){
 133
 134     double s        =  exp(smin+i*step);
 135     double sumEWwt  =  0.;
 136     double sumEWwt0 =  0.;
 137     double costheta = -1.;
 138     int    NCOS     =  21;
 139
 140     // Calculate sum of w/w0 for all cosTheta
 141     for(int j=0; j<NCOS; j++){
 142
 143       costheta = -1. + 1.0/NCOS + j*2./NCOS;
 144
 145       t_pair.recalculateRij(m_incoming_pdg_id,15,s,costheta);
 146
 147       sumEWwt +=Tauola::getEWwt();
 148       sumEWwt0+=Tauola::getEWwt0();
 149     }
 150
 151     f1<<sqrt(s)<<" "<<sumEWwt/sumEWwt0/m_n_plot_points<<endl;
 152     f2<<sqrt(s)<<" "<< 2./NCOS * sumEWwt  <<endl;
 153     f3<<sqrt(s)<<" "<< 2./NCOS * sumEWwt0 <<endl;
 154   }
 155
 156   f1.close();
 157   f2.close();
 158   f3.close();
 159 }
 160
 161 void Plots::setSancVariables(int incoming, double cosTheta) {
 162   m_incoming_pdg_id = incoming;
 163   m_cosTheta        = cosTheta;
 164 }
 165
 166 } // namespace Tauolapp