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.
EvtSLBaryonAmp.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 
22 
23 #include "EvtGenBase/EvtAmp.hh"
27 #include "EvtGenBase/EvtGenKine.hh"
28 #include "EvtGenBase/EvtId.hh"
29 #include "EvtGenBase/EvtPDL.hh"
31 #include "EvtGenBase/EvtPatches.hh"
33 #include "EvtGenBase/EvtReport.hh"
37 
38 #include <stdlib.h>
39 
40 using std::endl;
41 
43 {
44 }
45 
47  EvtSemiLeptonicFF* FormFactors )
48 {
49  static EvtId EM = EvtPDL::getId( "e-" );
50  static EvtId MUM = EvtPDL::getId( "mu-" );
51  static EvtId TAUM = EvtPDL::getId( "tau-" );
52  static EvtId EP = EvtPDL::getId( "e+" );
53  static EvtId MUP = EvtPDL::getId( "mu+" );
54  static EvtId TAUP = EvtPDL::getId( "tau+" );
55 
56  //Add the lepton and neutrino 4 momenta to find q2
57 
58  EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();
59  double q2 = ( q.mass2() );
60 
61  double f1v, f1a, f2v, f2a;
62  double m_meson = parent->getDaug( 0 )->mass();
63 
64  FormFactors->getbaryonff( parent->getId(), parent->getDaug( 0 )->getId(),
65  q2, m_meson, &f1v, &f1a, &f2v, &f2a );
66 
67  EvtVector4R p4b;
68  p4b.set( parent->mass(), 0.0, 0.0, 0.0 );
69 
70  EvtVector4C temp_00_term1;
71  EvtVector4C temp_00_term2;
72 
73  EvtVector4C temp_01_term1;
74  EvtVector4C temp_01_term2;
75 
76  EvtVector4C temp_10_term1;
77  EvtVector4C temp_10_term2;
78 
79  EvtVector4C temp_11_term1;
80  EvtVector4C temp_11_term2;
81 
82  EvtDiracSpinor p0 = parent->sp( 0 );
83  EvtDiracSpinor p1 = parent->sp( 1 );
84 
85  EvtDiracSpinor d0 = parent->getDaug( 0 )->spParent( 0 );
86  EvtDiracSpinor d1 = parent->getDaug( 0 )->spParent( 1 );
87 
88  temp_00_term1.set( 0, f1v * ( d0 * ( EvtGammaMatrix::g0() * p0 ) ) );
89  temp_00_term2.set(
90  0,
91  f1a * ( d0 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p0 ) ) );
92  temp_01_term1.set( 0, f1v * ( d0 * ( EvtGammaMatrix::g0() * p1 ) ) );
93  temp_01_term2.set(
94  0,
95  f1a * ( d0 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p1 ) ) );
96  temp_10_term1.set( 0, f1v * ( d1 * ( EvtGammaMatrix::g0() * p0 ) ) );
97  temp_10_term2.set(
98  0,
99  f1a * ( d1 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p0 ) ) );
100  temp_11_term1.set( 0, f1v * ( d1 * ( EvtGammaMatrix::g0() * p1 ) ) );
101  temp_11_term2.set(
102  0,
103  f1a * ( d1 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p1 ) ) );
104 
105  temp_00_term1.set( 1, f1v * ( d0 * ( EvtGammaMatrix::g1() * p0 ) ) );
106  temp_00_term2.set(
107  1,
108  f1a * ( d0 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p0 ) ) );
109  temp_01_term1.set( 1, f1v * ( d0 * ( EvtGammaMatrix::g1() * p1 ) ) );
110  temp_01_term2.set(
111  1,
112  f1a * ( d0 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p1 ) ) );
113  temp_10_term1.set( 1, f1v * ( d1 * ( EvtGammaMatrix::g1() * p0 ) ) );
114  temp_10_term2.set(
115  1,
116  f1a * ( d1 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p0 ) ) );
117  temp_11_term1.set( 1, f1v * ( d1 * ( EvtGammaMatrix::g1() * p1 ) ) );
118  temp_11_term2.set(
119  1,
120  f1a * ( d1 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p1 ) ) );
121 
122  temp_00_term1.set( 2, f1v * ( d0 * ( EvtGammaMatrix::g2() * p0 ) ) );
123  temp_00_term2.set(
124  2,
125  f1a * ( d0 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p0 ) ) );
126  temp_01_term1.set( 2, f1v * ( d0 * ( EvtGammaMatrix::g2() * p1 ) ) );
127  temp_01_term2.set(
128  2,
129  f1a * ( d0 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p1 ) ) );
130  temp_10_term1.set( 2, f1v * ( d1 * ( EvtGammaMatrix::g2() * p0 ) ) );
131  temp_10_term2.set(
132  2,
133  f1a * ( d1 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p0 ) ) );
134  temp_11_term1.set( 2, f1v * ( d1 * ( EvtGammaMatrix::g2() * p1 ) ) );
135  temp_11_term2.set(
136  2,
137  f1a * ( d1 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p1 ) ) );
138 
139  temp_00_term1.set( 3, f1v * ( d0 * ( EvtGammaMatrix::g3() * p0 ) ) );
140  temp_00_term2.set(
141  3,
142  f1a * ( d0 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p0 ) ) );
143  temp_01_term1.set( 3, f1v * ( d0 * ( EvtGammaMatrix::g3() * p1 ) ) );
144  temp_01_term2.set(
145  3,
146  f1a * ( d0 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p1 ) ) );
147  temp_10_term1.set( 3, f1v * ( d1 * ( EvtGammaMatrix::g3() * p0 ) ) );
148  temp_10_term2.set(
149  3,
150  f1a * ( d1 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p0 ) ) );
151  temp_11_term1.set( 3, f1v * ( d1 * ( EvtGammaMatrix::g3() * p1 ) ) );
152  temp_11_term2.set(
153  3,
154  f1a * ( d1 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p1 ) ) );
155 
156  EvtVector4C l1, l2;
157 
158  EvtId l_num = parent->getDaug( 1 )->getId();
159  if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
160  l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
161  parent->getDaug( 2 )->spParentNeutrino() );
162  l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
163  parent->getDaug( 2 )->spParentNeutrino() );
164  } else {
165  if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
166  l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
167  parent->getDaug( 1 )->spParent( 0 ) );
168  l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
169  parent->getDaug( 1 )->spParent( 1 ) );
170  } else {
171  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
172  << "Wrong lepton number" << endl;
173  }
174  }
175 
176  amp.vertex( 0, 0, 0, l1.cont( temp_00_term1 + temp_00_term2 ) );
177  amp.vertex( 0, 0, 1, l2.cont( temp_00_term1 + temp_00_term2 ) );
178 
179  amp.vertex( 0, 1, 0, l1.cont( temp_01_term1 + temp_01_term2 ) );
180  amp.vertex( 0, 1, 1, l2.cont( temp_01_term1 + temp_01_term2 ) );
181 
182  amp.vertex( 1, 0, 0, l1.cont( temp_10_term1 + temp_10_term2 ) );
183  amp.vertex( 1, 0, 1, l2.cont( temp_10_term1 + temp_10_term2 ) );
184 
185  amp.vertex( 1, 1, 0, l1.cont( temp_11_term1 + temp_11_term2 ) );
186  amp.vertex( 1, 1, 1, l2.cont( temp_11_term1 + temp_11_term2 ) );
187 
188  return;
189 }
190 
191 double EvtSLBaryonAmp::CalcMaxProb( EvtId parent, EvtId baryon, EvtId lepton,
192  EvtId nudaug, EvtSemiLeptonicFF* FormFactors,
193  EvtComplex r00, EvtComplex r01,
194  EvtComplex r10, EvtComplex r11 )
195 {
196  //This routine takes the arguements parent, baryon, and lepton
197  //number, and a form factor model, and returns a maximum
198  //probability for this semileptonic form factor model. A
199  //brute force method is used. The 2D cos theta lepton and
200  //q2 phase space is probed.
201 
202  //Start by declaring a particle at rest.
203 
204  //It only makes sense to have a scalar parent. For now.
205  //This should be generalized later.
206 
207  // EvtScalarParticle *scalar_part;
208  // scalar_part=new EvtScalarParticle;
209 
210  EvtDiracParticle* dirac_part;
211  EvtParticle* root_part;
212  dirac_part = new EvtDiracParticle;
213 
214  //cludge to avoid generating random numbers!
215  // scalar_part->noLifeTime();
216  dirac_part->noLifeTime();
217 
218  EvtVector4R p_init;
219 
220  p_init.set( EvtPDL::getMass( parent ), 0.0, 0.0, 0.0 );
221  // scalar_part->init(parent,p_init);
222  // root_part=(EvtParticle *)scalar_part;
223  // root_part->set_type(EvtSpinType::SCALAR);
224 
225  dirac_part->init( parent, p_init );
226  root_part = (EvtParticle*)dirac_part;
227  root_part->setDiagonalSpinDensity();
228 
229  EvtParticle *daughter, *lep, *trino;
230 
231  EvtAmp amp;
232 
233  EvtId listdaug[3];
234  listdaug[0] = baryon;
235  listdaug[1] = lepton;
236  listdaug[2] = nudaug;
237 
238  amp.init( parent, 3, listdaug );
239 
240  root_part->makeDaughters( 3, listdaug );
241  daughter = root_part->getDaug( 0 );
242  lep = root_part->getDaug( 1 );
243  trino = root_part->getDaug( 2 );
244 
245  //cludge to avoid generating random numbers!
246  daughter->noLifeTime();
247  lep->noLifeTime();
248  trino->noLifeTime();
249 
250  //Initial particle is unpolarized, well it is a scalar so it is
251  //trivial
252  EvtSpinDensity rho;
253  rho.setDiag( root_part->getSpinStates() );
254 
255  double mass[3];
256 
257  double m = root_part->mass();
258 
259  EvtVector4R p4baryon, p4lepton, p4nu, p4w;
260  double q2max;
261 
262  double q2, elepton, plepton;
263  int i, j;
264  double erho, prho, costl;
265 
266  double maxfoundprob = 0.0;
267  double prob = -10.0;
268  int massiter;
269 
270  for ( massiter = 0; massiter < 3; massiter++ ) {
271  mass[0] = EvtPDL::getMass( baryon );
272  mass[1] = EvtPDL::getMass( lepton );
273  mass[2] = EvtPDL::getMass( nudaug );
274  if ( massiter == 1 ) {
275  mass[0] = EvtPDL::getMinMass( baryon );
276  }
277  if ( massiter == 2 ) {
278  mass[0] = EvtPDL::getMaxMass( baryon );
279  }
280 
281  q2max = ( m - mass[0] ) * ( m - mass[0] );
282 
283  //loop over q2
284 
285  for ( i = 0; i < 25; i++ ) {
286  q2 = ( ( i + 0.5 ) * q2max ) / 25.0;
287 
288  erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m );
289 
290  prho = sqrt( erho * erho - mass[0] * mass[0] );
291 
292  p4baryon.set( erho, 0.0, 0.0, -1.0 * prho );
293  p4w.set( m - erho, 0.0, 0.0, prho );
294 
295  //This is in the W rest frame
296  elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) );
297  plepton = sqrt( elepton * elepton - mass[1] * mass[1] );
298 
299  double probctl[3];
300 
301  for ( j = 0; j < 3; j++ ) {
302  costl = 0.99 * ( j - 1.0 );
303 
304  //These are in the W rest frame. Need to boost out into
305  //the B frame.
306  p4lepton.set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ),
307  plepton * costl );
308  p4nu.set( plepton, 0.0,
309  -1.0 * plepton * sqrt( 1.0 - costl * costl ),
310  -1.0 * plepton * costl );
311 
312  EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho );
313  p4lepton = boostTo( p4lepton, boost );
314  p4nu = boostTo( p4nu, boost );
315 
316  //Now initialize the daughters...
317 
318  daughter->init( baryon, p4baryon );
319  lep->init( lepton, p4lepton );
320  trino->init( nudaug, p4nu );
321 
322  CalcAmp( root_part, amp, FormFactors, r00, r01, r10, r11 );
323 
324  //Now find the probability at this q2 and cos theta lepton point
325  //and compare to maxfoundprob.
326 
327  //Do a little magic to get the probability!!
328  prob = rho.normalizedProb( amp.getSpinDensity() );
329 
330  probctl[j] = prob;
331  }
332 
333  //probclt contains prob at ctl=-1,0,1.
334  //prob=a+b*ctl+c*ctl^2
335 
336  double a = probctl[1];
337  double b = 0.5 * ( probctl[2] - probctl[0] );
338  double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1];
339 
340  prob = probctl[0];
341  if ( probctl[1] > prob )
342  prob = probctl[1];
343  if ( probctl[2] > prob )
344  prob = probctl[2];
345 
346  if ( fabs( c ) > 1e-20 ) {
347  double ctlx = -0.5 * b / c;
348  if ( fabs( ctlx ) < 1.0 ) {
349  double probtmp = a + b * ctlx + c * ctlx * ctlx;
350  if ( probtmp > prob )
351  prob = probtmp;
352  }
353  }
354 
355  //EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
356  // << probctl[0]<<" "
357  // << probctl[1]<<" "
358  // << probctl[2]<<std::endl;
359 
360  if ( prob > maxfoundprob ) {
361  maxfoundprob = prob;
362  }
363  }
364  if ( EvtPDL::getWidth( baryon ) <= 0.0 ) {
365  //if the particle is narrow dont bother with changing the mass.
366  massiter = 4;
367  }
368  }
369  root_part->deleteTree();
370 
371  maxfoundprob *= 1.1;
372  return maxfoundprob;
373 }
375  EvtSemiLeptonicFF* FormFactors, EvtComplex r00,
376  EvtComplex r01, EvtComplex r10, EvtComplex r11 )
377 {
378  // Leptons
379  static EvtId EM = EvtPDL::getId( "e-" );
380  static EvtId MUM = EvtPDL::getId( "mu-" );
381  static EvtId TAUM = EvtPDL::getId( "tau-" );
382  // Anti-Leptons
383  static EvtId EP = EvtPDL::getId( "e+" );
384  static EvtId MUP = EvtPDL::getId( "mu+" );
385  static EvtId TAUP = EvtPDL::getId( "tau+" );
386 
387  // Baryons
388  static EvtId LAMCP = EvtPDL::getId( "Lambda_c+" );
389  static EvtId LAMC1P = EvtPDL::getId( "Lambda_c(2593)+" );
390  static EvtId LAMC2P = EvtPDL::getId( "Lambda_c(2625)+" );
391  static EvtId LAMB = EvtPDL::getId( "Lambda_b0" );
392  static EvtId PRO = EvtPDL::getId( "p+" );
393  static EvtId N1440 = EvtPDL::getId( "N(1440)+" );
394  static EvtId N1520 = EvtPDL::getId( "N(1520)+" );
395  static EvtId N1535 = EvtPDL::getId( "N(1535)+" );
396  static EvtId N1720 = EvtPDL::getId( "N(1720)+" );
397  static EvtId N1650 = EvtPDL::getId( "N(1650)+" );
398  static EvtId N1700 = EvtPDL::getId( "N(1700)+" );
399  static EvtId N1710 = EvtPDL::getId( "N(1710)+" );
400  static EvtId N1875 = EvtPDL::getId( "N(1875)+" );
401  static EvtId N1900 = EvtPDL::getId( "N(1900)+" );
402 
403  // Anti-Baryons
404  static EvtId LAMCM = EvtPDL::getId( "anti-Lambda_c-" );
405  static EvtId LAMC1M = EvtPDL::getId( "anti-Lambda_c(2593)-" );
406  static EvtId LAMC2M = EvtPDL::getId( "anti-Lambda_c(2625)-" );
407  static EvtId LAMBB = EvtPDL::getId( "anti-Lambda_b0" );
408  static EvtId PROB = EvtPDL::getId( "anti-p-" );
409  static EvtId N1440B = EvtPDL::getId( "anti-N(1440)-" );
410  static EvtId N1520B = EvtPDL::getId( "anti-N(1520)-" );
411  static EvtId N1535B = EvtPDL::getId( "anti-N(1535)-" );
412  static EvtId N1720B = EvtPDL::getId( "anti-N(1720)-" );
413  static EvtId N1650B = EvtPDL::getId( "anti-N(1650)-" );
414  static EvtId N1700B = EvtPDL::getId( "anti-N(1700)-" );
415  static EvtId N1710B = EvtPDL::getId( "anti-N(1710)-" );
416  static EvtId N1875B = EvtPDL::getId( "anti-N(1875)-" );
417  static EvtId N1900B = EvtPDL::getId( "anti-N(1900)-" );
418 
419  // Set the spin density matrix of the parent baryon
420  EvtSpinDensity rho;
421  rho.setDim( 2 );
422  rho.set( 0, 0, r00 );
423  rho.set( 0, 1, r01 );
424 
425  rho.set( 1, 0, r10 );
426  rho.set( 1, 1, r11 );
427 
428  EvtVector4R vector4P = parent->getP4Lab();
429  double pmag = vector4P.d3mag();
430  double cosTheta = vector4P.get( 3 ) / pmag;
431 
432  double theta = acos( cosTheta );
433  double phi = atan2( vector4P.get( 2 ), vector4P.get( 1 ) );
434 
435  parent->setSpinDensityForwardHelicityBasis( rho, phi, theta, 0.0 );
436  //parent->setSpinDensityForward(rho);
437 
438  // Set the four momentum of the parent baryon in it's rest frame
439  EvtVector4R p4b;
440  p4b.set( parent->mass(), 0.0, 0.0, 0.0 );
441 
442  // Get the four momentum of the daughter baryon in the parent's rest frame
443  EvtVector4R p4daught = parent->getDaug( 0 )->getP4();
444 
445  // Add the lepton and neutrino 4 momenta to find q (q^2)
446  EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();
447 
448  double q2 = q.mass2();
449 
450  EvtId l_num = parent->getDaug( 1 )->getId();
451  EvtId bar_num = parent->getDaug( 0 )->getId();
452  EvtId par_num = parent->getId();
453 
454  double baryonmass = parent->getDaug( 0 )->mass();
455 
456  // Handle spin-1/2 daughter baryon Dirac spinor cases
457  if ( EvtPDL::getSpinType( parent->getDaug( 0 )->getId() ) ==
459  // Set the form factors
460  double f1, f2, f3, g1, g2, g3;
461  FormFactors->getdiracff( par_num, bar_num, q2, baryonmass, &f1, &f2,
462  &f3, &g1, &g2, &g3 );
463 
464  const double form_fact[6] = {f1, f2, f3, g1, g2, g3};
465 
466  EvtVector4C b11, b12, b21, b22, l1, l2;
467 
468  // Lepton Current
469  if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
470  l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
471  parent->getDaug( 2 )->spParentNeutrino() );
472  l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
473  parent->getDaug( 2 )->spParentNeutrino() );
474 
475  } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
476  l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
477  parent->getDaug( 1 )->spParent( 0 ) );
478  l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
479  parent->getDaug( 1 )->spParent( 1 ) );
480 
481  } else {
482  EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number \n";
483  ::abort();
484  }
485 
486  // Baryon current
487 
488  // Flag for particle/anti-particle parent, daughter with same/opp. parity
489  // pflag = 0 => particle, same parity parent, daughter
490  // pflag = 1 => particle, opp. parity parent, daughter
491  // pflag = 2 => anti-particle, same parity parent, daughter
492  // pflag = 3 => anti-particle, opp. parity parent, daughter
493 
494  int pflag = 0;
495 
496  // Handle 1/2+ -> 1/2+ first
497  if ( ( par_num == LAMB && bar_num == LAMCP ) ||
498  ( par_num == LAMBB && bar_num == LAMCM ) ||
499  ( par_num == LAMB && bar_num == PRO ) ||
500  ( par_num == LAMBB && bar_num == PROB ) ||
501  ( par_num == LAMB && bar_num == N1440 ) ||
502  ( par_num == LAMBB && bar_num == N1440B ) ||
503  ( par_num == LAMB && bar_num == N1710 ) ||
504  ( par_num == LAMBB && bar_num == N1710B )
505 
506  ) {
507  // Set particle/anti-particle flag
508  if ( bar_num == LAMCP || bar_num == PRO || bar_num == N1440 ||
509  bar_num == N1710 )
510  pflag = 0;
511  else if ( bar_num == LAMCM || bar_num == PROB ||
512  bar_num == N1440B || bar_num == N1710B )
513  pflag = 2;
514 
515  b11 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 0 ),
516  parent->sp( 0 ), p4b, p4daught, form_fact,
517  pflag );
518  b21 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 0 ),
519  parent->sp( 1 ), p4b, p4daught, form_fact,
520  pflag );
521  b12 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 1 ),
522  parent->sp( 0 ), p4b, p4daught, form_fact,
523  pflag );
524  b22 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 1 ),
525  parent->sp( 1 ), p4b, p4daught, form_fact,
526  pflag );
527  }
528 
529  // Handle 1/2+ -> 1/2- second
530  else if ( ( par_num == LAMB && bar_num == LAMC1P ) ||
531  ( par_num == LAMBB && bar_num == LAMC1M ) ||
532  ( par_num == LAMB && bar_num == N1535 ) ||
533  ( par_num == LAMBB && bar_num == N1535B ) ||
534  ( par_num == LAMB && bar_num == N1650 ) ||
535  ( par_num == LAMBB && bar_num == N1650B ) ) {
536  // Set particle/anti-particle flag
537  if ( bar_num == LAMC1P || bar_num == N1535 || bar_num == N1650 )
538  pflag = 1;
539  else if ( bar_num == LAMC1M || bar_num == N1535B || bar_num == N1650B )
540  pflag = 3;
541 
542  b11 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 0 ) ),
543  ( EvtGammaMatrix::g5() * parent->sp( 0 ) ),
544  p4b, p4daught, form_fact, pflag );
545  b21 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 0 ) ),
546  ( EvtGammaMatrix::g5() * parent->sp( 1 ) ),
547  p4b, p4daught, form_fact, pflag );
548  b12 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 1 ) ),
549  ( EvtGammaMatrix::g5() * parent->sp( 0 ) ),
550  p4b, p4daught, form_fact, pflag );
551  b22 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 1 ) ),
552  ( EvtGammaMatrix::g5() * parent->sp( 1 ) ),
553  p4b, p4daught, form_fact, pflag );
554 
555  }
556 
557  else {
558  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
559  << "Dirac semilep. baryon current "
560  << "not implemented for this decay sequence." << std::endl;
561  ::abort();
562  }
563 
564  amp.vertex( 0, 0, 0, l1 * b11 );
565  amp.vertex( 0, 0, 1, l2 * b11 );
566 
567  amp.vertex( 1, 0, 0, l1 * b21 );
568  amp.vertex( 1, 0, 1, l2 * b21 );
569 
570  amp.vertex( 0, 1, 0, l1 * b12 );
571  amp.vertex( 0, 1, 1, l2 * b12 );
572 
573  amp.vertex( 1, 1, 0, l1 * b22 );
574  amp.vertex( 1, 1, 1, l2 * b22 );
575 
576  }
577 
578  // Need special handling for the spin-3/2 daughter baryon
579  // Rarita-Schwinger spinor cases
580  else if ( EvtPDL::getSpinType( parent->getDaug( 0 )->getId() ) ==
582  // Set the form factors
583  double f1, f2, f3, f4, g1, g2, g3, g4;
584  FormFactors->getraritaff( par_num, bar_num, q2, baryonmass, &f1, &f2,
585  &f3, &f4, &g1, &g2, &g3, &g4 );
586 
587  const double form_fact[8] = {f1, f2, f3, f4, g1, g2, g3, g4};
588 
589  EvtId l_num = parent->getDaug( 1 )->getId();
590 
591  EvtVector4C b11, b12, b21, b22, b13, b23, b14, b24, l1, l2;
592 
593  // Lepton Current
594  if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
595  // Lepton Current
596  l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
597  parent->getDaug( 2 )->spParentNeutrino() );
598  l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
599  parent->getDaug( 2 )->spParentNeutrino() );
600  } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
601  l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
602  parent->getDaug( 1 )->spParent( 0 ) );
603  l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
604  parent->getDaug( 1 )->spParent( 1 ) );
605 
606  } else {
607  EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number \n";
608  }
609 
610  // Baryon Current
611  // Declare particle, anti-particle flag, same/opp. parity
612  // pflag = 0 => particle
613  // pflag = 1 => anti-particle
614  int pflag = 0;
615 
616  // Handle cases of 1/2+ -> 3/2- or 3/2+
617  if ( ( par_num == LAMB && bar_num == LAMC2P ) ||
618  ( par_num == LAMB && bar_num == N1720 ) ||
619  ( par_num == LAMB && bar_num == N1520 ) ||
620  ( par_num == LAMB && bar_num == N1700 ) ||
621  ( par_num == LAMB && bar_num == N1875 ) ||
622  ( par_num == LAMB && bar_num == N1900 ) ) {
623  // Set flag for particle case
624  pflag = 0;
625  } else if ( ( par_num == LAMBB && bar_num == LAMC2M ) ||
626  ( par_num == LAMBB && bar_num == N1520B ) ||
627  ( par_num == LAMBB && bar_num == N1700B ) ||
628  ( par_num == LAMBB && bar_num == N1875B ) ) {
629  // Set flag for anti-particle opposite parity case
630  pflag = 1;
631  }
632  // Handle anti-particle case for 1/2+ -> 3/2+
633  else if ( ( par_num == LAMBB && bar_num == N1720B ) ||
634  ( par_num == LAMBB && bar_num == N1900B ) ) {
635  pflag = 2;
636  } else {
637  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
638  << "Rarita-Schwinger semilep. baryon current "
639  << "not implemented for this decay sequence." << std::endl;
640  ::abort();
641  }
642 
643  // Baryon current
644  b11 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 0 ),
645  parent->sp( 0 ), p4b, p4daught,
646  form_fact, pflag );
647  b21 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 0 ),
648  parent->sp( 1 ), p4b, p4daught,
649  form_fact, pflag );
650 
651  b12 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 1 ),
652  parent->sp( 0 ), p4b, p4daught,
653  form_fact, pflag );
654  b22 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 1 ),
655  parent->sp( 1 ), p4b, p4daught,
656  form_fact, pflag );
657 
658  b13 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 2 ),
659  parent->sp( 0 ), p4b, p4daught,
660  form_fact, pflag );
661  b23 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 2 ),
662  parent->sp( 1 ), p4b, p4daught,
663  form_fact, pflag );
664 
665  b14 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 3 ),
666  parent->sp( 0 ), p4b, p4daught,
667  form_fact, pflag );
668  b24 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 3 ),
669  parent->sp( 1 ), p4b, p4daught,
670  form_fact, pflag );
671 
672  amp.vertex( 0, 0, 0, l1 * b11 );
673  amp.vertex( 0, 0, 1, l2 * b11 );
674 
675  amp.vertex( 1, 0, 0, l1 * b21 );
676  amp.vertex( 1, 0, 1, l2 * b21 );
677 
678  amp.vertex( 0, 1, 0, l1 * b12 );
679  amp.vertex( 0, 1, 1, l2 * b12 );
680 
681  amp.vertex( 1, 1, 0, l1 * b22 );
682  amp.vertex( 1, 1, 1, l2 * b22 );
683 
684  amp.vertex( 0, 2, 0, l1 * b13 );
685  amp.vertex( 0, 2, 1, l2 * b13 );
686 
687  amp.vertex( 1, 2, 0, l1 * b23 );
688  amp.vertex( 1, 2, 1, l2 * b23 );
689 
690  amp.vertex( 0, 3, 0, l1 * b14 );
691  amp.vertex( 0, 3, 1, l2 * b14 );
692 
693  amp.vertex( 1, 3, 0, l1 * b24 );
694  amp.vertex( 1, 3, 1, l2 * b24 );
695  }
696 }
697 
699  const EvtDiracSpinor& Bi,
700  EvtVector4R parent,
701  EvtVector4R daught,
702  const double* ff, int pflag )
703 {
704  // flag == 0 => particle
705  // flag == 1 => particle, opposite parity
706  // flag == 2 => anti-particle, same parity
707  // flag == 3 => anti-particle, opposite parity
708 
709  // particle
710  EvtComplex cv = EvtComplex( 1.0, 0. );
711  EvtComplex ca = EvtComplex( 1.0, 0. );
712 
713  EvtComplex cg0 = EvtComplex( 1.0, 0. );
714  EvtComplex cg5 = EvtComplex( 1.0, 0. );
715 
716  // antiparticle- same parity parent & daughter
717  if ( pflag == 2 ) {
718  cv = EvtComplex( -1.0, 0. );
719  ca = EvtComplex( 1.0, 0. );
720 
721  cg0 = EvtComplex( 1.0, 0.0 );
722  // Changed cg5 from -i to -1 as appears to fix particle - anti-particle discrepency
723  cg5 = EvtComplex( -1.0, 0.0 );
724  }
725  // antiparticle- opposite parity parent & daughter
726  else if ( pflag == 3 ) {
727  cv = EvtComplex( 1.0, 0. );
728  ca = EvtComplex( -1.0, 0. );
729 
730  // Changed cg0 from -i to -1 as appears to fix particle - anti-particle discrepency
731  cg0 = EvtComplex( -1.0, 0.0 );
732  cg5 = EvtComplex( 1.0, 0.0 );
733  }
734 
735  EvtVector4C t[6];
736 
737  // Term 1 = \bar{u}(p',s')*(F_1(q^2)*\gamma_{mu})*u(p,s)
738  t[0] = cv * EvtLeptonVCurrent( Bf, Bi );
739 
740  // Term 2 = \bar{u}(p',s')*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
741  t[1] = cg0 * EvtLeptonSCurrent( Bf, Bi ) * ( parent / parent.mass() );
742 
743  // Term 3 = \bar{u}(p',s')*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
744  t[2] = cg0 * EvtLeptonSCurrent( Bf, Bi ) * ( daught / daught.mass() );
745 
746  // Term 4 = \bar{u}(p',s')*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
747  t[3] = ca * EvtLeptonACurrent( Bf, Bi );
748 
749  // Term 5 = \bar{u}(p',s')*(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
750  t[4] = cg5 * EvtLeptonPCurrent( Bf, Bi ) * ( parent / parent.mass() );
751 
752  // Term 6 = \bar{u}(p',s')*(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
753  t[5] = cg5 * EvtLeptonPCurrent( Bf, Bi ) * ( daught / daught.mass() );
754 
755  // Sum the individual terms
756  EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] -
757  ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] );
758 
759  return current;
760 }
761 
763  const EvtRaritaSchwinger& Bf, const EvtDiracSpinor& Bi, EvtVector4R parent,
764  EvtVector4R daught, const double* ff, int pflag )
765 {
766  // flag == 0 => particle
767  // flag == 1 => anti-particle
768 
769  // particle
770  EvtComplex cv = EvtComplex( 1.0, 0. );
771  EvtComplex ca = EvtComplex( 1.0, 0. );
772 
773  EvtComplex cg0 = EvtComplex( 1.0, 0. );
774  EvtComplex cg5 = EvtComplex( 1.0, 0. );
775 
776  // antiparticle opposite parity
777  if ( pflag == 1 ) {
778  cv = EvtComplex( -1.0, 0. );
779  ca = EvtComplex( 1.0, 0. );
780 
781  cg0 = EvtComplex( 1.0, 0.0 );
782  cg5 = EvtComplex( -1.0, 0.0 );
783  }
784  // antiparticle same parity
785  else if ( pflag == 2 ) {
786  cv = EvtComplex( 1.0, 0. );
787  ca = EvtComplex( -1.0, 0. );
788 
789  cg0 = EvtComplex( -1.0, 0.0 );
790  cg5 = EvtComplex( 1.0, 0.0 );
791  }
792 
793  EvtVector4C t[8];
794  EvtTensor4C id;
795  id.setdiag( 1.0, 1.0, 1.0, 1.0 );
796 
797  EvtDiracSpinor tmp;
798  for ( int i = 0; i < 4; i++ ) {
799  tmp.set_spinor( i, Bf.getVector( i ) * parent );
800  }
801 
802  EvtVector4C v1, v2;
803  for ( int i = 0; i < 4; i++ ) {
804  v1.set( i, EvtLeptonSCurrent( Bf.getSpinor( i ), Bi ) );
805  v2.set( i, EvtLeptonPCurrent( Bf.getSpinor( i ), Bi ) );
806  }
807 
808  // Term 1 = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_1(q^2)*\gamma_{mu})*u(p,s)
809  t[0] = ( cv / parent.mass() ) * EvtLeptonVCurrent( tmp, Bi );
810 
811  // Term 2
812  // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
813  t[1] = ( ( cg0 / parent.mass() ) * EvtLeptonSCurrent( tmp, Bi ) ) *
814  ( parent / parent.mass() );
815 
816  // Term 3
817  // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
818  t[2] = ( ( cg0 / parent.mass() ) * EvtLeptonSCurrent( tmp, Bi ) ) *
819  ( daught / daught.mass() );
820 
821  // Term 4 = \bar{u}^{\alpha}(p',s')*(F_4(q^2)*g_{\alpha,\mu})*u(p,s)
822  t[3] = cg0 * ( id.cont2( v1 ) );
823 
824  // Term 5
825  // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
826  t[4] = ( ca / parent.mass() ) * EvtLeptonACurrent( tmp, Bi );
827 
828  // Term 6
829  // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
830  // *(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
831  t[5] = ( ( cg5 / parent.mass() ) * EvtLeptonPCurrent( tmp, Bi ) ) *
832  ( parent / parent.mass() );
833 
834  // Term 7
835  // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
836  // *(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
837  t[6] = ( ( cg5 / parent.mass() ) * EvtLeptonPCurrent( tmp, Bi ) ) *
838  ( daught / daught.mass() );
839 
840  // Term 8 = \bar{u}^{\alpha}(p',s')*(G_4(q^2)*g_{\alpha,\mu}*\gamma_5))*u(p,s)
841  t[7] = cg5 * ( id.cont2( v2 ) );
842 
843  // Sum the individual terms
844  EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] +
845  ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] -
846  ff[6] * t[6] - ff[7] * t[7] );
847 
848  return current;
849 }
virtual EvtRaritaSchwinger spRSParent(int) const
void setdiag(double t00, double t11, double t22, double t33)
virtual void getraritaff(EvtId parent, EvtId daught, double q2, double mass, double *f1, double *f2, double *f3, double *f4, double *g1, double *g2, double *g3, double *g4)=0
void setDiag(int n)
virtual EvtDiracSpinor sp(int) const
static EvtSpinType::spintype getSpinType(EvtId i)
Definition: EvtPDL.cpp:377
EvtVector4R getP4Lab() const
void set(int, const EvtComplex &)
Definition: EvtVector4C.hh:98
static const EvtGammaMatrix & g0()
EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
double mass() const
Definition: EvtVector4R.cpp:49
void set_spinor(int i, const EvtComplex &sp)
int getSpinStates() const
EvtId getId() const
void makeDaughters(unsigned int ndaug, EvtId *id)
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtRaritaSchwinger boostTo(const EvtRaritaSchwinger &rs, const EvtVector4R p4)
void init(EvtId part_n, const EvtVector4R &p4) override
EvtSpinDensity getSpinDensity()
Definition: EvtAmp.cpp:142
virtual EvtDiracSpinor spParent(int) const
Definition: EvtId.hh:27
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
double mass2() const
Definition: EvtVector4R.hh:100
static const EvtGammaMatrix & g1()
double get(int i) const
Definition: EvtVector4R.hh:162
void deleteTree()
Definition: EvtAmp.hh:30
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtSemiLeptonicFF *FormFactors, EvtComplex r00, EvtComplex r01, EvtComplex r10, EvtComplex r11)
static double getMinMass(EvtId i)
Definition: EvtPDL.cpp:342
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
const EvtVector4R & getP4() const
static const EvtGammaMatrix & g2()
double mass() const
void setSpinDensityForwardHelicityBasis(const EvtSpinDensity &rho)
static double getWidth(EvtId i)
Definition: EvtPDL.cpp:352
static const EvtGammaMatrix & g3()
virtual EvtDiracSpinor spParentNeutrino() const
double d3mag() const
EvtVector4C EvtBaryonVARaritaCurrent(const EvtRaritaSchwinger &Bf_vect, const EvtDiracSpinor &Bi, EvtVector4R parent, EvtVector4R daught, const double *ff, int pflag)
EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtSemiLeptonicFF *FormFactors) override
EvtVector4C getVector(int i) const
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void noLifeTime()
Definition: EvtParticle.hh:382
void set(int i, int j, const EvtComplex &rhoij)
EvtDiracSpinor getSpinor(int i) const
EvtVector4C EvtBaryonVACurrent(const EvtDiracSpinor &Bf, const EvtDiracSpinor &Bi, EvtVector4R parent, EvtVector4R daught, const double *ff, int pflag)
EvtComplex cont(const EvtVector4C &v4) const
Definition: EvtVector4C.hh:140
static double getMass(EvtId i)
Definition: EvtPDL.cpp:319
EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void setDim(int n)
void init(EvtId p, int ndaug, EvtId *daug)
Definition: EvtAmp.cpp:68
static const EvtGammaMatrix & g5()
virtual void getbaryonff(EvtId parent, EvtId daught, double t, double m_meson, double *f1v, double *f1a, double *f2v, double *f2a)=0
static double getMaxMass(EvtId i)
Definition: EvtPDL.cpp:337
virtual void getdiracff(EvtId parent, EvtId daught, double q2, double mass, double *f1, double *f2, double *f3, double *g1, double *g2, double *g3)=0
double normalizedProb(const EvtSpinDensity &d)