]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //---------------------------------------------------------------------------------- |
2 | // | |
3 | // Module: EvtLb2Lll.cpp | |
4 | // | |
5 | // Desription: Routine to implement Lambda_b0 -> Lambda_0 l+ l- decays accroding to | |
6 | // several models: Chen. Geng. | |
7 | // Aliev. Ozpineci. Savci. | |
8 | // | |
9 | // Modification history: | |
10 | // | |
11 | // 09/02/2009 PR Commented check for (anti-)Lambda0 names | |
12 | // 15/09/2004 PR Module created according to PHSP model | |
13 | // 20/02/2005 PR Added parameters, created matrix element (without polarization) | |
14 | // 04/03/2005 PR LD contrib., corrected WC eff. according to Chen. Geng. | |
15 | // | |
16 | //---------------------------------------------------------------------------------- | |
17 | ||
18 | #ifdef WIN32 | |
19 | #pragma warning( disable : 4786 ) | |
20 | // Disable anoying warning about symbol size | |
21 | #endif | |
22 | ||
23 | #include "EvtGenModels/EvtLb2Lll.hh" | |
24 | #include "EvtGenBase/EvtParticle.hh" | |
25 | #include "EvtGenBase/EvtDiracParticle.hh" | |
26 | #include "EvtGenBase/EvtPDL.hh" | |
27 | #include "EvtGenBase/EvtDiracSpinor.hh" | |
28 | #include "EvtGenBase/EvtTensor4C.hh" | |
29 | #include "EvtGenBase/EvtVector4C.hh" | |
30 | #include "EvtGenBase/EvtVector4R.hh" | |
31 | #include "EvtGenBase/EvtComplex.hh" | |
32 | #include "EvtGenBase/EvtGammaMatrix.hh" | |
33 | #include <stdio.h> | |
34 | #include <string.h> | |
35 | ||
36 | EvtLb2Lll::~EvtLb2Lll() {} | |
37 | ||
38 | EvtDecayBase* EvtLb2Lll::clone(){ | |
39 | return new EvtLb2Lll; | |
40 | } | |
41 | ||
42 | std::string EvtLb2Lll::getName(){ | |
43 | return "Lb2Lll"; | |
44 | } | |
45 | ||
46 | void EvtLb2Lll::init(){ | |
47 | ||
48 | if(getNArg()>8){ // Decay parameters | |
49 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected max. 8 arguments but found: " << getNArg() << std::endl; | |
50 | report(INFO ,"EvtGen") << " 1. Lambda_b0 polarization - zero is default" << std::endl; | |
51 | report(INFO ,"EvtGen") << " 2. Model type - \"SM\" for Standard Model is default" << std::endl; | |
52 | report(INFO ,"EvtGen") << " 3. Form-Factors - \"HQET\" is used by default" << std::endl; | |
53 | report(INFO ,"EvtGen") << " 4. How to set polarization - \"ModifiedSpinors\" is default" << std::endl; | |
54 | report(INFO ,"EvtGen") << " 5. Include long distance (LD) effects - \"SD\" (no) is default" << std::endl; | |
55 | report(INFO ,"EvtGen") << " 6. NonFactorizable contribution (omega) to b->sg decay at q2=0 " << std::endl; | |
56 | report(INFO ,"EvtGen") << " 7. Note on every x-th decay" << std::endl; | |
57 | report(INFO ,"EvtGen") << " 8. Maximum probability - automatic by default" << std::endl; | |
58 | ::abort(); | |
59 | } | |
60 | ||
61 | if(getNDaug()!=3){ // Check that there are 3 daughters only | |
62 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected 3 daughters but found: " << getNDaug() << std::endl; | |
63 | ::abort(); | |
64 | } | |
65 | ||
66 | // TODO: better check based on spin and falvour is needed to allow usage of aliases ! | |
67 | if(EvtPDL::name(getParentId())=="Lambda_b0"){ // Check daughters of Lambda_b0 | |
68 | report(INFO,"EvtGen") << " EvtLb2Lll generator found Lambda_b0" << std::endl; | |
69 | //if(EvtPDL::name(getDaug(0))!="Lambda0"){ | |
70 | // report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl; | |
71 | // ::abort(); | |
72 | //} | |
73 | if(EvtPDL::name(getDaug(1))=="e-" && EvtPDL::name(getDaug(2))=="e+"){ | |
74 | m_decayName="Lambda_b0 -> Lambda0 e- e+"; | |
75 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 e- e+" << std::endl; | |
76 | }else if(EvtPDL::name(getDaug(1))=="mu-" && EvtPDL::name(getDaug(2))=="mu+"){ | |
77 | m_decayName="Lambda_b0 -> Lambda0 mu- mu+"; | |
78 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 mu- mu+" << std::endl; | |
79 | }else if(EvtPDL::name(getDaug(1))=="tau-" && EvtPDL::name(getDaug(2))=="tau+"){ | |
80 | m_decayName="Lambda_b0 -> Lambda0 tau- tau+"; | |
81 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 tau- tau+" << std::endl; | |
82 | }else{ | |
83 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found: " << EvtPDL::name(getDaug(1)) << " " << EvtPDL::name(getDaug(2)) << std::endl; | |
84 | ::abort(); | |
85 | } | |
86 | //TODO: The model is known not to work correctly for anti-Lambda_b0 (A_FB does not change its sign) | |
87 | }else if(EvtPDL::name(getParentId())=="anti-Lambda_b0"){ // Check daughters of anti-Lambda_b0 | |
88 | report(INFO,"EvtGen") << " EvtLb2Lll generator found anti-Lambda_b0" << std::endl; | |
89 | //if(EvtPDL::name(getDaug(0))!="anti-Lambda0"){ | |
90 | // report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected anti-Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl; | |
91 | // ::abort(); | |
92 | //} | |
93 | if(EvtPDL::name(getDaug(1))=="e+" && EvtPDL::name(getDaug(2))=="e-"){ | |
94 | m_decayName="anti-Lambda_b0 -> anti-Lambda0 e+ e-"; | |
95 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 e+ e-" << std::endl; | |
96 | }else if(EvtPDL::name(getDaug(1))=="mu+" && EvtPDL::name(getDaug(2))=="mu-"){ | |
97 | m_decayName="anti-Lambda_b0 -> anti-Lambda0 mu+ mu-"; | |
98 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 mu+ mu-" << std::endl; | |
99 | }else if(EvtPDL::name(getDaug(1))=="tau-" && EvtPDL::name(getDaug(2))=="tau+"){ | |
100 | m_decayName="anti-Lambda_b0 -> anti-Lambda0 tau+ tau-"; | |
101 | report(INFO,"EvtGen") << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 tau+ tau-" << std::endl; | |
102 | }else{ | |
103 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found: " << EvtPDL::name(getDaug(1)) << " " << EvtPDL::name(getDaug(2)) << std::endl; | |
104 | ::abort(); | |
105 | } | |
106 | }else{ // This model is not intended for decay of anything else than (anti-)Lambda_b0 | |
107 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected (anti-)Lambda_b0 parent but found: " << EvtPDL::name(getParentId()) << std::endl; | |
108 | ::abort(); | |
109 | } | |
110 | ||
111 | // Read and check all parameters | |
112 | if(getNArg()>0){ | |
113 | if(getArg(0)>1. || getArg(0)<-1.){ | |
114 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll expects polarization to be in interval <-1,1>, not " << getArg(0) << std::endl; | |
115 | ::abort(); | |
116 | } | |
117 | m_polarizationLambdab0 = getArg(0); | |
118 | }else{ | |
119 | m_polarizationLambdab0 = 0; | |
120 | } | |
121 | report(INFO,"EvtGen") << " EvtLb2Lll set Lambda_b0 polarization to " << m_polarizationLambdab0 << std::endl; | |
122 | ||
123 | if(getNArg()>1){ | |
124 | if(getArgStr(1).substr(1,getArgStr(1).size()-2)!="SM" && | |
125 | getArgStr(1).substr(1,getArgStr(1).size()-2)!="-C7_SM" && | |
126 | getArgStr(1).substr(1,getArgStr(1).size()-2)!="SUSY-ChenGeng"){ | |
127 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know this physics model: " << getArgStr(1) << std::endl; | |
128 | ::abort(); | |
129 | } | |
130 | m_HEPmodel = getArgStr(1).substr(1,getArgStr(1).size()-2); | |
131 | }else{ | |
132 | m_HEPmodel = "SM"; | |
133 | } | |
134 | report(INFO,"EvtGen") << " EvtLb2Lll will use this physics model: " << m_HEPmodel << std::endl; | |
135 | ||
136 | if(getNArg()>2){ | |
137 | if(getArgStr(2).substr(1,getArgStr(2).size()-2)!="HQET" && | |
138 | getArgStr(2).substr(1,getArgStr(2).size()-2)!="HQET-noF2" && | |
139 | getArgStr(2).substr(1,11) !="HQET-delta="){ | |
140 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know this Form-Factors model: " << getArgStr(2) << std::endl; | |
141 | ::abort(); | |
142 | } | |
143 | m_FFtype = getArgStr(2).substr(1,getArgStr(2).size()-2); | |
144 | }else{ | |
145 | m_FFtype = "HQET"; | |
146 | } | |
147 | report(INFO,"EvtGen") << " EvtLb2Lll will use this Form-Factors model: " << m_FFtype << std::endl; | |
148 | ||
149 | if(getNArg()>3){ | |
150 | if(getArgStr(3).substr(1,getArgStr(3).size()-2)!="Unpolarized"){ | |
151 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll doesn't know kind of introducing polarization: " << getArgStr(3) << std::endl; | |
152 | ::abort(); | |
153 | } | |
154 | m_polarizationIntroduction = getArgStr(3).substr(1,getArgStr(3).size()-2); | |
155 | }else{ | |
156 | m_polarizationIntroduction = "Unpolarized"; | |
157 | } | |
158 | report(INFO,"EvtGen") << " EvtLb2Lll will use this kind of introducing polarization: " << m_polarizationIntroduction << std::endl; | |
159 | ||
160 | if(getNArg()>4){ | |
161 | if(getArgStr(4).substr(1,getArgStr(4).size()-2)!="SD" && getArgStr(4).substr(1,getArgStr(4).size()-2)!="LD"){ | |
162 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll didn't find SD or LD parameter: " << getArgStr(4) << std::endl; | |
163 | ::abort(); | |
164 | } | |
165 | m_effectContribution = getArgStr(5).substr(1,getArgStr(4).size()-2); | |
166 | }else{ | |
167 | m_effectContribution = "SD"; | |
168 | } | |
169 | report(INFO,"EvtGen") << " EvtLb2Lll will include contribution from these effects: " << m_effectContribution << std::endl; | |
170 | ||
171 | if(getNArg()>5){ | |
172 | if(fabs(getArg(5))>0.15){ | |
173 | report(WARNING,"EvtGen") << " WARNING: EvtLb2Lll found very high contribution to b->sg decay at q2=0: " << getArg(5) << std::endl; | |
174 | } | |
175 | m_omega = getArg(5); | |
176 | }else{ | |
177 | m_omega = 0; | |
178 | } | |
179 | report(INFO,"EvtGen") << " EvtLb2Lll will use this contribution to b->sg decay at q2=0: " << m_omega << std::endl; | |
180 | ||
181 | if(getNArg()>6) m_noTries = (long)(getArg(6)); | |
182 | else m_noTries = 0; | |
183 | ||
184 | if(getNArg()>7){ | |
185 | if(getArg(7)<0.){ | |
186 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll expects positive maximum probability not : " << getArg(7) << std::endl; | |
187 | ::abort(); | |
188 | } | |
189 | m_maxProbability = getArg(7); | |
190 | }else{ | |
191 | m_maxProbability = 0.; | |
192 | } | |
193 | report(INFO,"EvtGen") << " EvtLb2Lll maximum probability was set to " << m_maxProbability << std::endl; | |
194 | m_poleSize=0; | |
195 | ||
196 | // Initialize Wilson coeficients by Buras and Munz | |
197 | // TODO: should have common W.C. source for all decays in EvtGen | |
198 | m_WC.CalculateAllCoeficients(); | |
199 | ||
200 | } | |
201 | ||
202 | void EvtLb2Lll::initProbMax(){ | |
203 | ||
204 | report(INFO,"EvtGen") << " EvtLb2Lll is finding maximum probability ... " << std::endl; | |
205 | ||
206 | if(m_maxProbability==0){ | |
207 | ||
208 | EvtDiracParticle *parent = new EvtDiracParticle; | |
209 | parent->noLifeTime(); | |
210 | parent->init(getParentId(),EvtVector4R(EvtPDL::getMass(getParentId()),0,0,0)); | |
211 | parent->setDiagonalSpinDensity(); | |
212 | ||
213 | EvtAmp amp; | |
214 | EvtId daughters[3] = {getDaug(0),getDaug(1),getDaug(2)}; | |
215 | amp.init(getParentId(),3,daughters); | |
216 | parent->makeDaughters(3,daughters); | |
217 | EvtParticle *lambda = parent->getDaug(0); | |
218 | EvtParticle *lep1 = parent->getDaug(1); | |
219 | EvtParticle *lep2 = parent->getDaug(2); | |
220 | lambda -> noLifeTime(); | |
221 | lep1 -> noLifeTime(); | |
222 | lep2 -> noLifeTime(); | |
223 | ||
224 | EvtSpinDensity rho; | |
225 | rho.setDiag(parent->getSpinStates()); | |
226 | ||
227 | double M0 = EvtPDL::getMass(getParentId()); | |
228 | double mL = EvtPDL::getMass(getDaug(0)); | |
229 | double m1 = EvtPDL::getMass(getDaug(1)); | |
230 | double m2 = EvtPDL::getMass(getDaug(2)); | |
231 | ||
232 | double q2,pstar,elambda,theta; | |
233 | double q2min = (m1+m2)*(m1+m2); | |
234 | double q2max = (M0-mL)*(M0-mL); | |
235 | ||
236 | EvtVector4R p4lambda,p4lep1,p4lep2,boost; | |
237 | ||
238 | report(INFO,"EvtGen") << " EvtLb2Lll is probing whole phase space ..." << std::endl; | |
239 | ||
240 | int i,j; | |
241 | double prob=0; | |
242 | for(i=0;i<=100;i++){ | |
243 | q2 = q2min+i*(q2max-q2min)/100.; | |
244 | elambda = (M0*M0+mL*mL-q2)/2/M0; | |
245 | if(i==0) pstar = 0; | |
246 | else pstar = sqrt(q2-(m1+m2)*(m1+m2))*sqrt(q2-(m1-m2)*(m1-m2))/2/sqrt(q2); | |
247 | boost.set(M0-elambda,0,0,+sqrt(elambda*elambda-mL*mL)); | |
248 | p4lambda.set(elambda,0,0,-sqrt(elambda*elambda-mL*mL)); | |
249 | for(j=0;j<=45;j++){ | |
250 | theta = j*EvtConst::pi/45; | |
251 | p4lep1.set(sqrt(pstar*pstar+m1*m1),0,+pstar*sin(theta),+pstar*cos(theta)); | |
252 | p4lep2.set(sqrt(pstar*pstar+m2*m2),0,-pstar*sin(theta),-pstar*cos(theta)); | |
253 | //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " pstar: " << pstar << std::endl; | |
254 | p4lep1 = boostTo(p4lep1,boost); | |
255 | p4lep2 = boostTo(p4lep1,boost); | |
256 | lambda -> init(getDaug(0),p4lambda); | |
257 | lep1 -> init(getDaug(1),p4lep1 ); | |
258 | lep2 -> init(getDaug(2),p4lep2 ); | |
259 | calcAmp(&,parent); | |
260 | prob = rho.normalizedProb(amp.getSpinDensity()); | |
261 | //std::cout << "q2: " << q2 << " \t theta: " << theta << " \t prob: " << prob << std::endl; | |
262 | //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl; | |
263 | if(prob>m_maxProbability){ | |
264 | report(INFO,"EvtGen") << " - probability " << prob << " found at q2 = " << q2 << " (" << 100*(q2-q2min)/(q2max-q2min) << " %) and theta = " << theta*180/EvtConst::pi << std::endl; | |
265 | m_maxProbability=prob; | |
266 | } | |
267 | } | |
268 | //::abort(); | |
269 | } | |
270 | ||
271 | //m_poleSize = 0.04*q2min; | |
272 | m_maxProbability *= 1.2; | |
273 | delete parent; | |
274 | } | |
275 | ||
276 | setProbMax(m_maxProbability); | |
277 | report(INFO,"EvtGen") << " EvtLb2Lll set up maximum probability to " << m_maxProbability << std::endl; | |
278 | ||
279 | } | |
280 | ||
281 | void EvtLb2Lll::decay(EvtParticle* parent){ | |
282 | ||
283 | //setWeight(parent->initializePhaseSpace(getNDaug(),getDaugs(),m_poleSize,1,2)); | |
284 | parent->initializePhaseSpace(getNDaug(),getDaugs()); | |
285 | calcAmp(&_amp2,parent); | |
286 | ||
287 | } | |
288 | ||
289 | void EvtLb2Lll::calcAmp(EvtAmp *amp,EvtParticle *parent){ | |
290 | ||
291 | static long noTries=0; | |
292 | static double delta=0; | |
293 | ||
294 | EvtComplex Matrix[2][2][2][2]; | |
295 | ||
296 | EvtComplex i1(0,1); | |
297 | ||
298 | int i,j,spins[4]; | |
299 | char ch; | |
300 | ||
301 | double r,M_L,M_Lb,M_s,M_c,M_b,q2,alpha,alpha_s,eta,M_W,M_t; | |
302 | double M_psi[2]={0,0},Gamma_psi[2]={0,0},k_psi[2]={0,0}; | |
303 | double F0_1,F0_2,a_F1,a_F2,b_F1,b_F2,F1,F2; | |
304 | double f_1,f_2,f_3,g_1,g_2,g_3,f_1T,f_2T,f_3T,g_1T,g_2T,g_3T,f_TV,f_TS,g_TV,g_TS,f_T,g_T; | |
305 | EvtComplex A1,A2,A3,B1,B2,B3,D1,D2,D3,E1,E2,E3,N1,N2,H1,H2; | |
306 | EvtComplex C_SL,C_BR,C_LLtot,C_LRtot,C_LL,C_LR,C_RL,C_RR,C_LRLR,C_RLLR,C_LRRL,C_RLRL,C_T,C_TE; | |
307 | EvtComplex Yld,C_7eff,C_9eff; | |
308 | EvtComplex V_ts,V_tb; | |
309 | ||
310 | EvtVector4C lbar_Gmu_l[2][2],lbar_GmuG5_l[2][2],hbar_GmuPlusG5_h[2][2],hbar_GmuMinusG5_h[2][2],hbar_Gmu_h[2][2]; | |
311 | EvtComplex lbar_l[2][2],lbar_G5_l[2][2],hbar_1PlusG5_h[2][2],hbar_1MinusG5_h[2][2],hbar_G5_h[2][2],hbar_h[2][2]; | |
312 | EvtTensor4C lbar_Smunu_l[2][2],lbar_ESmunu_l[2][2],hbar_SmunuPlusG5_h[2][2],hbar_SmunuMinusG5_h[2][2],hbar_Smunu_h[2][2]; | |
313 | EvtVector4R q_mu,P_mu; | |
314 | ||
315 | EvtDiracSpinor parent__spParent[2]; | |
316 | ||
317 | M_Lb = parent->mass(); | |
318 | M_L = parent->getDaug(0)->mass(); | |
319 | M_s = 0.13; | |
320 | M_c = 1.35; | |
321 | M_b = 4.8; | |
322 | alpha = 1./137.036; | |
323 | alpha_s = 0.217484; | |
324 | eta = 0.556289; | |
325 | M_W = 80.425; | |
326 | M_t = 174.3; | |
327 | M_psi[0] = 3.096916; | |
328 | M_psi[1] = 3.686093; | |
329 | if(m_decayName=="Lambda_b0 -> Lambda0 e- e+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 e+ e-"){ | |
330 | Gamma_psi[0] = 5.40; | |
331 | Gamma_psi[1] = 2.12; | |
332 | } | |
333 | if(m_decayName=="Lambda_b0 -> Lambda0 mu- mu+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 mu+ mu-"){ | |
334 | Gamma_psi[0] = 5.35; | |
335 | Gamma_psi[1] = 2.05; | |
336 | } | |
337 | if(m_decayName=="Lambda_b0 -> Lambda0 tau- tau+" || m_decayName=="anti-Lambda_b0 -> anti-Lambda0 tau+ tau-"){ | |
338 | Gamma_psi[0] = 0.00; | |
339 | Gamma_psi[1] = 0.79; | |
340 | } | |
341 | if(m_effectContribution=="LD"){ | |
342 | k_psi[0] = 1.65; | |
343 | k_psi[1] = 1.65; | |
344 | } | |
345 | //G_F = 1.16637e-5; | |
346 | //V_tb = sqrt(1-pow(0.0413,2))*sqrt(1-pow(0.0037,2)); | |
347 | //V_ts = -sqrt(1-pow(0.2243,2))*0.0413-0.2243*sqrt(1-pow(0.0413,2))*0.0037*(cos(60*EvtConst::pi/180)+i1*sin(60*EvtConst::pi/180)); | |
348 | ||
349 | P_mu = parent->getP4Restframe()+parent->getDaug(0)->getP4(); | |
350 | q_mu = parent->getP4Restframe()-parent->getDaug(0)->getP4(); | |
351 | q2 = q_mu.mass2(); | |
352 | ||
353 | if(m_noTries>0) if(!((++noTries)%m_noTries)) report(DEBUG,"EvtGen") << " EvtLb2Lll already finished " << noTries << " matrix element calculations" << std::endl; | |
354 | ||
355 | if(m_FFtype=="HQET"){ | |
356 | r = M_L*M_L/M_Lb/M_Lb; | |
357 | F0_1 = +0.462; | |
358 | F0_2 = -0.077; | |
359 | a_F1 = -0.0182; | |
360 | a_F2 = -0.0685; | |
361 | b_F1 = -0.000176; | |
362 | b_F2 = +0.001460; | |
363 | F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb))); | |
364 | F2 = F0_2/(1.0-(q2/M_Lb/M_Lb)*(a_F2-b_F2*(q2/M_Lb/M_Lb))); | |
365 | g_1 = f_1 = f_2T = g_2T = F1+sqrt(r)*F2; | |
366 | //std::cout << " F1: " << F1 << " F2: " << F2 << " r: " << r << " M_L: " << M_L << " M_Lb: " << M_Lb << std::endl; | |
367 | //std::cout << " sqrt(q2): " << sqrt(q2) << " q2: " << q2 << " M_Lb^2" << M_Lb*M_Lb << std::endl; | |
368 | g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2/M_Lb; | |
369 | g_TS = f_TS = 0; | |
370 | g_1T = f_1T = F2/M_Lb*q2; | |
371 | g_3T = +F2/M_Lb*(M_Lb+M_L); | |
372 | f_3T = -F2/M_Lb*(M_Lb-M_L); | |
373 | f_T = f_2T-f_TS*q2; | |
374 | g_T = g_2T-g_TS*q2; | |
375 | }else if(strstr(m_FFtype.c_str(),"HQET-delta=")==m_FFtype.c_str()){ | |
376 | //report(WARNING,"EvtGen") << " WARNING: HQET-delta FF model should be checked for correctness" << std::endl; | |
377 | if(delta==0) sscanf(m_FFtype.c_str(),"%c%c%c%c%c%c%c%c%c%c%c%lf",&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&ch,&delta); | |
378 | r = M_L*M_L/M_Lb/M_Lb; | |
379 | F0_1 = +0.462; | |
380 | F0_2 = -0.077; | |
381 | a_F1 = -0.0182; | |
382 | a_F2 = -0.0685; | |
383 | b_F1 = -0.000176; | |
384 | b_F2 = +0.001460; | |
385 | F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb))); | |
386 | F2 = F0_2/(1.0-(q2/M_Lb/M_Lb)*(a_F2-b_F2*(q2/M_Lb/M_Lb))); | |
387 | g_1 = f_1 = f_2T = g_2T = F1+sqrt(r)*F2; | |
388 | g_1 += delta*g_1; | |
389 | f_1 -= delta*f_1; | |
390 | g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2/M_Lb; | |
391 | g_TS = f_TS = 0; | |
392 | g_1T = f_1T = F2/M_Lb*q2; | |
393 | g_3T = +F2/M_Lb*(M_Lb+M_L); | |
394 | f_3T = -F2/M_Lb*(M_Lb-M_L); | |
395 | f_T = f_2T-f_TS*q2; | |
396 | g_T = g_2T-g_TS*q2; | |
397 | }else if(m_FFtype=="HQET-noF2"){ | |
398 | //report(WARNING,"EvtGen") << " WARNING: HQET-noF2 FF model should be checked for correctness" << std::endl; | |
399 | r = M_L*M_L/M_Lb/M_Lb; | |
400 | F0_1 = +0.462; | |
401 | a_F1 = -0.0182; | |
402 | b_F1 = -0.000176; | |
403 | F1 = F0_1/(1.0-(q2/M_Lb/M_Lb)*(a_F1-b_F1*(q2/M_Lb/M_Lb))); | |
404 | g_1 = f_1 = f_2T = g_2T = F1; | |
405 | g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = 0; | |
406 | g_TS = f_TS = 0; | |
407 | g_1T = f_1T = 0; | |
408 | g_3T = 0; | |
409 | f_3T = 0; | |
410 | f_T = f_2T-f_TS*q2; | |
411 | g_T = g_2T-g_TS*q2; | |
412 | }else{ // general relations for Form-Factors | |
413 | f_1 = f_2 = f_3 = g_1 = g_2 = g_3 = f_3T = g_3T = f_TS = g_TS = f_T = g_T = f_TV = 0; | |
414 | f_2T = f_T+f_TS*q2; | |
415 | f_1T = (f_TV+f_TS*(M_L+M_Lb))*q2; | |
416 | f_1T = -q2/(M_Lb-M_L)*f_3T; | |
417 | g_2T = g_T+g_TS*q2; | |
418 | g_1T = (g_TV-g_TS*(M_L-M_Lb))*q2; | |
419 | g_1T = +q2/(M_Lb+M_L)*g_3T; | |
420 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown Form-Factors model: " << m_FFtype << " - this should never happen !" << std::endl; | |
421 | ::abort(); | |
422 | } | |
423 | ||
424 | if(m_HEPmodel=="SM"){ | |
425 | C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0); | |
426 | Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha); | |
427 | C_7eff = m_WC.GetC7eff0() + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld); | |
428 | C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda(),m_WC.GetRenormSchemePar()); | |
429 | C_SL = -2*M_s*C_7eff; | |
430 | C_BR = -2*M_b*C_7eff; | |
431 | C_LLtot = C_9eff-m_WC.GetC10tilda()+C_LL; | |
432 | C_LRtot = C_9eff+m_WC.GetC10tilda()+C_LR; | |
433 | //std::cout << "Yld: " << Yld << " C7eff: " << C_7eff << " C_9eff: " << C_9eff << " Diff7: " << C_7eff-m_WC.GetC7eff0() << " Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl; | |
434 | }else if(m_HEPmodel=="-C7_SM"){ | |
435 | C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0); | |
436 | Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha); | |
437 | C_7eff = m_WC.GetC7eff0() + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld); | |
438 | C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda(),m_WC.GetRenormSchemePar()); | |
439 | C_SL = +2*M_s*C_7eff; | |
440 | C_BR = +2*M_b*C_7eff; | |
441 | C_LLtot = C_9eff-m_WC.GetC10tilda()+C_LL; | |
442 | C_LRtot = C_9eff+m_WC.GetC10tilda()+C_LR; | |
443 | //std::cout << "Yld: " << Yld << " C7eff: " << C_7eff << " C_9eff: " << C_9eff << " Diff7: " << C_7eff-m_WC.GetC7eff0() << " Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl; | |
444 | }else if(m_HEPmodel=="SUSY-ChenGeng"){ | |
445 | //report(WARNING,"EvtGen") << " WARNING: SUSY-ChenGeng model should be checked for correctness" << std::endl; | |
446 | C_LL=C_LR=C_RL=C_RR=C_LRLR=C_RLLR=C_LRRL=C_RLRL=C_T=C_TE=EvtComplex(0,0); | |
447 | EvtComplex d_u23LL = 0.1; | |
448 | EvtComplex d_u33RL = 0.65; | |
449 | EvtComplex d_d23LR = 0.03*exp(i1*EvtConst::pi*2/5); | |
450 | EvtComplex d_u23LR = -0.8*exp(i1*EvtConst::pi/4); | |
451 | EvtComplex C_7susy = -1.75*d_u23LL - 0.25*d_u23LR - 10.3*d_d23LR; | |
452 | EvtComplex C_9susy = 0.82*d_u23LR; | |
453 | EvtComplex C_10susy = -9.37*d_u23LR + 1.4*d_u23LR*d_u33RL + 2.7*d_u23LL; | |
454 | Yld = m_WC.Yld(q2,k_psi,Gamma_psi,M_psi,2,m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),1./alpha); | |
455 | C_7eff = m_WC.GetC7eff0() + C_7susy*pow(m_WC.GetEta(),16./23.) + m_WC.C7b2sg(m_WC.GetStrongCouplingConst(),m_WC.GetEta(),m_WC.GetC2(),M_t,M_W) + m_omega*(m_WC.hzs(M_c/M_b,q2/M_b/M_b,M_b,M_b) + Yld); | |
456 | C_9eff = Yld + m_WC.C9efftilda(M_c/M_b,q2/M_b/M_b,m_WC.GetStrongCouplingConst(),m_WC.GetC1(),m_WC.GetC2(),m_WC.GetC3(),m_WC.GetC4(),m_WC.GetC5(),m_WC.GetC6(),m_WC.GetC9tilda()+C_9susy,m_WC.GetRenormSchemePar()); | |
457 | C_SL = -2*M_s*C_7eff; | |
458 | C_BR = -2*M_b*C_7eff; | |
459 | C_LLtot = C_9eff-m_WC.GetC10tilda()-C_10susy+C_LL; | |
460 | C_LRtot = C_9eff+m_WC.GetC10tilda()+C_10susy+C_LR; | |
461 | //std::cout << "Yld: " << Yld << " C7eff: " << C_7eff << " C_9eff: " << C_9eff << " Diff7: " << C_7eff-m_WC.GetC7eff0() << " Diff9: " << C_9eff-m_WC.GetC9tilda() << std::endl; | |
462 | }else{ | |
463 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown physics model: " << m_HEPmodel << " - this should never happen !" << std::endl; | |
464 | ::abort(); | |
465 | } | |
466 | ||
467 | A1 = (f_1T-g_1T)*C_SL/q2 + (f_1T+g_1T)*C_BR/q2 + 0.5*(f_1-g_1)*(C_LLtot+C_LRtot) + 0.5*(f_1+g_1)*(C_RL+C_RR); | |
468 | //std::cout << "f_1T: " << f_1T << " g_1T: " << g_1T << " C_SL: " << C_SL << " C_BR: " << C_BR << " f_1: " << f_1 << " g_1: " << g_1 << " C_LLtot: " << C_LLtot << " C_LRtot: " << C_LRtot << " C_RL: " << C_RL << " C_RR: " << C_RR << std::endl; | |
469 | A2 = (f_2T-g_2T)*C_SL/q2 + (f_2T+g_2T)*C_BR/q2 + 0.5*(f_2-g_2)*(C_LLtot+C_LRtot) + 0.5*(f_2+g_2)*(C_RL+C_RR); | |
470 | A3 = (f_3T-g_3T)*C_SL/q2 + (f_3T+g_3T)*C_BR/q2 + 0.5*(f_3-g_3)*(C_LLtot+C_LRtot) + 0.5*(f_3+g_3)*(C_RL+C_RR); | |
471 | ||
472 | B1 = (f_1T+g_1T)*C_SL/q2 + (f_1T-g_1T)*C_BR/q2 + 0.5*(f_1+g_1)*(C_LLtot+C_LRtot) + 0.5*(f_1-g_1)*(C_RL+C_RR); | |
473 | B2 = (f_2T+g_2T)*C_SL/q2 + (f_2T-g_2T)*C_BR/q2 + 0.5*(f_2+g_2)*(C_LLtot+C_LRtot) + 0.5*(f_2-g_2)*(C_RL+C_RR); | |
474 | B3 = (f_3T+g_3T)*C_SL/q2 + (f_3T-g_3T)*C_BR/q2 + 0.5*(f_3+g_3)*(C_LLtot+C_LRtot) + 0.5*(f_3-g_3)*(C_RL+C_RR); | |
475 | ||
476 | D1 = 0.5*(C_RR-C_RL)*(f_1+g_1) + 0.5*(C_LRtot-C_LLtot)*(f_1-g_1); | |
477 | D2 = 0.5*(C_RR-C_RL)*(f_2+g_2) + 0.5*(C_LRtot-C_LLtot)*(f_2-g_2); | |
478 | D3 = 0.5*(C_RR-C_RL)*(f_3+g_3) + 0.5*(C_LRtot-C_LLtot)*(f_3-g_3); | |
479 | ||
480 | E1 = 0.5*(C_RR-C_RL)*(f_1-g_1) + 0.5*(C_LRtot-C_LLtot)*(f_1+g_1); | |
481 | E2 = 0.5*(C_RR-C_RL)*(f_2-g_2) + 0.5*(C_LRtot-C_LLtot)*(f_2+g_2); | |
482 | E3 = 0.5*(C_RR-C_RL)*(f_3-g_3) + 0.5*(C_LRtot-C_LLtot)*(f_3+g_3); | |
483 | ||
484 | N1 = (f_1*(M_Lb+M_L)+f_3*q2)/M_b*(C_LRLR+C_RLLR+C_LRRL+C_RLRL); | |
485 | N2 = (f_1*(M_Lb+M_L)+f_3*q2)/M_b*(C_LRLR+C_RLLR-C_LRRL-C_RLRL); | |
486 | ||
487 | H1 = (g_1*(M_Lb+M_L)-g_3*q2)/M_b*(C_LRLR-C_RLLR+C_LRRL-C_RLRL); | |
488 | H2 = (g_1*(M_Lb+M_L)-g_3*q2)/M_b*(C_LRLR-C_RLLR-C_LRRL+C_RLRL); | |
489 | ||
490 | for(i=0;i<4;i++){ | |
491 | lbar_Gmu_l [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)); | |
492 | lbar_GmuG5_l [i/2][i%2] = EvtLeptonACurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)); | |
493 | lbar_l [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)); | |
494 | lbar_G5_l [i/2][i%2] = EvtLeptonPCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)); | |
495 | lbar_Smunu_l [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2)); | |
496 | lbar_ESmunu_l[i/2][i%2] = dual(EvtLeptonTCurrent(parent->getDaug(1)->spParent(i/2),parent->getDaug(2)->spParent(i%2))); | |
497 | } | |
498 | ||
499 | // TODO: polarization not yet introduced | |
500 | if(m_polarizationIntroduction=="SpinDensityMatrix"){ | |
501 | //parent->setSpinDensityForward(); | |
502 | parent__spParent[0]=parent->sp(0); | |
503 | parent__spParent[1]=parent->sp(1); | |
504 | }else if(m_polarizationIntroduction=="ModifiedSpinors"){ | |
505 | parent__spParent[0]=parent->sp(0); | |
506 | parent__spParent[1]=parent->sp(1); | |
507 | }else if(m_polarizationIntroduction=="Unpolarized"){ | |
508 | parent__spParent[0]=parent->sp(0); | |
509 | parent__spParent[1]=parent->sp(1); | |
510 | }else{ | |
511 | report(ERROR,"EvtGen") << " ERROR: EvtLb2Lll - unknown polarization: " << m_polarizationIntroduction << " - this should never happen !" << std::endl; | |
512 | ::abort(); | |
513 | } | |
514 | ||
515 | for(i=0;i<4;i++){ | |
516 | hbar_GmuPlusG5_h [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonACurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
517 | hbar_GmuMinusG5_h [i/2][i%2] = EvtLeptonVACurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
518 | hbar_SmunuPlusG5_h [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonTG5Current(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
519 | hbar_SmunuMinusG5_h[i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) - EvtLeptonTG5Current(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
520 | hbar_1PlusG5_h [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) + EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
521 | hbar_1MinusG5_h [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]) - EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
522 | hbar_G5_h [i/2][i%2] = EvtLeptonPCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
523 | hbar_h [i/2][i%2] = EvtLeptonSCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
524 | hbar_Smunu_h [i/2][i%2] = EvtLeptonTCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
525 | hbar_Gmu_h [i/2][i%2] = EvtLeptonVCurrent(parent->getDaug(0)->spParent(i/2),parent__spParent[i%2]); | |
526 | } | |
527 | ||
528 | for(i=0;i<4;i++) for(j=0;j<4;j++) { | |
529 | //std::cout << "--------------------------------------------------" << std::endl; | |
530 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
531 | Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2] * (A1*hbar_GmuPlusG5_h[j/2][j%2]+B1*hbar_GmuMinusG5_h[j/2][j%2]); | |
532 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
533 | //std::cout << "A1: " << A1 << " B1: " << B1 << " lbar_Gmu_l: " << lbar_Gmu_l[i/2][i%2] << | |
534 | // " hbar_GmuPlusG5_h: " << hbar_GmuPlusG5_h[j/2][j%2] << " hbar_GmuMinusG5_h: " << hbar_GmuMinusG5_h[j/2][j%2] << | |
535 | // " sp1: " << parent->getDaug(1)->spParent(i/2) << " sp2: " << parent->getDaug(1)->spParent(i%2) << std::endl; | |
536 | Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2] * (i1*A2*(hbar_SmunuPlusG5_h[j/2][j%2].cont2(q_mu))+B2*(hbar_SmunuMinusG5_h[j/2][j%2].cont2(q_mu))); | |
537 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
538 | Matrix[j/2][j%2][i/2][i%2] += lbar_Gmu_l[i/2][i%2] * ((A3*hbar_1PlusG5_h[j/2][j%2]+B3*hbar_1MinusG5_h[j/2][j%2])*q_mu); | |
539 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
540 | Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2] * (D1*hbar_GmuPlusG5_h[j/2][j%2]+E1*hbar_GmuMinusG5_h[j/2][j%2]); | |
541 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
542 | Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2] * (i1*D2*(hbar_SmunuPlusG5_h[j/2][j%2].cont2(q_mu))+E2*(hbar_SmunuMinusG5_h[j/2][j%2].cont2(q_mu))); | |
543 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
544 | Matrix[j/2][j%2][i/2][i%2] += lbar_GmuG5_l[i/2][i%2] * ((D3*hbar_1PlusG5_h[j/2][j%2]+E3*hbar_1MinusG5_h[j/2][j%2])*q_mu); | |
545 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
546 | Matrix[j/2][j%2][i/2][i%2] += lbar_l[i/2][i%2] * (N1*hbar_h[j/2][j%2]+H1*hbar_G5_h[j/2][j%2]); | |
547 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
548 | Matrix[j/2][j%2][i/2][i%2] += lbar_G5_l[i/2][i%2] * (N2*hbar_h[j/2][j%2]+H2*hbar_G5_h[j/2][j%2]); | |
549 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
550 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_Smunu_l[i/2][i%2] , 4*C_T*f_T*hbar_Smunu_h[j/2][j%2]); | |
551 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
552 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_Smunu_l[i/2][i%2] , -4*C_T*f_TV*i1*(directProd(q_mu,hbar_Gmu_h[j/2][j%2])-directProd(hbar_Gmu_h[j/2][j%2],q_mu))); | |
553 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
554 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_Smunu_l[i/2][i%2] , -4*C_T*f_TS*i1*(directProd(P_mu,q_mu)-directProd(q_mu,P_mu))*hbar_h[j/2][j%2]); | |
555 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
556 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_T*i1*hbar_Smunu_h[j/2][j%2]); | |
557 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
558 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_TV*(directProd(q_mu,hbar_Gmu_h[j/2][j%2])-directProd(hbar_Gmu_h[j/2][j%2],q_mu))); | |
559 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
560 | Matrix[j/2][j%2][i/2][i%2] += cont(lbar_ESmunu_l[i/2][i%2] , 4*C_TE*f_TS*(directProd(P_mu,q_mu)-directProd(q_mu,P_mu))*hbar_h[j/2][j%2]); | |
561 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
562 | //Matrix[j/2][j%2][i/2][i%2] *= G_F*alpha/4/sqrt(2)/EvtConst::pi*V_tb*conj(V_ts); | |
563 | //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl; | |
564 | //std::cout << "--------------------------------------------------" << std::endl; | |
565 | spins[0]=j%2; | |
566 | spins[1]=j/2; | |
567 | spins[2]=i/2; | |
568 | spins[3]=i%2; | |
569 | amp->vertex(spins,Matrix[j/2][j%2][i/2][i%2]); | |
570 | } | |
571 | ||
572 | //std::cout << "==================================================" << std::endl; | |
573 | //std::cout << "Lambda_b0: " << parent->getP4Restframe() << std::endl; | |
574 | //std::cout << "Lambda0: " << parent->getDaug(0)->getP4() << std::endl; | |
575 | //std::cout << "mu-: " << parent->getDaug(1)->getP4() << std::endl; | |
576 | //std::cout << "mu+: " << parent->getDaug(2)->getP4() << std::endl; | |
577 | //std::cout << "P_mu: " << P_mu << std::endl; | |
578 | //std::cout << "q_mu: " << q_mu << std::endl; | |
579 | //std::cout << "q2: " << q2 << std::endl; | |
580 | //std::cout << "==================================================" << std::endl; | |
581 | ||
582 | return; | |
583 | } | |
584 | ||
585 | EvtTensor4C EvtLb2Lll::EvtLeptonTG5Current(const EvtDiracSpinor &d,const EvtDiracSpinor &dp){ | |
586 | // <u|sigma^munu*gamma^5|v> | |
587 | ||
588 | EvtTensor4C temp; | |
589 | temp.zero(); | |
590 | EvtComplex i2(0,0.5); | |
591 | ||
592 | static EvtGammaMatrix mat01=EvtGammaMatrix::g0()* | |
593 | (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()- | |
594 | EvtGammaMatrix::g1()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5(); | |
595 | static EvtGammaMatrix mat02=EvtGammaMatrix::g0()* | |
596 | (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()- | |
597 | EvtGammaMatrix::g2()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5(); | |
598 | static EvtGammaMatrix mat03=EvtGammaMatrix::g0()* | |
599 | (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()- | |
600 | EvtGammaMatrix::g3()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5(); | |
601 | static EvtGammaMatrix mat12=EvtGammaMatrix::g0()* | |
602 | (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()- | |
603 | EvtGammaMatrix::g2()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5(); | |
604 | static EvtGammaMatrix mat13=EvtGammaMatrix::g0()* | |
605 | (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()- | |
606 | EvtGammaMatrix::g3()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5(); | |
607 | static EvtGammaMatrix mat23=EvtGammaMatrix::g0()* | |
608 | (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()- | |
609 | EvtGammaMatrix::g3()*EvtGammaMatrix::g2())*EvtGammaMatrix::g5(); | |
610 | ||
611 | temp.set(0,1,i2*(d*(mat01*dp))); | |
612 | temp.set(1,0,-temp.get(0,1)); | |
613 | ||
614 | temp.set(0,2,i2*(d*(mat02*dp))); | |
615 | temp.set(2,0,-temp.get(0,2)); | |
616 | ||
617 | temp.set(0,3,i2*(d*(mat03*dp))); | |
618 | temp.set(3,0,-temp.get(0,3)); | |
619 | ||
620 | temp.set(1,2,i2*(d*(mat12*dp))); | |
621 | temp.set(2,1,-temp.get(1,2)); | |
622 | ||
623 | temp.set(1,3,i2*(d*(mat13*dp))); | |
624 | temp.set(3,1,-temp.get(1,3)); | |
625 | ||
626 | temp.set(2,3,i2*(d*(mat23*dp))); | |
627 | temp.set(3,2,-temp.get(2,3)); | |
628 | ||
629 | return temp; | |
630 | } |