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