C++InterfacetoTauola
nonSM.cxx
1 #include <iostream>
2 #include "TauSpinner/nonSM.h"
3 #include "TauSpinner/EWtables.h"
4 using std::cout;
5 using std::endl;
6 
7 namespace TauSpinner {
8 
9 /***** GLOBAL VARIABLES *****/
10 
11 double (*nonSM_bornZ)(int, double, double, int, int, int) = default_nonSM_born;
12 double (*nonSM_bornH)(int, double, double, int, int, int) = default_nonSM_bornH;
13 
14 extern int nonSM2;
15 extern int nonSMN;
16 extern double WTnonSM;
17 extern bool IfHiggs;
18 /********************************/
19 
20 double 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 
27 double 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 
90 double 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 
100 int 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 
106 int 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 *******************************************************************************/
115 double 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 
175 double 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 
189 double 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 
216 void 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 default_nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
Definition: nonSM.cxx:27
double plzap2(int ide, int idf, double svar, double costhe)
Definition: nonSM.cxx:115
int set_nonSM_bornH(double(*fun)(int, double, double, int, int, int))
Definition: nonSM.cxx:106
double default_nonSM_bornH(int ID, double S, double cost, int H1, int H2, int key)
Definition: nonSM.cxx:90
double nonSM_born(int ID, double S, double cost, int H1, int H2, int key)
Definition: nonSM.cxx:20
double plweight(int ide, double svar, double costhe)
Definition: nonSM.cxx:175
void nonSMHcorrPol(double S, SimpleParticle &tau1, SimpleParticle &tau2, double *corrX2, double *polX2)
Definition: nonSM.cxx:216
int initEWff(int ID, double S, double cost, int key)
Definition: EWtables.cxx:563
double plnorm(int ide, double svar)
Definition: nonSM.cxx:189
int set_nonSM_born(double(*fun)(int, double, double, int, int, int))
Definition: nonSM.cxx:100