C++ Interface to Tauola
nonSM.cxx
1#include <iostream>
2#include "TauSpinner/nonSM.h"
3#include "TauSpinner/EWtables.h"
4using std::cout;
5using std::endl;
6
7namespace TauSpinner {
8
9/***** GLOBAL VARIABLES *****/
10
11double (*nonSM_bornZ)(int, double, double, int, int, int) = default_nonSM_born;
12double (*nonSM_bornH)(int, double, double, int, int, int) = default_nonSM_bornH;
13
14extern int nonSM2;
15extern int nonSMN;
16extern double WTnonSM;
17extern bool IfHiggs;
18/********************************/
19
20double nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
21{
22 if(IfHiggs) return nonSM_bornH(ID,S,cost,H1,H2,key);
23 else return nonSM_bornZ(ID,S,cost,H1,H2,key);
24}
25
26
27double default_nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
28{
29 if (ID < 0)
30 {
31 ID = -ID;
32 cost = -cost;
33 }
34
35 if (ID == 3 || ID == 5)
36 {
37 ID = 1;
38 }
39 else if (ID == 4)
40 {
41 ID = 2;
42 }
43
44 if (ID >= 3 )
45 {
46 cout << "WARNING: default_nonSM_born: unexpected ID continue debugging: " << ID << endl;
47 }
48
49 // WARNING: code below compiles/links but was never tested fir numerical
50 // results. Comments in polish are of the work plan
51 if (initEWff(ID,S,cost,key) == 1){ // initEWff check if initialization
52 // init() was executed
53 // and only if ID,S,cost changed
54 // interpolate what bornew needs.
55 // chyba juz wiem co ze zmiennymi
56 // smiecoiwymi sigbornswdelt,
57 // lokalnie inicjalizowanymi jak:
58 // AMZi GAMMZi, SWeff, DeltSQ
59 // DeltV, Gmu, alfinv
60 // ktore przynajmniej po czesci sa pod
61 // innymi nazwami i z (prawie) takimi
62 // samymi wartosciami.
63 // If no intialization for bornew then
64 // exit(-1) and program stop.
65 int Bmode=1; // locally defined const., to change user may pointer-activate his own nonSM_born()
66 int keyGSW=1; // locally defined const, see T_BORNEW_ header, note that keGSW 1/2 is distinguished only by init
67 // for comments and meaning of these two parameters
68 // at present used for some specialized test plots.
69
70 double H1d = H1;
71 double H2d = H2;
72 int zero = 0;
73 if(key!=0){ //cout<<" key!0 " << H1<<" " << H2 << " " << t_bornew_(&Bmode, &keyGSW, &S, &cost, &H1d , &H2d) <<endl;
74 return t_bornew_(&Bmode, &keyGSW, &S, &cost, &H1d , &H2d);}
75 else { // cout<<" key=0 " << H1<< H2 << " " << t_born_ (&zero, &S, &cost, &H1d , &H2d) <<endl;
76 return t_born_ (&zero, &S, &cost, &H1d , &H2d);}
77 }
78 else
79 {
80 cout<<"TauSpinner::default_nonSM_born: EW tables were not initialize\n"<<endl;
81 cout<<" and user has not provided his own nonSM_born too"<<endl;
82 cout<<" see: nonSM_adopt() and set_nonSM_born( NULL )"<<endl;
83 cout<<" in TauSpinner/examples/tau-reweight-test.cxx for details"<<endl;
84 exit(-1);
85 return 1.0;
86 }
87}
88
89
90double default_nonSM_bornH(int ID, double S, double cost, int H1, int H2, int key)
91{
92 cout<<"TauSpinner::default_nonSM_born: this function is dummy\n"<<endl;
93 cout<<" user must provide his own nonSM_born"<<endl;
94 cout<<" see: nonSM_adopt() and set_nonSM_born( NULL )"<<endl;
95 cout<<" in TauSpinner/examples/tau-reweight-test.cxx for details"<<endl;
96 exit(-1);
97 return 1.0;
98}
99
100int set_nonSM_born( double (*fun)(int, double, double, int, int, int) )
101{
102 if(fun==NULL) {nonSM_bornZ = default_nonSM_born; return 0;}
103 else {nonSM_bornZ = fun; return 1;}
104}
105
106int set_nonSM_bornH( double (*fun)(int, double, double, int, int, int) )
107{
108 if(fun==NULL) {nonSM_bornH = default_nonSM_bornH; return 0;}
109 else {nonSM_bornH = fun; return 1;}
110}
111
112/*******************************************************************************
113 Probability of the helicity state
114*******************************************************************************/
115double plzap2(int ide, int idf, double svar, double costhe)
116{
117 if(ide!=0)
118 {
119 if(idf > 0) initwk_(&ide,&idf,&svar);
120 else
121 {
122 int mide=-ide;
123 int midf=-idf;
124
125 initwk_(&mide,&midf,&svar);
126 }
127 }
128
129 int zero = 0;
130 double one = 1.0;
131 double mone = -1.0;
132 double ret = 0.0;
133
134 if(nonSM2==0)
135 {
136 ret = t_born_(&zero,&svar,&costhe, &one, &one)
137 /(t_born_(&zero,&svar,&costhe, &one, &one)
138 +t_born_(&zero,&svar,&costhe,&mone,&mone));
139 }
140 else if(nonSM2==1)
141 {
142 ret = nonSM_born(ide,svar,costhe, 1, 1, nonSM2)
143 /(nonSM_born(ide,svar,costhe, 1, 1, nonSM2)
144 +nonSM_born(ide,svar,costhe,-1,-1, nonSM2));
145
146 // test of user prepared born cross section can be prepared here.
147 // convention between t_born and nonSM_born in choice of flavours
148 // sign of costhe helicity signs etc have to be performed using SM version
149 // of nonSM_born. Matching (up to may be overall s-dependent factor) between
150 // t_born and nonSM_born must be achieved, see Section 4 for details
151 // on technical tests.
152 DEBUG(
153 double sm= t_born_(&zero,&svar,&costhe, &one, &one)
154 /(t_born_(&zero,&svar,&costhe, &one, &one)
155 +t_born_(&zero,&svar,&costhe,&mone,&mone));
156 double nsm= nonSM_born(ide,svar,costhe, 1, 1, nonSM2)
157 /(nonSM_born(ide,svar,costhe, 1, 1, nonSM2)
158 +nonSM_born(ide,svar,costhe,-1,-1, nonSM2));
159 double smn= nonSM_born(ide,svar,costhe, 1, 1, 0 )
160 /(nonSM_born(ide,svar,costhe, 1, 1, 0 )
161 +nonSM_born(ide,svar,costhe,-1,-1, 0 ));
162
163 cout<<"test of nonSM Born nonsm2="<<nonSM2 << endl;
164 cout<<"ide,svar,costhe="<<ide <<" " << svar <<" " << costhe << endl;
165 cout<<"sm="<<sm <<" sm (new)="<<smn <<" nsm="<<nsm << endl;
166 cout<<"sm and sm (new) should be essentially equal" << endl << endl;
167 if (IfHiggs) cout << "(for Higgs we need to improve algorithm)" << endl;
168
169 )
170
171 }
172 return ret;
173}
174
175double plweight(int ide, double svar, double costhe)
176{
177 if(nonSM2==0) return 1.0;
178 if (ide==0 && !IfHiggs) return 1.0;
179
180 double ret =( nonSM_born(ide,svar,costhe, 1, 1, nonSM2)/svar
181 +nonSM_born(ide,svar,costhe,-1,-1, nonSM2)/svar);
182 double retm =( nonSM_born(ide,svar,costhe, 1, 1, 0)/svar
183 +nonSM_born(ide,svar,costhe,-1,-1, 0)/svar); // svar is introduced for future use
184 // another option angular dependence only. Effect on x-section removed
185 if(nonSMN==1) ret = ret * plnorm(ide,svar);
186 return ret/retm;
187}
188
189double plnorm(int ide, double svar)
190{
191 if(nonSMN==0) return 1.0;
192
193 double c1 = 1.0/sqrt(3.0);
194 double c2 = sqrt(2.0/3.0);
195
196 double alpha = 2*nonSM_born(ide,svar,0.0, 1, 1,0)+
197 2*nonSM_born(ide,svar,0.0,-1,-1,0);
198 double beta = nonSM_born(ide,svar, c1, 1, 1,0) + nonSM_born(ide,svar,-c1, 1, 1,0)+
199 nonSM_born(ide,svar, c1,-1,-1,0) + nonSM_born(ide,svar,-c1,-1,-1,0);
200 double gamma = nonSM_born(ide,svar, c2, 1, 1,0) + nonSM_born(ide,svar,-c2, 1, 1,0)+
201 nonSM_born(ide,svar, c2,-1,-1,0) + nonSM_born(ide,svar,-c2,-1,-1,0);
202 double ret = ( alpha + 0.9*(gamma+alpha-2*beta) + 0.5*(4*beta-3*alpha-gamma) );
203
204 alpha = 2*nonSM_born(ide,svar,0.0, 1, 1,nonSM2)+
205 2*nonSM_born(ide,svar,0.0,-1,-1,nonSM2);
206 beta = nonSM_born(ide,svar, c1, 1, 1,nonSM2) + nonSM_born(ide,svar,-c1, 1, 1,nonSM2)+
207 nonSM_born(ide,svar, c1,-1,-1,nonSM2) + nonSM_born(ide,svar,-c1,-1,-1,nonSM2);
208 gamma = nonSM_born(ide,svar, c2, 1, 1,nonSM2) + nonSM_born(ide,svar,-c2, 1, 1,nonSM2)+
209 nonSM_born(ide,svar, c2,-1,-1,nonSM2) + nonSM_born(ide,svar,-c2,-1,-1,nonSM2);
210
211 ret = ret / ( alpha + 0.9*(gamma+alpha-2*beta) + 0.5*(4*beta-3*alpha-gamma) );
212
213 return ret;
214}
215
216void nonSMHcorrPol(double S, SimpleParticle &tau1, SimpleParticle &tau2,
217 double *corrX2, double *polX2)
218{ // tau+ and tau- in lab frame
219 Particle tau_plus ( tau1.px(), tau1.py(), tau1.pz(), tau1.e(), tau1.pdgid() );
220 Particle tau_minus( tau2.px(), tau2.py(), tau2.pz(), tau2.e(), tau2.pdgid() );
221
222 // P_QQ = sum of tau+ and tau- in lab frame
223 Particle P_QQ( tau_plus.px()+tau_minus.px(), tau_plus.py()+tau_minus.py(), tau_plus.pz()+tau_minus.pz(), tau_plus.e()+tau_minus.e(), 0 );
224
225 Particle P_B1(0, 0, 1, 1, 0);
226 Particle P_B2(0, 0,-1, 1, 0);
227
228 tau_plus. boostToRestFrame(P_QQ);
229 tau_minus.boostToRestFrame(P_QQ);
230 P_B1. boostToRestFrame(P_QQ);
231 P_B2. boostToRestFrame(P_QQ);
232
233 double costheta1 = (tau_plus.px()*P_B1.px() +tau_plus.py()*P_B1.py() +tau_plus.pz()*P_B1.pz() ) /
234 sqrt(tau_plus.px()*tau_plus.px()+tau_plus.py()*tau_plus.py()+tau_plus.pz()*tau_plus.pz()) /
235 sqrt(P_B1.px() *P_B1.px() +P_B1.py() *P_B1.py() +P_B1.pz() *P_B1.pz() );
236
237 double costheta2 = (tau_minus.px()*P_B2.px() +tau_minus.py()*P_B2.py() +tau_minus.pz()*P_B2.pz() ) /
238 sqrt(tau_minus.px()*tau_minus.px()+tau_minus.py()*tau_minus.py()+tau_minus.pz()*tau_minus.pz()) /
239 sqrt(P_B2.px() *P_B2.px() +P_B2.py() *P_B2.py() +P_B2.pz() *P_B2.pz() );
240
241 double sintheta1 = sqrt(1-costheta1*costheta1);
242 double sintheta2 = sqrt(1-costheta2*costheta2);
243
244 // Cosine of hard scattering
245 double costhe = (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
246
247 // Invariant mass of tau+tau- pair
248 double SS = S; // other option is for tests: P_QQ.recalculated_mass()*P_QQ.recalculated_mass();
249
250
251 // corrX2 calculation
252 // polX2 calculation
253 // WTnonSM calculation
254 *corrX2 = -1.0; // default
255 *polX2 = 0.0; // default
256 WTnonSM = 1.0; // default
257}
258
259} // namespace TauSpinner
double plweight(int ide, double svar, double costhe)
Definition nonSM.cxx:175
double plzap2(int ide, int idf, double svar, double costhe)
Definition nonSM.cxx:115
double plnorm(int ide, double svar)
Definition nonSM.cxx:189
double default_nonSM_bornH(int ID, double S, double cost, int H1, int H2, int key)
Definition nonSM.cxx:90
void nonSMHcorrPol(double S, SimpleParticle &tau1, SimpleParticle &tau2, double *corrX2, double *polX2)
Definition nonSM.cxx:216
int set_nonSM_bornH(double(*fun)(int, double, double, int, int, int))
Definition nonSM.cxx:106
int set_nonSM_born(double(*fun)(int, double, double, int, int, int))
Definition nonSM.cxx:100
int initEWff(int ID, double S, double cost, int key)
Definition EWtables.cxx:563
double default_nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
Definition nonSM.cxx:27
double nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
Definition nonSM.cxx:20