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