evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtLb2Lll.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
21 #include "EvtGenModels/EvtLb2Lll.hh"
22 
23 #include "EvtGenBase/EvtComplex.hh"
24 #include "EvtGenBase/EvtDiracParticle.hh"
25 #include "EvtGenBase/EvtDiracSpinor.hh"
26 #include "EvtGenBase/EvtGammaMatrix.hh"
27 #include "EvtGenBase/EvtPDL.hh"
28 #include "EvtGenBase/EvtParticle.hh"
29 #include "EvtGenBase/EvtTensor4C.hh"
30 #include "EvtGenBase/EvtVector4C.hh"
31 #include "EvtGenBase/EvtVector4R.hh"
32 
33 #include <stdio.h>
34 #include <string.h>
35 
36 EvtDecayBase* EvtLb2Lll::clone()
37 {
38  return new EvtLb2Lll;
39 }
40 
41 std::string EvtLb2Lll::getName()
42 {
43  return "Lb2Lll";
44 }
45 
46 void EvtLb2Lll::init()
47 {
48  if ( getNArg() > 8 ) { // Decay parameters
49  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
50  << " ERROR: EvtLb2Lll generator expected max. 8 arguments but found: "
51  << getNArg() << std::endl;
52  EvtGenReport( EVTGEN_INFO, "EvtGen" )
53  << " 1. Lambda_b0 polarization - zero is default" << std::endl;
54  EvtGenReport( EVTGEN_INFO, "EvtGen" )
55  << " 2. Model type - \"SM\" for Standard Model is default"
56  << std::endl;
57  EvtGenReport( EVTGEN_INFO, "EvtGen" )
58  << " 3. Form-Factors - \"HQET\" is used by default" << std::endl;
59  EvtGenReport( EVTGEN_INFO, "EvtGen" )
60  << " 4. How to set polarization - \"ModifiedSpinors\" is default"
61  << std::endl;
62  EvtGenReport( EVTGEN_INFO, "EvtGen" )
63  << " 5. Include long distance (LD) effects - \"SD\" (no) is default"
64  << std::endl;
65  EvtGenReport( EVTGEN_INFO, "EvtGen" )
66  << " 6. NonFactorizable contribution (omega) to b->sg decay at q2=0 "
67  << std::endl;
68  EvtGenReport( EVTGEN_INFO, "EvtGen" )
69  << " 7. Note on every x-th decay" << std::endl;
70  EvtGenReport( EVTGEN_INFO, "EvtGen" )
71  << " 8. Maximum probability - automatic by default" << std::endl;
72  ::abort();
73  }
74 
75  if ( getNDaug() != 3 ) { // Check that there are 3 daughters only
76  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
77  << " ERROR: EvtLb2Lll generator expected 3 daughters but found: "
78  << getNDaug() << std::endl;
79  ::abort();
80  }
81 
82  EvtId LbID = EvtPDL::getId( std::string( "Lambda_b0" ) );
83  EvtId aLbID = EvtPDL::getId( std::string( "anti-Lambda_b0" ) );
84  EvtId eID = EvtPDL::getId( std::string( "e-" ) );
85  EvtId aeID = EvtPDL::getId( std::string( "e+" ) );
86  EvtId muID = EvtPDL::getId( std::string( "mu-" ) );
87  EvtId amuID = EvtPDL::getId( std::string( "mu+" ) );
88  EvtId tauID = EvtPDL::getId( std::string( "tau-" ) );
89  EvtId atauID = EvtPDL::getId( std::string( "tau+" ) );
90 
91  // TODO: better check based on spin and falvour is needed to allow usage of aliases !
92  if ( getParentId() == LbID ) { // Check daughters of Lambda_b0
93  EvtGenReport( EVTGEN_INFO, "EvtGen" )
94  << " EvtLb2Lll generator found Lambda_b0" << std::endl;
95  //if(EvtPDL::name(getDaug(0))!="Lambda0"){
96  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl;
97  // ::abort();
98  //}
99  if ( getDaug( 1 ) == eID && getDaug( 2 ) == aeID ) {
100  m_decayName = "Lambda_b0 -> Lambda0 e- e+";
101  EvtGenReport( EVTGEN_INFO, "EvtGen" )
102  << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 e- e+"
103  << std::endl;
104  } else if ( getDaug( 1 ) == muID && getDaug( 2 ) == amuID ) {
105  m_decayName = "Lambda_b0 -> Lambda0 mu- mu+";
106  EvtGenReport( EVTGEN_INFO, "EvtGen" )
107  << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 mu- mu+"
108  << std::endl;
109  } else if ( getDaug( 1 ) == tauID && getDaug( 2 ) == atauID ) {
110  m_decayName = "Lambda_b0 -> Lambda0 tau- tau+";
111  EvtGenReport( EVTGEN_INFO, "EvtGen" )
112  << " EvtLb2Lll generator found decay: Lambda_b0 -> Lambda0 tau- tau+"
113  << std::endl;
114  } else {
115  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
116  << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found: "
117  << EvtPDL::name( getDaug( 1 ) ) << " "
118  << EvtPDL::name( getDaug( 2 ) ) << std::endl;
119  ::abort();
120  }
121  //TODO: The model is known not to work correctly for anti-Lambda_b0 (A_FB does not change its sign)
122  } else if ( getParentId() == aLbID ) { // Check daughters of anti-Lambda_b0
123  EvtGenReport( EVTGEN_INFO, "EvtGen" )
124  << " EvtLb2Lll generator found anti-Lambda_b0" << std::endl;
125  //if(EvtPDL::name(getDaug(0))!="anti-Lambda0"){
126  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << " ERROR: EvtLb2Lll generator expected anti-Lambda0 daughter but found: " << EvtPDL::name(getDaug(0)) << std::endl;
127  // ::abort();
128  //}
129  if ( getDaug( 1 ) == aeID && getDaug( 2 ) == eID ) {
130  m_decayName = "anti-Lambda_b0 -> anti-Lambda0 e+ e-";
131  EvtGenReport( EVTGEN_INFO, "EvtGen" )
132  << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 e+ e-"
133  << std::endl;
134  } else if ( getDaug( 1 ) == amuID && getDaug( 2 ) == muID ) {
135  m_decayName = "anti-Lambda_b0 -> anti-Lambda0 mu+ mu-";
136  EvtGenReport( EVTGEN_INFO, "EvtGen" )
137  << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 mu+ mu-"
138  << std::endl;
139  } else if ( getDaug( 1 ) == atauID && getDaug( 2 ) == tauID ) {
140  m_decayName = "anti-Lambda_b0 -> anti-Lambda0 tau+ tau-";
141  EvtGenReport( EVTGEN_INFO, "EvtGen" )
142  << " EvtLb2Lll generator found decay: anti-Lambda_b0 -> anti-Lambda0 tau+ tau-"
143  << std::endl;
144  } else {
145  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
146  << " ERROR: EvtLb2Lll generator expected lepton pair daughters but found: "
147  << EvtPDL::name( getDaug( 1 ) ) << " "
148  << EvtPDL::name( getDaug( 2 ) ) << std::endl;
149  ::abort();
150  }
151  } else { // This model is not intended for decay of anything else than (anti-)Lambda_b0
152  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
153  << " ERROR: EvtLb2Lll generator expected (anti-)Lambda_b0 parent but found: "
154  << EvtPDL::name( getParentId() ) << std::endl;
155  ::abort();
156  }
157 
158  // Read and check all parameters
159  if ( getNArg() > 0 ) {
160  if ( getArg( 0 ) > 1. || getArg( 0 ) < -1. ) {
161  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
162  << " ERROR: EvtLb2Lll expects polarization to be in interval <-1,1>, not "
163  << getArg( 0 ) << std::endl;
164  ::abort();
165  }
166  m_polarizationLambdab0 = getArg( 0 );
167  } else {
168  m_polarizationLambdab0 = 0;
169  }
170  EvtGenReport( EVTGEN_INFO, "EvtGen" )
171  << " EvtLb2Lll set Lambda_b0 polarization to " << m_polarizationLambdab0
172  << std::endl;
173 
174  if ( getNArg() > 1 ) {
175  if ( getArgStr( 1 ).substr( 1, getArgStr( 1 ).size() - 2 ) != "SM" &&
176  getArgStr( 1 ).substr( 1, getArgStr( 1 ).size() - 2 ) != "-C7_SM" &&
177  getArgStr( 1 ).substr( 1, getArgStr( 1 ).size() - 2 ) !=
178  "SUSY-ChenGeng" ) {
179  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
180  << " ERROR: EvtLb2Lll doesn't know this physics model: "
181  << getArgStr( 1 ) << std::endl;
182  ::abort();
183  }
184  m_HEPmodel = getArgStr( 1 ).substr( 1, getArgStr( 1 ).size() - 2 );
185  } else {
186  m_HEPmodel = "SM";
187  }
188  EvtGenReport( EVTGEN_INFO, "EvtGen" )
189  << " EvtLb2Lll will use this physics model: " << m_HEPmodel << std::endl;
190 
191  if ( getNArg() > 2 ) {
192  if ( getArgStr( 2 ).substr( 1, getArgStr( 2 ).size() - 2 ) != "HQET" &&
193  getArgStr( 2 ).substr( 1, getArgStr( 2 ).size() - 2 ) != "HQET-noF2" &&
194  getArgStr( 2 ).substr( 1, 11 ) != "HQET-delta=" ) {
195  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
196  << " ERROR: EvtLb2Lll doesn't know this Form-Factors model: "
197  << getArgStr( 2 ) << std::endl;
198  ::abort();
199  }
200  m_FFtype = getArgStr( 2 ).substr( 1, getArgStr( 2 ).size() - 2 );
201  } else {
202  m_FFtype = "HQET";
203  }
204  EvtGenReport( EVTGEN_INFO, "EvtGen" )
205  << " EvtLb2Lll will use this Form-Factors model: " << m_FFtype
206  << std::endl;
207 
208  if ( getNArg() > 3 ) {
209  if ( getArgStr( 3 ).substr( 1, getArgStr( 3 ).size() - 2 ) !=
210  "Unpolarized" ) {
211  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
212  << " ERROR: EvtLb2Lll doesn't know kind of introducing polarization: "
213  << getArgStr( 3 ) << std::endl;
214  ::abort();
215  }
216  m_polarizationIntroduction =
217  getArgStr( 3 ).substr( 1, getArgStr( 3 ).size() - 2 );
218  } else {
219  m_polarizationIntroduction = "Unpolarized";
220  }
221  EvtGenReport( EVTGEN_INFO, "EvtGen" )
222  << " EvtLb2Lll will use this kind of introducing polarization: "
223  << m_polarizationIntroduction << std::endl;
224 
225  if ( getNArg() > 4 ) {
226  if ( getArgStr( 4 ).substr( 1, getArgStr( 4 ).size() - 2 ) != "SD" &&
227  getArgStr( 4 ).substr( 1, getArgStr( 4 ).size() - 2 ) != "LD" ) {
228  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
229  << " ERROR: EvtLb2Lll didn't find SD or LD parameter: "
230  << getArgStr( 4 ) << std::endl;
231  ::abort();
232  }
233  m_effectContribution = getArgStr( 5 ).substr( 1, getArgStr( 4 ).size() -
234  2 );
235  } else {
236  m_effectContribution = "SD";
237  }
238  EvtGenReport( EVTGEN_INFO, "EvtGen" )
239  << " EvtLb2Lll will include contribution from these effects: "
240  << m_effectContribution << std::endl;
241 
242  if ( getNArg() > 5 ) {
243  if ( fabs( getArg( 5 ) ) > 0.15 ) {
244  EvtGenReport( EVTGEN_WARNING, "EvtGen" )
245  << " WARNING: EvtLb2Lll found very high contribution to b->sg decay at q2=0: "
246  << getArg( 5 ) << std::endl;
247  }
248  m_omega = getArg( 5 );
249  } else {
250  m_omega = 0;
251  }
252  EvtGenReport( EVTGEN_INFO, "EvtGen" )
253  << " EvtLb2Lll will use this contribution to b->sg decay at q2=0: "
254  << m_omega << std::endl;
255 
256  if ( getNArg() > 6 )
257  m_noTries = (long)( getArg( 6 ) );
258  else
259  m_noTries = 0;
260 
261  if ( getNArg() > 7 ) {
262  if ( getArg( 7 ) < 0. ) {
263  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
264  << " ERROR: EvtLb2Lll expects positive maximum probability not : "
265  << getArg( 7 ) << std::endl;
266  ::abort();
267  }
268  m_maxProbability = getArg( 7 );
269  } else {
270  m_maxProbability = 0.;
271  }
272  EvtGenReport( EVTGEN_INFO, "EvtGen" )
273  << " EvtLb2Lll maximum probability was set to " << m_maxProbability
274  << std::endl;
275  m_poleSize = 0;
276 
277  // Initialize Wilson coefficients by Buras and Munz
278  // TODO: should have common W.C. source for all decays in EvtGen
279  m_WC.CalculateAllCoefficients();
280 }
281 
282 void EvtLb2Lll::initProbMax()
283 {
284  EvtGenReport( EVTGEN_INFO, "EvtGen" )
285  << " EvtLb2Lll is finding maximum probability ... " << std::endl;
286 
287  if ( m_maxProbability == 0 ) {
288  EvtDiracParticle* parent = new EvtDiracParticle;
289  parent->noLifeTime();
290  parent->init( getParentId(),
291  EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) );
292  parent->setDiagonalSpinDensity();
293 
294  EvtAmp amp;
295  EvtId daughters[3] = {getDaug( 0 ), getDaug( 1 ), getDaug( 2 )};
296  amp.init( getParentId(), 3, daughters );
297  parent->makeDaughters( 3, daughters );
298  EvtParticle* lambda = parent->getDaug( 0 );
299  EvtParticle* lep1 = parent->getDaug( 1 );
300  EvtParticle* lep2 = parent->getDaug( 2 );
301  lambda->noLifeTime();
302  lep1->noLifeTime();
303  lep2->noLifeTime();
304 
305  EvtSpinDensity rho;
306  rho.setDiag( parent->getSpinStates() );
307 
308  double M0 = EvtPDL::getMass( getParentId() );
309  double mL = EvtPDL::getMass( getDaug( 0 ) );
310  double m1 = EvtPDL::getMass( getDaug( 1 ) );
311  double m2 = EvtPDL::getMass( getDaug( 2 ) );
312 
313  double q2, pstar, elambda, theta;
314  double q2min = ( m1 + m2 ) * ( m1 + m2 );
315  double q2max = ( M0 - mL ) * ( M0 - mL );
316 
317  EvtVector4R p4lambda, p4lep1, p4lep2, boost;
318 
319  EvtGenReport( EVTGEN_INFO, "EvtGen" )
320  << " EvtLb2Lll is probing whole phase space ..." << std::endl;
321 
322  int i, j;
323  double prob = 0;
324  for ( i = 0; i <= 100; i++ ) {
325  q2 = q2min + i * ( q2max - q2min ) / 100.;
326  elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0;
327  if ( i == 0 )
328  pstar = 0;
329  else
330  pstar = sqrt( q2 - ( m1 + m2 ) * ( m1 + m2 ) ) *
331  sqrt( q2 - ( m1 - m2 ) * ( m1 - m2 ) ) / 2 / sqrt( q2 );
332  boost.set( M0 - elambda, 0, 0, +sqrt( elambda * elambda - mL * mL ) );
333  if ( i != 100 ) {
334  p4lambda.set( elambda, 0, 0,
335  -sqrt( elambda * elambda - mL * mL ) );
336  } else {
337  p4lambda.set( mL, 0, 0, 0 );
338  }
339  for ( j = 0; j <= 45; j++ ) {
340  theta = j * EvtConst::pi / 45;
341  p4lep1.set( sqrt( pstar * pstar + m1 * m1 ), 0,
342  +pstar * sin( theta ), +pstar * cos( theta ) );
343  p4lep2.set( sqrt( pstar * pstar + m2 * m2 ), 0,
344  -pstar * sin( theta ), -pstar * cos( theta ) );
345  if ( i != 100 ) // At maximal q2 we are already in correct frame as Lambda and W/Zvirtual are at rest
346  {
347  p4lep1 = boostTo( p4lep1, boost );
348  p4lep2 = boostTo( p4lep2, boost );
349  }
350  calcAmp( &amp, parent );
351  prob = rho.normalizedProb( amp.getSpinDensity() );
352  //std::cout << "q2: " << q2 << " \t theta: " << theta << " \t prob: " << prob << std::endl;
353  //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl;
354  if ( prob > m_maxProbability ) {
355  EvtGenReport( EVTGEN_INFO, "EvtGen" )
356  << " - probability " << prob << " found at q2 = " << q2
357  << " (" << 100 * ( q2 - q2min ) / ( q2max - q2min )
358  << " %) and theta = " << theta * 180 / EvtConst::pi
359  << std::endl;
360  m_maxProbability = prob;
361  }
362  }
363  //::abort();
364  }
365 
366  //m_poleSize = 0.04*q2min;
367  m_maxProbability *= 1.2;
368  delete parent;
369  }
370 
371  setProbMax( m_maxProbability );
372  EvtGenReport( EVTGEN_INFO, "EvtGen" )
373  << " EvtLb2Lll set up maximum probability to " << m_maxProbability
374  << std::endl;
375 }
376 
377 void EvtLb2Lll::decay( EvtParticle* parent )
378 {
379  //setWeight(parent->initializePhaseSpace(getNDaug(),getDaugs(),m_poleSize,1,2));
380  parent->initializePhaseSpace( getNDaug(), getDaugs() );
381  calcAmp( &_amp2, parent );
382 }
383 
384 void EvtLb2Lll::calcAmp( EvtAmp* amp, EvtParticle* parent )
385 {
386  static long noTries = 0;
387  static double delta = 0;
388 
389  EvtComplex Matrix[2][2][2][2];
390 
391  EvtComplex i1( 0, 1 );
392 
393  int i, j, spins[4];
394  char ch;
395 
396  double r, M_L, M_Lb, M_s, M_c, M_b, q2, alpha, M_W, M_t;
397  double M_psi[2] = {0, 0}, Gamma_psi[2] = {0, 0}, k_psi[2] = {0, 0};
398  double F0_1, F0_2, a_F1, a_F2, b_F1, b_F2, F1, F2;
399  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,
400  f_TV, f_TS, g_TV( 0.0 ), g_TS, f_T, g_T;
401  EvtComplex A1, A2, A3, B1, B2, B3, D1, D2, D3, E1, E2, E3, N1, N2, H1, H2;
402  EvtComplex C_SL, C_BR, C_LLtot, C_LRtot, C_LL, C_LR, C_RL, C_RR, C_LRLR,
403  C_RLLR, C_LRRL, C_RLRL, C_T, C_TE;
404  EvtComplex Yld, C_7eff, C_9eff;
405  EvtComplex V_ts, V_tb;
406 
407  EvtVector4C lbar_Gmu_l[2][2], lbar_GmuG5_l[2][2], hbar_GmuPlusG5_h[2][2],
408  hbar_GmuMinusG5_h[2][2], hbar_Gmu_h[2][2];
409  EvtComplex lbar_l[2][2], lbar_G5_l[2][2], hbar_1PlusG5_h[2][2],
410  hbar_1MinusG5_h[2][2], hbar_G5_h[2][2], hbar_h[2][2];
411  EvtTensor4C lbar_Smunu_l[2][2], lbar_ESmunu_l[2][2],
412  hbar_SmunuPlusG5_h[2][2], hbar_SmunuMinusG5_h[2][2], hbar_Smunu_h[2][2];
413  EvtVector4R q_mu, P_mu;
414 
415  EvtDiracSpinor parent__spParent[2];
416 
417  M_Lb = parent->mass();
418  M_L = parent->getDaug( 0 )->mass();
419  M_s = 0.13;
420  M_c = 1.35;
421  M_b = 4.8;
422  alpha = 1. / 137.036;
423  M_W = 80.425;
424  M_t = 174.3;
425  M_psi[0] = 3.096916;
426  M_psi[1] = 3.686093;
427  if ( m_decayName == "Lambda_b0 -> Lambda0 e- e+" ||
428  m_decayName == "anti-Lambda_b0 -> anti-Lambda0 e+ e-" ) {
429  Gamma_psi[0] = 5.40;
430  Gamma_psi[1] = 2.12;
431  }
432  if ( m_decayName == "Lambda_b0 -> Lambda0 mu- mu+" ||
433  m_decayName == "anti-Lambda_b0 -> anti-Lambda0 mu+ mu-" ) {
434  Gamma_psi[0] = 5.35;
435  Gamma_psi[1] = 2.05;
436  }
437  if ( m_decayName == "Lambda_b0 -> Lambda0 tau- tau+" ||
438  m_decayName == "anti-Lambda_b0 -> anti-Lambda0 tau+ tau-" ) {
439  Gamma_psi[0] = 0.00;
440  Gamma_psi[1] = 0.79;
441  }
442  if ( m_effectContribution == "LD" ) {
443  k_psi[0] = 1.65;
444  k_psi[1] = 1.65;
445  }
446  //G_F = 1.16637e-5;
447  //V_tb = sqrt(1-pow(0.0413,2))*sqrt(1-pow(0.0037,2));
448  //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));
449 
450  P_mu = parent->getP4Restframe() + parent->getDaug( 0 )->getP4();
451  q_mu = parent->getP4Restframe() - parent->getDaug( 0 )->getP4();
452  q2 = q_mu.mass2();
453 
454  if ( m_noTries > 0 )
455  if ( !( ( ++noTries ) % m_noTries ) )
456  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
457  << " EvtLb2Lll already finished " << noTries
458  << " matrix element calculations" << std::endl;
459 
460  if ( m_FFtype == "HQET" ) {
461  r = M_L * M_L / M_Lb / M_Lb;
462  F0_1 = +0.462;
463  F0_2 = -0.077;
464  a_F1 = -0.0182;
465  a_F2 = -0.0685;
466  b_F1 = -0.000176;
467  b_F2 = +0.001460;
468  F1 = F0_1 / ( 1.0 - ( q2 / M_Lb / M_Lb ) *
469  ( a_F1 - b_F1 * ( q2 / M_Lb / M_Lb ) ) );
470  F2 = F0_2 / ( 1.0 - ( q2 / M_Lb / M_Lb ) *
471  ( a_F2 - b_F2 * ( q2 / M_Lb / M_Lb ) ) );
472  g_1 = f_1 = f_2T = g_2T = F1 + sqrt( r ) * F2;
473  //std::cout << " F1: " << F1 << " F2: " << F2 << " r: " << r << " M_L: " << M_L << " M_Lb: " << M_Lb << std::endl;
474  //std::cout << " sqrt(q2): " << sqrt(q2) << " q2: " << q2 << " M_Lb^2" << M_Lb*M_Lb << std::endl;
475  g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2 / M_Lb;
476  g_TS = f_TS = 0;
477  g_1T = f_1T = F2 / M_Lb * q2;
478  g_3T = +F2 / M_Lb * ( M_Lb + M_L );
479  f_3T = -F2 / M_Lb * ( M_Lb - M_L );
480  f_T = f_2T - f_TS * q2;
481  g_T = g_2T - g_TS * q2;
482  } else if ( strstr( m_FFtype.c_str(), "HQET-delta=" ) == m_FFtype.c_str() ) {
483  //EvtGenReport(EVTGEN_WARNING,"EvtGen") << " WARNING: HQET-delta FF model should be checked for correctness" << std::endl;
484  if ( delta == 0 )
485  sscanf( m_FFtype.c_str(), "%c%c%c%c%c%c%c%c%c%c%c%lf", &ch, &ch,
486  &ch, &ch, &ch, &ch, &ch, &ch, &ch, &ch, &ch, &delta );
487  r = M_L * M_L / M_Lb / M_Lb;
488  F0_1 = +0.462;
489  F0_2 = -0.077;
490  a_F1 = -0.0182;
491  a_F2 = -0.0685;
492  b_F1 = -0.000176;
493  b_F2 = +0.001460;
494  F1 = F0_1 / ( 1.0 - ( q2 / M_Lb / M_Lb ) *
495  ( a_F1 - b_F1 * ( q2 / M_Lb / M_Lb ) ) );
496  F2 = F0_2 / ( 1.0 - ( q2 / M_Lb / M_Lb ) *
497  ( a_F2 - b_F2 * ( q2 / M_Lb / M_Lb ) ) );
498  g_1 = f_1 = f_2T = g_2T = F1 + sqrt( r ) * F2;
499  g_1 += delta * g_1;
500  f_1 -= delta * f_1;
501  g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = F2 / M_Lb;
502  g_TS = f_TS = 0;
503  g_1T = f_1T = F2 / M_Lb * q2;
504  g_3T = +F2 / M_Lb * ( M_Lb + M_L );
505  f_3T = -F2 / M_Lb * ( M_Lb - M_L );
506  f_T = f_2T - f_TS * q2;
507  g_T = g_2T - g_TS * q2;
508  } else if ( m_FFtype == "HQET-noF2" ) {
509  //EvtGenReport(EVTGEN_WARNING,"EvtGen") << " WARNING: HQET-noF2 FF model should be checked for correctness" << std::endl;
510  r = M_L * M_L / M_Lb / M_Lb;
511  F0_1 = +0.462;
512  a_F1 = -0.0182;
513  b_F1 = -0.000176;
514  F1 = F0_1 / ( 1.0 - ( q2 / M_Lb / M_Lb ) *
515  ( a_F1 - b_F1 * ( q2 / M_Lb / M_Lb ) ) );
516  g_1 = f_1 = f_2T = g_2T = F1;
517  g_2 = f_2 = g_3 = f_3 = g_TV = f_TV = 0;
518  g_TS = f_TS = 0;
519  g_1T = f_1T = 0;
520  g_3T = 0;
521  f_3T = 0;
522  f_T = f_2T - f_TS * q2;
523  g_T = g_2T - g_TS * q2;
524  } else { // general relations for Form-Factors
525  f_1 = f_2 = f_3 = g_1 = g_2 = g_3 = f_3T = g_3T = f_TS = g_TS = f_T =
526  g_T = f_TV = 0;
527  f_2T = f_T + f_TS * q2;
528  f_1T = ( f_TV + f_TS * ( M_L + M_Lb ) ) * q2;
529  f_1T = -q2 / ( M_Lb - M_L ) * f_3T;
530  g_2T = g_T + g_TS * q2;
531  g_1T = ( g_TV - g_TS * ( M_L - M_Lb ) ) * q2;
532  g_1T = +q2 / ( M_Lb + M_L ) * g_3T;
533  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
534  << " ERROR: EvtLb2Lll - unknown Form-Factors model: " << m_FFtype
535  << " - this should never happen !" << std::endl;
536  ::abort();
537  }
538 
539  if ( m_HEPmodel == "SM" ) {
540  C_LL = C_LR = C_RL = C_RR = C_LRLR = C_RLLR = C_LRRL = C_RLRL = C_T =
541  C_TE = EvtComplex( 0, 0 );
542  Yld = m_WC.Yld( q2, k_psi, Gamma_psi, M_psi, 2, m_WC.GetC1(),
543  m_WC.GetC2(), m_WC.GetC3(), m_WC.GetC4(), m_WC.GetC5(),
544  m_WC.GetC6(), 1. / alpha );
545  C_7eff = m_WC.GetC7eff0() +
546  m_WC.C7b2sg( m_WC.GetStrongCouplingConst(), m_WC.GetEta(),
547  m_WC.GetC2(), M_t, M_W ) +
548  m_omega *
549  ( m_WC.hzs( M_c / M_b, q2 / M_b / M_b, M_b, M_b ) + Yld );
550  C_9eff = Yld +
551  m_WC.C9efftilda( M_c / M_b, q2 / M_b / M_b,
552  m_WC.GetStrongCouplingConst(), m_WC.GetC1(),
553  m_WC.GetC2(), m_WC.GetC3(), m_WC.GetC4(),
554  m_WC.GetC5(), m_WC.GetC6(), m_WC.GetC9tilda(),
555  m_WC.GetRenormSchemePar() );
556  C_SL = -2 * M_s * C_7eff;
557  C_BR = -2 * M_b * C_7eff;
558  C_LLtot = C_9eff - m_WC.GetC10tilda() + C_LL;
559  C_LRtot = C_9eff + m_WC.GetC10tilda() + C_LR;
560  //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;
561  } else if ( m_HEPmodel == "-C7_SM" ) {
562  C_LL = C_LR = C_RL = C_RR = C_LRLR = C_RLLR = C_LRRL = C_RLRL = C_T =
563  C_TE = EvtComplex( 0, 0 );
564  Yld = m_WC.Yld( q2, k_psi, Gamma_psi, M_psi, 2, m_WC.GetC1(),
565  m_WC.GetC2(), m_WC.GetC3(), m_WC.GetC4(), m_WC.GetC5(),
566  m_WC.GetC6(), 1. / alpha );
567  C_7eff = m_WC.GetC7eff0() +
568  m_WC.C7b2sg( m_WC.GetStrongCouplingConst(), m_WC.GetEta(),
569  m_WC.GetC2(), M_t, M_W ) +
570  m_omega *
571  ( m_WC.hzs( M_c / M_b, q2 / M_b / M_b, M_b, M_b ) + Yld );
572  C_9eff = Yld +
573  m_WC.C9efftilda( M_c / M_b, q2 / M_b / M_b,
574  m_WC.GetStrongCouplingConst(), m_WC.GetC1(),
575  m_WC.GetC2(), m_WC.GetC3(), m_WC.GetC4(),
576  m_WC.GetC5(), m_WC.GetC6(), m_WC.GetC9tilda(),
577  m_WC.GetRenormSchemePar() );
578  C_SL = +2 * M_s * C_7eff;
579  C_BR = +2 * M_b * C_7eff;
580  C_LLtot = C_9eff - m_WC.GetC10tilda() + C_LL;
581  C_LRtot = C_9eff + m_WC.GetC10tilda() + C_LR;
582  //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;
583  } else if ( m_HEPmodel == "SUSY-ChenGeng" ) {
584  //EvtGenReport(EVTGEN_WARNING,"EvtGen") << " WARNING: SUSY-ChenGeng model should be checked for correctness" << std::endl;
585  C_LL = C_LR = C_RL = C_RR = C_LRLR = C_RLLR = C_LRRL = C_RLRL = C_T =
586  C_TE = EvtComplex( 0, 0 );
587  EvtComplex d_u23LL = 0.1;
588  EvtComplex d_u33RL = 0.65;
589  EvtComplex d_d23LR = 0.03 * exp( i1 * EvtConst::pi * 2 / 5 );
590  EvtComplex d_u23LR = -0.8 * exp( i1 * EvtConst::pi / 4 );
591  EvtComplex C_7susy = -1.75 * d_u23LL - 0.25 * d_u23LR - 10.3 * d_d23LR;
592  EvtComplex C_9susy = 0.82 * d_u23LR;
593  EvtComplex C_10susy = -9.37 * d_u23LR + 1.4 * d_u23LR * d_u33RL +
594  2.7 * d_u23LL;
595  Yld = m_WC.Yld( q2, k_psi, Gamma_psi, M_psi, 2, m_WC.GetC1(),
596  m_WC.GetC2(), m_WC.GetC3(), m_WC.GetC4(), m_WC.GetC5(),
597  m_WC.GetC6(), 1. / alpha );
598  C_7eff = m_WC.GetC7eff0() + C_7susy * pow( m_WC.GetEta(), 16. / 23. ) +
599  m_WC.C7b2sg( m_WC.GetStrongCouplingConst(), m_WC.GetEta(),
600  m_WC.GetC2(), M_t, M_W ) +
601  m_omega *
602  ( m_WC.hzs( M_c / M_b, q2 / M_b / M_b, M_b, M_b ) + Yld );
603  C_9eff = Yld + m_WC.C9efftilda( M_c / M_b, q2 / M_b / M_b,
604  m_WC.GetStrongCouplingConst(),
605  m_WC.GetC1(), m_WC.GetC2(), m_WC.GetC3(),
606  m_WC.GetC4(), m_WC.GetC5(), m_WC.GetC6(),
607  m_WC.GetC9tilda() + C_9susy,
608  m_WC.GetRenormSchemePar() );
609  C_SL = -2 * M_s * C_7eff;
610  C_BR = -2 * M_b * C_7eff;
611  C_LLtot = C_9eff - m_WC.GetC10tilda() - C_10susy + C_LL;
612  C_LRtot = C_9eff + m_WC.GetC10tilda() + C_10susy + C_LR;
613  //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;
614  } else {
615  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
616  << " ERROR: EvtLb2Lll - unknown physics model: " << m_HEPmodel
617  << " - this should never happen !" << std::endl;
618  ::abort();
619  }
620 
621  A1 = ( f_1T - g_1T ) * C_SL / q2 + ( f_1T + g_1T ) * C_BR / q2 +
622  0.5 * ( f_1 - g_1 ) * ( C_LLtot + C_LRtot ) +
623  0.5 * ( f_1 + g_1 ) * ( C_RL + C_RR );
624  //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;
625  A2 = ( f_2T - g_2T ) * C_SL / q2 + ( f_2T + g_2T ) * C_BR / q2 +
626  0.5 * ( f_2 - g_2 ) * ( C_LLtot + C_LRtot ) +
627  0.5 * ( f_2 + g_2 ) * ( C_RL + C_RR );
628  A3 = ( f_3T - g_3T ) * C_SL / q2 + ( f_3T + g_3T ) * C_BR / q2 +
629  0.5 * ( f_3 - g_3 ) * ( C_LLtot + C_LRtot ) +
630  0.5 * ( f_3 + g_3 ) * ( C_RL + C_RR );
631 
632  B1 = ( f_1T + g_1T ) * C_SL / q2 + ( f_1T - g_1T ) * C_BR / q2 +
633  0.5 * ( f_1 + g_1 ) * ( C_LLtot + C_LRtot ) +
634  0.5 * ( f_1 - g_1 ) * ( C_RL + C_RR );
635  B2 = ( f_2T + g_2T ) * C_SL / q2 + ( f_2T - g_2T ) * C_BR / q2 +
636  0.5 * ( f_2 + g_2 ) * ( C_LLtot + C_LRtot ) +
637  0.5 * ( f_2 - g_2 ) * ( C_RL + C_RR );
638  B3 = ( f_3T + g_3T ) * C_SL / q2 + ( f_3T - g_3T ) * C_BR / q2 +
639  0.5 * ( f_3 + g_3 ) * ( C_LLtot + C_LRtot ) +
640  0.5 * ( f_3 - g_3 ) * ( C_RL + C_RR );
641 
642  D1 = 0.5 * ( C_RR - C_RL ) * ( f_1 + g_1 ) +
643  0.5 * ( C_LRtot - C_LLtot ) * ( f_1 - g_1 );
644  D2 = 0.5 * ( C_RR - C_RL ) * ( f_2 + g_2 ) +
645  0.5 * ( C_LRtot - C_LLtot ) * ( f_2 - g_2 );
646  D3 = 0.5 * ( C_RR - C_RL ) * ( f_3 + g_3 ) +
647  0.5 * ( C_LRtot - C_LLtot ) * ( f_3 - g_3 );
648 
649  E1 = 0.5 * ( C_RR - C_RL ) * ( f_1 - g_1 ) +
650  0.5 * ( C_LRtot - C_LLtot ) * ( f_1 + g_1 );
651  E2 = 0.5 * ( C_RR - C_RL ) * ( f_2 - g_2 ) +
652  0.5 * ( C_LRtot - C_LLtot ) * ( f_2 + g_2 );
653  E3 = 0.5 * ( C_RR - C_RL ) * ( f_3 - g_3 ) +
654  0.5 * ( C_LRtot - C_LLtot ) * ( f_3 + g_3 );
655 
656  N1 = ( f_1 * ( M_Lb - M_L ) + f_3 * q2 ) / M_b *
657  ( C_LRLR + C_RLLR + C_LRRL + C_RLRL ); // Should be mLb - mL
658  N2 = ( f_1 * ( M_Lb - M_L ) + f_3 * q2 ) / M_b *
659  ( C_LRLR + C_RLLR - C_LRRL - C_RLRL );
660 
661  H1 = ( g_1 * ( M_Lb + M_L ) - g_3 * q2 ) / M_b *
662  ( C_LRLR - C_RLLR + C_LRRL - C_RLRL );
663  H2 = ( g_1 * ( M_Lb + M_L ) - g_3 * q2 ) / M_b *
664  ( C_LRLR - C_RLLR - C_LRRL + C_RLRL );
665 
666  for ( i = 0; i < 4; i++ ) {
667  lbar_Gmu_l[i / 2][i % 2] =
668  EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( i / 2 ),
669  parent->getDaug( 2 )->spParent( i % 2 ) );
670  lbar_GmuG5_l[i / 2][i % 2] =
671  EvtLeptonACurrent( parent->getDaug( 1 )->spParent( i / 2 ),
672  parent->getDaug( 2 )->spParent( i % 2 ) );
673  lbar_l[i / 2][i % 2] =
674  EvtLeptonSCurrent( parent->getDaug( 1 )->spParent( i / 2 ),
675  parent->getDaug( 2 )->spParent( i % 2 ) );
676  lbar_G5_l[i / 2][i % 2] =
677  EvtLeptonPCurrent( parent->getDaug( 1 )->spParent( i / 2 ),
678  parent->getDaug( 2 )->spParent( i % 2 ) );
679  lbar_Smunu_l[i / 2][i % 2] =
680  EvtLeptonTCurrent( parent->getDaug( 1 )->spParent( i / 2 ),
681  parent->getDaug( 2 )->spParent( i % 2 ) );
682  lbar_ESmunu_l[i / 2][i % 2] = dual(
683  EvtLeptonTCurrent( parent->getDaug( 1 )->spParent( i / 2 ),
684  parent->getDaug( 2 )->spParent( i % 2 ) ) );
685  }
686 
687  // TODO: polarization not yet introduced
688  if ( m_polarizationIntroduction == "SpinDensityMatrix" ) {
689  //parent->setSpinDensityForward();
690  parent__spParent[0] = parent->sp( 0 );
691  parent__spParent[1] = parent->sp( 1 );
692  } else if ( m_polarizationIntroduction == "ModifiedSpinors" ) {
693  parent__spParent[0] = parent->sp( 0 );
694  parent__spParent[1] = parent->sp( 1 );
695  } else if ( m_polarizationIntroduction == "Unpolarized" ) {
696  parent__spParent[0] = parent->sp( 0 );
697  parent__spParent[1] = parent->sp( 1 );
698  } else {
699  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
700  << " ERROR: EvtLb2Lll - unknown polarization: "
701  << m_polarizationIntroduction << " - this should never happen !"
702  << std::endl;
703  ::abort();
704  }
705 
706  for ( i = 0; i < 4; i++ ) {
707  hbar_GmuPlusG5_h[i / 2][i % 2] =
708  EvtLeptonVCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
709  parent__spParent[i % 2] ) +
710  EvtLeptonACurrent( parent->getDaug( 0 )->spParent( i / 2 ),
711  parent__spParent[i % 2] );
712  hbar_GmuMinusG5_h[i / 2][i % 2] = EvtLeptonVACurrent(
713  parent->getDaug( 0 )->spParent( i / 2 ), parent__spParent[i % 2] );
714  hbar_SmunuPlusG5_h[i / 2][i % 2] =
715  EvtLeptonTCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
716  parent__spParent[i % 2] ) +
717  EvtLeptonTG5Current( parent->getDaug( 0 )->spParent( i / 2 ),
718  parent__spParent[i % 2] );
719  hbar_SmunuMinusG5_h[i / 2][i % 2] =
720  EvtLeptonTCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
721  parent__spParent[i % 2] ) -
722  EvtLeptonTG5Current( parent->getDaug( 0 )->spParent( i / 2 ),
723  parent__spParent[i % 2] );
724  hbar_1PlusG5_h[i / 2][i % 2] =
725  EvtLeptonSCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
726  parent__spParent[i % 2] ) +
727  EvtLeptonPCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
728  parent__spParent[i % 2] );
729  hbar_1MinusG5_h[i / 2][i % 2] =
730  EvtLeptonSCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
731  parent__spParent[i % 2] ) -
732  EvtLeptonPCurrent( parent->getDaug( 0 )->spParent( i / 2 ),
733  parent__spParent[i % 2] );
734  hbar_G5_h[i / 2][i % 2] = EvtLeptonPCurrent(
735  parent->getDaug( 0 )->spParent( i / 2 ), parent__spParent[i % 2] );
736  hbar_h[i / 2][i % 2] = EvtLeptonSCurrent(
737  parent->getDaug( 0 )->spParent( i / 2 ), parent__spParent[i % 2] );
738  hbar_Smunu_h[i / 2][i % 2] = EvtLeptonTCurrent(
739  parent->getDaug( 0 )->spParent( i / 2 ), parent__spParent[i % 2] );
740  hbar_Gmu_h[i / 2][i % 2] = EvtLeptonVCurrent(
741  parent->getDaug( 0 )->spParent( i / 2 ), parent__spParent[i % 2] );
742  }
743 
744  for ( i = 0; i < 4; i++ )
745  for ( j = 0; j < 4; j++ ) {
746  //std::cout << "--------------------------------------------------" << std::endl;
747  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
748  Matrix[j / 2][j % 2][i / 2][i % 2] +=
749  lbar_Gmu_l[i / 2][i % 2] *
750  ( A1 * hbar_GmuPlusG5_h[j / 2][j % 2] +
751  B1 * hbar_GmuMinusG5_h[j / 2][j % 2] );
752  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
753  //std::cout << "A1: " << A1 << " B1: " << B1 << " lbar_Gmu_l: " << lbar_Gmu_l[i/2][i%2] <<
754  // " hbar_GmuPlusG5_h: " << hbar_GmuPlusG5_h[j/2][j%2] << " hbar_GmuMinusG5_h: " << hbar_GmuMinusG5_h[j/2][j%2] <<
755  // " sp1: " << parent->getDaug(1)->spParent(i/2) << " sp2: " << parent->getDaug(1)->spParent(i%2) << std::endl;
756  Matrix[j / 2][j % 2][i / 2][i % 2] +=
757  lbar_Gmu_l[i / 2][i % 2] *
758  ( i1 * A2 * ( hbar_SmunuPlusG5_h[j / 2][j % 2].cont2( q_mu ) ) +
759  B2 * ( hbar_SmunuMinusG5_h[j / 2][j % 2].cont2( q_mu ) ) );
760  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
761  Matrix[j / 2][j % 2][i / 2][i % 2] +=
762  lbar_Gmu_l[i / 2][i % 2] *
763  ( ( A3 * hbar_1PlusG5_h[j / 2][j % 2] +
764  B3 * hbar_1MinusG5_h[j / 2][j % 2] ) *
765  q_mu );
766  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
767  Matrix[j / 2][j % 2][i / 2][i % 2] +=
768  lbar_GmuG5_l[i / 2][i % 2] *
769  ( D1 * hbar_GmuPlusG5_h[j / 2][j % 2] +
770  E1 * hbar_GmuMinusG5_h[j / 2][j % 2] );
771  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
772  Matrix[j / 2][j % 2][i / 2][i % 2] +=
773  lbar_GmuG5_l[i / 2][i % 2] *
774  ( i1 * D2 * ( hbar_SmunuPlusG5_h[j / 2][j % 2].cont2( q_mu ) ) +
775  E2 * ( hbar_SmunuMinusG5_h[j / 2][j % 2].cont2( q_mu ) ) );
776  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
777  Matrix[j / 2][j % 2][i / 2][i % 2] +=
778  lbar_GmuG5_l[i / 2][i % 2] *
779  ( ( D3 * hbar_1PlusG5_h[j / 2][j % 2] +
780  E3 * hbar_1MinusG5_h[j / 2][j % 2] ) *
781  q_mu );
782  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
783  Matrix[j / 2][j % 2][i / 2][i % 2] +=
784  lbar_l[i / 2][i % 2] *
785  ( N1 * hbar_h[j / 2][j % 2] + H1 * hbar_G5_h[j / 2][j % 2] );
786  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
787  Matrix[j / 2][j % 2][i / 2][i % 2] +=
788  lbar_G5_l[i / 2][i % 2] *
789  ( N2 * hbar_h[j / 2][j % 2] + H2 * hbar_G5_h[j / 2][j % 2] );
790  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
791  Matrix[j / 2][j % 2][i / 2][i % 2] +=
792  cont( lbar_Smunu_l[i / 2][i % 2],
793  4 * C_T * f_T * hbar_Smunu_h[j / 2][j % 2] );
794  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
795  Matrix[j / 2][j % 2][i / 2][i % 2] += cont(
796  lbar_Smunu_l[i / 2][i % 2],
797  -4 * C_T * f_TV * i1 *
798  ( EvtGenFunctions::directProd( q_mu, hbar_Gmu_h[j / 2][j % 2] ) -
799  EvtGenFunctions::directProd( hbar_Gmu_h[j / 2][j % 2],
800  q_mu ) ) );
801  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
802  Matrix[j / 2][j % 2][i / 2][i % 2] +=
803  cont( lbar_Smunu_l[i / 2][i % 2],
804  -4 * C_T * f_TS * i1 *
805  ( EvtGenFunctions::directProd( P_mu, q_mu ) -
806  EvtGenFunctions::directProd( q_mu, P_mu ) ) *
807  hbar_h[j / 2][j % 2] );
808  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
809  Matrix[j / 2][j % 2][i / 2][i % 2] +=
810  cont( lbar_ESmunu_l[i / 2][i % 2],
811  4 * C_TE * f_T * i1 * hbar_Smunu_h[j / 2][j % 2] );
812  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
813  Matrix[j / 2][j % 2][i / 2][i % 2] += cont(
814  lbar_ESmunu_l[i / 2][i % 2],
815  4 * C_TE * f_TV *
816  ( EvtGenFunctions::directProd( q_mu, hbar_Gmu_h[j / 2][j % 2] ) -
817  EvtGenFunctions::directProd( hbar_Gmu_h[j / 2][j % 2],
818  q_mu ) ) );
819  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
820  Matrix[j / 2][j % 2][i / 2][i % 2] +=
821  cont( lbar_ESmunu_l[i / 2][i % 2],
822  4 * C_TE * f_TS *
823  ( EvtGenFunctions::directProd( P_mu, q_mu ) -
824  EvtGenFunctions::directProd( q_mu, P_mu ) ) *
825  hbar_h[j / 2][j % 2] );
826  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
827  //Matrix[j/2][j%2][i/2][i%2] *= G_F*alpha/4/sqrt(2)/EvtConst::pi*V_tb*conj(V_ts);
828  //std::cout << "Matrix = " << Matrix[j/2][j%2][i/2][i%2] << std::endl;
829  //std::cout << "--------------------------------------------------" << std::endl;
830  spins[0] = j % 2;
831  spins[1] = j / 2;
832  spins[2] = i / 2;
833  spins[3] = i % 2;
834  amp->vertex( spins, Matrix[j / 2][j % 2][i / 2][i % 2] );
835  }
836 
837  //std::cout << "==================================================" << std::endl;
838  //std::cout << "Lambda_b0: " << parent->getP4Restframe() << std::endl;
839  //std::cout << "Lambda0: " << parent->getDaug(0)->getP4() << std::endl;
840  //std::cout << "mu-: " << parent->getDaug(1)->getP4() << std::endl;
841  //std::cout << "mu+: " << parent->getDaug(2)->getP4() << std::endl;
842  //std::cout << "P_mu: " << P_mu << std::endl;
843  //std::cout << "q_mu: " << q_mu << std::endl;
844  //std::cout << "q2: " << q2 << std::endl;
845  //std::cout << "==================================================" << std::endl;
846 
847  return;
848 }
849 
850 EvtTensor4C EvtLb2Lll::EvtLeptonTG5Current( const EvtDiracSpinor& d,
851  const EvtDiracSpinor& dp )
852 {
853  // <u|sigma^munu*gamma^5|v>
854 
855  EvtTensor4C temp;
856  temp.zero();
857  EvtComplex i2( 0, 0.5 );
858 
859  static EvtGammaMatrix mat01 = EvtGammaMatrix::g0() *
860  ( EvtGammaMatrix::g0() * EvtGammaMatrix::g1() -
861  EvtGammaMatrix::g1() * EvtGammaMatrix::g0() ) *
862  EvtGammaMatrix::g5();
863  static EvtGammaMatrix mat02 = EvtGammaMatrix::g0() *
864  ( EvtGammaMatrix::g0() * EvtGammaMatrix::g2() -
865  EvtGammaMatrix::g2() * EvtGammaMatrix::g0() ) *
866  EvtGammaMatrix::g5();
867  static EvtGammaMatrix mat03 = EvtGammaMatrix::g0() *
868  ( EvtGammaMatrix::g0() * EvtGammaMatrix::g3() -
869  EvtGammaMatrix::g3() * EvtGammaMatrix::g0() ) *
870  EvtGammaMatrix::g5();
871  static EvtGammaMatrix mat12 = EvtGammaMatrix::g0() *
872  ( EvtGammaMatrix::g1() * EvtGammaMatrix::g2() -
873  EvtGammaMatrix::g2() * EvtGammaMatrix::g1() ) *
874  EvtGammaMatrix::g5();
875  static EvtGammaMatrix mat13 = EvtGammaMatrix::g0() *
876  ( EvtGammaMatrix::g1() * EvtGammaMatrix::g3() -
877  EvtGammaMatrix::g3() * EvtGammaMatrix::g1() ) *
878  EvtGammaMatrix::g5();
879  static EvtGammaMatrix mat23 = EvtGammaMatrix::g0() *
880  ( EvtGammaMatrix::g2() * EvtGammaMatrix::g3() -
881  EvtGammaMatrix::g3() * EvtGammaMatrix::g2() ) *
882  EvtGammaMatrix::g5();
883 
884  temp.set( 0, 1, i2 * ( d * ( mat01 * dp ) ) );
885  temp.set( 1, 0, -temp.get( 0, 1 ) );
886 
887  temp.set( 0, 2, i2 * ( d * ( mat02 * dp ) ) );
888  temp.set( 2, 0, -temp.get( 0, 2 ) );
889 
890  temp.set( 0, 3, i2 * ( d * ( mat03 * dp ) ) );
891  temp.set( 3, 0, -temp.get( 0, 3 ) );
892 
893  temp.set( 1, 2, i2 * ( d * ( mat12 * dp ) ) );
894  temp.set( 2, 1, -temp.get( 1, 2 ) );
895 
896  temp.set( 1, 3, i2 * ( d * ( mat13 * dp ) ) );
897  temp.set( 3, 1, -temp.get( 1, 3 ) );
898 
899  temp.set( 2, 3, i2 * ( d * ( mat23 * dp ) ) );
900  temp.set( 3, 2, -temp.get( 2, 3 ) );
901 
902  return temp;
903 }
EvtLb2Lll::m_noTries
long m_noTries
Definition: EvtLb2Lll.hh:51
EvtGenFunctions::directProd
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
Definition: EvtTensor3C.cpp:154
EvtWilsonCoefficients::Yld
EvtComplex Yld(double q2, double ki[], double Gi[], double Mi[], int ni, EvtComplex c1, EvtComplex c2, EvtComplex c3, EvtComplex c4, EvtComplex c5, EvtComplex c6, double ialpha)
Definition: EvtWilsonCoefficients.cpp:563
EvtTensor4C::set
void set(int i, int j, const EvtComplex &c)
Definition: EvtTensor4C.hh:107
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.cpp:382
EvtLb2Lll::m_HEPmodel
std::string m_HEPmodel
Definition: EvtLb2Lll.hh:56
EvtSpinDensity::setDiag
void setDiag(int n)
Definition: EvtSpinDensity.cpp:119
EvtParticle::sp
virtual EvtDiracSpinor sp(int) const
Definition: EvtParticle.cpp:677
EvtLeptonTCurrent
EvtTensor4C EvtLeptonTCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:272
EvtWilsonCoefficients::hzs
EvtComplex hzs(double z, double shat, double mu, double M_b)
Definition: EvtWilsonCoefficients.cpp:480
EvtDecayBase::getArg
double getArg(unsigned int j)
Definition: EvtDecayBase.cpp:604
EvtDecayBase::getArgStr
std::string getArgStr(int j) const
Definition: EvtDecayBase.hh:78
EvtGammaMatrix::g0
static const EvtGammaMatrix & g0()
Definition: EvtGammaMatrix.cpp:232
EvtLeptonPCurrent
EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:260
dual
EvtTensor4C dual(const EvtTensor4C &t2)
Definition: EvtTensor4C.cpp:364
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtLb2Lll::clone
EvtDecayBase * clone() override
Definition: EvtLb2Lll.cpp:36
cont
EvtComplex cont(const EvtTensor4C &t1, const EvtTensor4C &t2)
Definition: EvtTensor4C.cpp:291
EvtLb2Lll::m_polarizationIntroduction
std::string m_polarizationIntroduction
Definition: EvtLb2Lll.hh:55
EvtTensor4C::cont2
EvtVector4C cont2(const EvtVector4C &v4) const
Definition: EvtTensor4C.cpp:461
EvtLb2Lll::m_WC
EvtWilsonCoefficients m_WC
Definition: EvtLb2Lll.hh:60
EvtLb2Lll::decay
void decay(EvtParticle *p) override
Definition: EvtLb2Lll.cpp:377
EvtLb2Lll::m_poleSize
double m_poleSize
Definition: EvtLb2Lll.hh:50
EvtParticle::getSpinStates
int getSpinStates() const
Definition: EvtParticle.cpp:145
EvtLb2Lll::m_decayName
std::string m_decayName
Definition: EvtLb2Lll.hh:54
EvtLb2Lll::m_FFtype
std::string m_FFtype
Definition: EvtLb2Lll.hh:57
EvtParticle::makeDaughters
void makeDaughters(unsigned int ndaug, EvtId *id)
Definition: EvtParticle.cpp:1239
EvtLb2Lll::getName
std::string getName() override
Definition: EvtLb2Lll.cpp:41
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:167
boostTo
EvtRaritaSchwinger boostTo(const EvtRaritaSchwinger &rs, const EvtVector4R p4)
Definition: EvtRaritaSchwinger.cpp:35
EvtDiracParticle::init
void init(EvtId part_n, const EvtVector4R &p4) override
Definition: EvtDiracParticle.cpp:36
EvtWilsonCoefficients::C9efftilda
EvtComplex C9efftilda(double z, double shat, double alpha_S, EvtComplex c1, EvtComplex c2, EvtComplex c3, EvtComplex c4, EvtComplex c5, EvtComplex c6, EvtComplex c9tilda, int ksi)
Definition: EvtWilsonCoefficients.cpp:535
EvtAmp::getSpinDensity
EvtSpinDensity getSpinDensity()
Definition: EvtAmp.cpp:142
EvtParticle::spParent
virtual EvtDiracSpinor spParent(int) const
Definition: EvtParticle.cpp:665
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:66
EvtLb2Lll::m_maxProbability
double m_maxProbability
Definition: EvtLb2Lll.hh:49
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cpp:321
EvtLb2Lll::init
void init() override
Definition: EvtLb2Lll.cpp:46
EvtId
Definition: EvtId.hh:27
EvtLeptonACurrent
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:224
EvtAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
EvtConst::pi
static const double pi
Definition: EvtConst.hh:26
EvtLb2Lll::calcAmp
void calcAmp(EvtAmp *amp, EvtParticle *parent)
Definition: EvtLb2Lll.cpp:384
EvtDecayBase::getParentId
EvtId getParentId() const
Definition: EvtDecayBase.hh:61
EvtVector4R::mass2
double mass2() const
Definition: EvtVector4R.hh:100
EvtGammaMatrix::g1
static const EvtGammaMatrix & g1()
Definition: EvtGammaMatrix.cpp:257
EvtParticle::initializePhaseSpace
double initializePhaseSpace(unsigned int numdaughter, EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cpp:1102
EvtAmp
Definition: EvtAmp.hh:30
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cpp:150
EvtGammaMatrix::g2
static const EvtGammaMatrix & g2()
Definition: EvtGammaMatrix.cpp:281
lambda
double lambda(double q, double m1, double m2)
Definition: EvtFlatQ2.cpp:31
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtGammaMatrix::g3
static const EvtGammaMatrix & g3()
Definition: EvtGammaMatrix.cpp:305
EvtDecayBase::getNDaug
int getNDaug() const
Definition: EvtDecayBase.hh:65
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:240
EvtLb2Lll::EvtLeptonTG5Current
EvtTensor4C EvtLeptonTG5Current(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtLb2Lll.cpp:850
EvtLb2Lll::m_omega
double m_omega
Definition: EvtLb2Lll.hh:52
EvtLb2Lll::m_effectContribution
std::string m_effectContribution
Definition: EvtLb2Lll.hh:58
EvtLb2Lll::m_polarizationLambdab0
double m_polarizationLambdab0
Definition: EvtLb2Lll.hh:48
EvtLeptonVACurrent
EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:179
EvtParticle::setDiagonalSpinDensity
void setDiagonalSpinDensity()
Definition: EvtParticle.cpp:170
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtLeptonVCurrent
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:209
EvtParticle::noLifeTime
void noLifeTime()
Definition: EvtParticle.hh:382
EvtDecayAmp::_amp2
EvtAmp _amp2
Definition: EvtDecayAmp.hh:73
EvtWilsonCoefficients::C7b2sg
EvtComplex C7b2sg(double alpha_S, double et, EvtComplex c2, double M_t, double M_W)
Definition: EvtWilsonCoefficients.cpp:552
EvtPDL::getMass
static double getMass(EvtId i)
Definition: EvtPDL.cpp:319
EvtLeptonSCurrent
EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:248
EvtLb2Lll::initProbMax
void initProbMax() override
Definition: EvtLb2Lll.cpp:282
EvtTensor4C::get
const EvtComplex & get(int i, int j) const
Definition: EvtTensor4C.hh:112
EvtParticle::getP4Restframe
EvtVector4R getP4Restframe() const
Definition: EvtParticle.cpp:793
EvtAmp::init
void init(EvtId p, int ndaug, EvtId *daug)
Definition: EvtAmp.cpp:68
EvtDecayBase::getNArg
int getNArg() const
Definition: EvtDecayBase.hh:68
EvtGammaMatrix::g5
static const EvtGammaMatrix & g5()
Definition: EvtGammaMatrix.cpp:329
EvtSpinDensity::normalizedProb
double normalizedProb(const EvtSpinDensity &d)
Definition: EvtSpinDensity.cpp:132
EvtDecayBase::getDaug
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67
Generated by   doxygen 1.8.15