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
EvtLambdacPHH.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/EvtLambdacPHH.hh"
22
23#include "EvtGenBase/EvtConst.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtPDL.hh"
26#include "EvtGenBase/EvtParticle.hh"
27#include "EvtGenBase/EvtRandom.hh"
28#include "EvtGenBase/EvtReport.hh"
29#include "EvtGenBase/EvtSpinType.hh"
30#include "EvtGenBase/EvtdFunction.hh"
31
32#include <algorithm>
33#include <cmath>
34#include <utility>
35
36EvtLambdacPHH::EvtLambdacPHH() :
37 m_d1( 0 ),
38 m_d2( 1 ),
39 m_d3( 3 ),
40 m_Nplusplus( 0.46 ),
41 m_Nplusminus( 1.0 ),
42 m_Nminusplus( 0.18 ),
43 m_Nminusminus( 0.94 ),
44 m_phiNplusplus( 3.48 ),
45 m_phiNplusminus( 0.00 ),
46 m_phiNminusplus( 0.75 ),
47 m_phiNminusminus( 1.13 ),
48 m_E1( 0.52 ),
49 m_phiE1( -1.01 ),
50 m_E2( 0.20 ),
51 m_phiE2( 2.35 ),
52 m_E3( 0.21 ),
53 m_phiE3( 3.46 ),
54 m_E4( 0.16 ),
55 m_phiE4( 5.29 ),
56 m_F1( 0.17 ),
57 m_phiF1( 4.98 ),
58 m_F2( 0.38 ),
59 m_phiF2( 4.88 ),
60 m_H1( 0.18 ),
61 m_phiH1( 5.93 ),
62 m_H2( 0.20 ),
63 m_phiH2( -0.06 ),
64 m_NRNorm( 1.0 ),
65 m_KstarNorm( 1.0 ),
66 m_DeltaNorm( 1.0 ),
67 m_LambdaNorm( 1.0 ),
68 m_KstarM( 0.890 ),
69 m_KstarW( 0.0498 ),
70 m_KstarR( 3.40 ),
71 m_DeltaM( 1.232 ),
72 m_DeltaW( 0.1120 ),
73 m_DeltaR( 5.22 ),
74 m_LambdaM( 1.520 ),
75 m_LambdaW( 0.0156 ),
76 m_LambdaR( 6.29 ),
77 m_Lambda_cR( 5.07 ),
78 m_zprime(),
79 m_p4_Lambda_c(),
80 m_zpMag( 0.0 ),
81 m_p4_Lambdac_Mag( 0.0 )
82{
83 // Fermilab E791 values from MINUIT fit arXiv:hep-ex/9912003v1
84}
85
86std::string EvtLambdacPHH::getName() const
87{
88 return "LAMBDAC_PHH";
89}
90
91EvtDecayBase* EvtLambdacPHH::clone() const
92{
93 return new EvtLambdacPHH;
94}
95
96bool compareId( const std::pair<EvtId, int>& left,
97 const std::pair<EvtId, int>& right )
98{
99 // Compare id numbers to achieve the ordering K-, pi+ and p
100 bool result( false );
101
102 int leftPDGid = EvtPDL::getStdHep( left.first );
103 int rightPDGid = EvtPDL::getStdHep( right.first );
104
105 if ( leftPDGid < rightPDGid ) {
106 result = true;
107 }
108
109 return result;
110}
111
112void EvtLambdacPHH::init()
113{
114 static const EvtId KM = EvtPDL::getId( "K-" );
115 static const EvtId PIP = EvtPDL::getId( "pi+" );
116 static const EvtId LAMBDAC = EvtPDL::getId( "Lambda_c+" );
117 static const EvtId LAMBDACB = EvtPDL::getId( "anti-Lambda_c-" );
118 static const EvtId PROTON = EvtPDL::getId( "p+" );
119
120 // check that there are 0 or 1 arguments and 3 daughters
121 checkNArg( 0, 1 );
122 checkNDaug( 3 );
123
124 EvtId parnum = getParentId();
125 checkSpinParent( EvtSpinType::DIRAC );
126 checkSpinDaughter( 0, EvtSpinType::DIRAC );
127 checkSpinDaughter( 1, EvtSpinType::SCALAR );
128 checkSpinDaughter( 2, EvtSpinType::SCALAR );
129
130 std::vector<std::pair<EvtId, int>> daughters;
131 if ( parnum == LAMBDAC ) {
132 for ( int i = 0; i < 3; ++i ) {
133 daughters.push_back( std::make_pair( getDaug( i ), i ) );
134 }
135 } else {
136 for ( int i = 0; i < 3; ++i ) {
137 daughters.push_back(
138 std::make_pair( EvtPDL::chargeConj( getDaug( i ) ), i ) );
139 }
140 }
141
142 // Sort daughters, they will end up in the order KM, PIP and PROTON
143 std::sort( daughters.begin(), daughters.end(), compareId );
144
145 if ( parnum == LAMBDAC || parnum == LAMBDACB ) {
146 if ( daughters[0].first == KM && daughters[1].first == PIP &&
147 daughters[2].first == PROTON ) {
148 m_d1 = daughters[0].second;
149 m_d2 = daughters[1].second;
150 m_d3 = daughters[2].second;
151 }
152 }
153
154 // Find resonance dynamics normalisations
155 calcNormalisations();
156
157 // Print out expected fit fractions
158 getFitFractions();
159}
160
161void EvtLambdacPHH::calcNormalisations()
162{
163 // Generate events uniform in the Lambda_c Dalitz plot and find the
164 // normalisation integrals of the Breit-Wigner lineshapes
165
166 // Lambda_c -> K- pi+ p
167 int nDaug( 3 );
168 EvtVector4R p4Daug[3];
169
170 double mDaug[3] = { EvtPDL::getMeanMass( EvtPDL::getId( "K-" ) ),
171 EvtPDL::getMeanMass( EvtPDL::getId( "pi+" ) ),
172 EvtPDL::getMeanMass( EvtPDL::getId( "p+" ) ) };
173
174 double norm[3] = { 0.0, 0.0, 0.0 };
175
176 // sample size
177 int N( 100000 );
178 for ( int i = 0; i < N; i++ ) {
179 double mParent = EvtPDL::getMass( EvtPDL::getId( "Lambda_c+" ) );
180 EvtVector4R p0( mParent, 0.0, 0.0, 0.0 );
181
182 // Generate uniform 4 momenta
183 EvtGenKine::PhaseSpace( nDaug, mDaug, p4Daug, mParent );
184
185 EvtResonance2 LambdacpKpi1( p0, p4Daug[0], p4Daug[1], 1.0, 0.0,
186 m_KstarW, m_KstarM, 1, true, m_KstarR,
187 m_Lambda_cR ); // K*0 -> K- and pi+; L = 1
188 EvtResonance2 LambdacpKpi2( p0, p4Daug[2], p4Daug[1], 1.0, 0.0,
189 m_DeltaW, m_DeltaM, 1, true, m_DeltaR,
190 m_Lambda_cR ); // Delta++ -> p and pi+; L = 1
191 EvtResonance2 LambdacpKpi3(
192 p0, p4Daug[2], p4Daug[0], 1.0, 0.0, m_LambdaW, m_LambdaM, 2, true,
193 m_LambdaR, m_Lambda_cR ); // Lambda(1520) -> K- and p; L = 2
194
195 // Sum amplitude magnitude squared
196 norm[0] += abs2( LambdacpKpi1.resAmpl() );
197 norm[1] += abs2( LambdacpKpi2.resAmpl() );
198 norm[2] += abs2( LambdacpKpi3.resAmpl() );
199 }
200
201 // Set normalisation lineshape multiplication factors
202 double N0( N * 1.0 );
203
204 // Scale NR to get sensible relative fit fractions
205 m_NRNorm = 1.0 / 3.0;
206 // Set this using a decay file parameter if required
207 if ( getNArg() > 1 ) {
208 m_NRNorm = getArg( 1 );
209 }
210
211 if ( norm[0] > 0.0 ) {
212 m_KstarNorm = sqrt( N0 / norm[0] );
213 }
214 if ( norm[1] > 0.0 ) {
215 m_DeltaNorm = sqrt( N0 / norm[1] );
216 }
217 if ( norm[2] > 0.0 ) {
218 m_LambdaNorm = sqrt( N0 / norm[2] );
219 }
220}
221
222void EvtLambdacPHH::getFitFractions()
223{
224 // Generate events uniform in the Lambda_c Dalitz plot and find the
225 // fit fractions for each resonance
226
227 // Lambda_c -> K- pi+ p
228 int nDaug( 3 );
229 EvtVector4R p4Daug[3];
230
231 double mDaug[3] = { EvtPDL::getMeanMass( EvtPDL::getId( "K-" ) ),
232 EvtPDL::getMeanMass( EvtPDL::getId( "pi+" ) ),
233 EvtPDL::getMeanMass( EvtPDL::getId( "p+" ) ) };
234
235 double FitFracTop[4] = { 0.0, 0.0, 0.0, 0.0 };
236 double FitFracDenom = 0.0;
237
238 // sample size
239 int N( 100000 );
240 for ( int i = 0; i < N; i++ ) {
241 double mParent = EvtPDL::getMass( EvtPDL::getId( "Lambda_c+" ) );
242 EvtVector4R p0( mParent, 0.0, 0.0, 0.0 );
243
244 // Generate uniform 4 momenta
245 EvtGenKine::PhaseSpace( nDaug, mDaug, p4Daug, mParent );
246
247 EvtResonance2 LambdacpKpi0( p0, p4Daug[0], p4Daug[1], 1.0, 0.0, 0.0, 0.0,
248 0, true, 0.0, 0.0 ); // Non resonant (NR)
249 EvtResonance2 LambdacpKpi1( p0, p4Daug[0], p4Daug[1], 1.0, 0.0,
250 m_KstarW, m_KstarM, 1, true, m_KstarR,
251 m_Lambda_cR ); // K*0 -> K- and pi+; L = 1
252 EvtResonance2 LambdacpKpi2( p0, p4Daug[2], p4Daug[1], 1.0, 0.0,
253 m_DeltaW, m_DeltaM, 1, true, m_DeltaR,
254 m_Lambda_cR ); // Delta++ -> p and pi+; L = 1
255 EvtResonance2 LambdacpKpi3(
256 p0, p4Daug[2], p4Daug[0], 1.0, 0.0, m_LambdaW, m_LambdaM, 2, true,
257 m_LambdaR, m_Lambda_cR ); // Lambda(1520) -> K- and p; L = 2
258
259 std::vector<EvtComplex> ampNonRes =
260 calcResAmpTerms( LcResLabel::NonReson, LambdacpKpi0, m_NRNorm );
261 std::vector<EvtComplex> ampKstar =
262 calcResAmpTerms( LcResLabel::Kstar, LambdacpKpi1, m_KstarNorm );
263 std::vector<EvtComplex> ampDelta =
264 calcResAmpTerms( LcResLabel::Delta, LambdacpKpi2, m_DeltaNorm );
265 std::vector<EvtComplex> ampLambda =
266 calcResAmpTerms( LcResLabel::Lambda, LambdacpKpi3, m_LambdaNorm );
267
268 // Combine resonance amplitudes for a given spin configuration
269 EvtComplex amp00 = ampNonRes[0] + ampKstar[0] + ampDelta[0] +
270 ampLambda[0];
271 EvtComplex amp01 = ampNonRes[1] + ampKstar[1] + ampDelta[1] +
272 ampLambda[1];
273 EvtComplex amp10 = ampNonRes[2] + ampKstar[2] + ampDelta[2] +
274 ampLambda[2];
275 EvtComplex amp11 = ampNonRes[3] + ampKstar[3] + ampDelta[3] +
276 ampLambda[3];
277
278 // Fit fraction numerator terms
279 FitFracTop[0] += abs2( ampNonRes[0] ) + abs2( ampNonRes[1] ) +
280 abs2( ampNonRes[2] ) + abs2( ampNonRes[3] );
281 FitFracTop[1] += abs2( ampKstar[0] ) + abs2( ampKstar[1] ) +
282 abs2( ampKstar[2] ) + abs2( ampKstar[3] );
283 FitFracTop[2] += abs2( ampDelta[0] ) + abs2( ampDelta[1] ) +
284 abs2( ampDelta[2] ) + abs2( ampDelta[3] );
285 FitFracTop[3] += abs2( ampLambda[0] ) + abs2( ampLambda[1] ) +
286 abs2( ampLambda[2] ) + abs2( ampLambda[3] );
287
288 // Fit fraction common denominator
289 FitFracDenom += abs2( amp00 ) + abs2( amp01 ) + abs2( amp10 ) +
290 abs2( amp11 );
291 }
292
293 EvtGenReport( EVTGEN_INFO, "EvtLambdacPHH" )
294 << "FitFracs: NR = " << FitFracTop[0] / FitFracDenom
295 << ", K* = " << FitFracTop[1] / FitFracDenom
296 << ", Del = " << FitFracTop[2] / FitFracDenom
297 << ", Lam = " << FitFracTop[3] / FitFracDenom << std::endl;
298}
299
300void EvtLambdacPHH::initProbMax()
301{
302 // Default value
303 setProbMax( 10.0 );
304
305 // Set probability using decay file parameter
306 if ( getNArg() > 0 ) {
307 setProbMax( getArg( 0 ) );
308 }
309}
310
311void EvtLambdacPHH::decay( EvtParticle* p )
312{
313 // Daughter order: 1 = K-, 2 = pi+, 3 = p
314 p->initializePhaseSpace( getNDaug(), getDaugs() );
315
316 // 4-momenta in the rest frame of the Lambda_c
317 EvtVector4R p4_p( p->mass(), 0.0, 0.0, 0.0 );
318 EvtVector4R moms1 = p->getDaug( m_d1 )->getP4();
319 EvtVector4R moms2 = p->getDaug( m_d2 )->getP4();
320 EvtVector4R moms3 = p->getDaug( m_d3 )->getP4();
321
322 // Lambda_c decay mode resonances. Spin L values from strong decay parity conservation:
323 // parity(resonance) = parity(daug1)*parity(daug2)*(-1)^L
324 EvtResonance2 LambdacpKpi0( p4_p, moms1, moms2, 1.0, 0.0, 0.0, 0.0, 0, true,
325 0.0, 0.0 ); // Non-resonant L = 0
326 EvtResonance2 LambdacpKpi1( p4_p, moms1, moms2, 1.0, 0.0, m_KstarW,
327 m_KstarM, 1, true, m_KstarR,
328 m_Lambda_cR ); // K*0 -> K- and pi+; L = 1
329 EvtResonance2 LambdacpKpi2( p4_p, moms3, moms2, 1.0, 0.0, m_DeltaW,
330 m_DeltaM, 1, true, m_DeltaR,
331 m_Lambda_cR ); // Delta++ -> p and pi+; L = 1
332 EvtResonance2 LambdacpKpi3( p4_p, moms3, moms1, 1.0, 0.0, m_LambdaW,
333 m_LambdaM, 2, true, m_LambdaR,
334 m_Lambda_cR ); // Lambda(1520) -> K- and p; L = 2
335
336 // Define the "beam" direction, used in Fig 1 of hep-ex/9912003v1
337 EvtVector4R beam( 0.0, 0.0, 0.0, 1.0 );
338 EvtParticle* parent = p->getParent();
339 if ( parent ) {
340 // If non prompt, the beam is along the direction of the mother
341 EvtVector4R p4_Lambda_c_mother = parent->getP4Lab();
342 p4_Lambda_c_mother.applyBoostTo( p->getP4Lab() );
343 beam = p4_Lambda_c_mother;
344 }
345
346 m_p4_Lambda_c = p->getP4Lab();
347 m_p4_Lambdac_Mag = m_p4_Lambda_c.d3mag();
348
349 // Define the unit vector denoting the "z" axis in Fig 1
350 m_zprime = -1.0 * m_p4_Lambda_c.cross( beam );
351 m_zprime.applyBoostTo( m_p4_Lambda_c, true ); // From lab frame to Lambda_c
352
353 m_zpMag = m_zprime.d3mag();
354 // Check if zprime magnitude is non-zero
355 if ( m_zpMag > 0.0 ) {
356 // Normalise
357 m_zprime /= m_zpMag;
358 } else {
359 // Set as the z direction
360 m_zprime.set( 0.0, 0.0, 0.0, 1.0 );
361 }
362 // Update normalised |z'|
363 m_zpMag = 1.0;
364
365 // Get the amplitudes: non-resonant, K*, Delta and Lambda
366 std::vector<EvtComplex> ampNonRes =
367 calcResAmpTerms( LcResLabel::NonReson, LambdacpKpi0, m_NRNorm );
368 std::vector<EvtComplex> ampKstar =
369 calcResAmpTerms( LcResLabel::Kstar, LambdacpKpi1, m_KstarNorm );
370 std::vector<EvtComplex> ampDelta =
371 calcResAmpTerms( LcResLabel::Delta, LambdacpKpi2, m_DeltaNorm );
372 std::vector<EvtComplex> ampLambda =
373 calcResAmpTerms( LcResLabel::Lambda, LambdacpKpi3, m_LambdaNorm );
374
375 // Combine resonance amplitudes for a given spin configuration
376 EvtComplex amp00 = ampNonRes[0] + ampKstar[0] + ampDelta[0] + ampLambda[0];
377 EvtComplex amp01 = ampNonRes[1] + ampKstar[1] + ampDelta[1] + ampLambda[1];
378 EvtComplex amp10 = ampNonRes[2] + ampKstar[2] + ampDelta[2] + ampLambda[2];
379 EvtComplex amp11 = ampNonRes[3] + ampKstar[3] + ampDelta[3] + ampLambda[3];
380
381 // Set the amplitude components
382 vertex( 0, 0, amp00 );
383 vertex( 0, 1, amp01 );
384 vertex( 1, 0, amp10 );
385 vertex( 1, 1, amp11 );
386}
387
388std::vector<EvtComplex> EvtLambdacPHH::calcResAmpTerms(
389 EvtLambdacPHH::LcResLabel resIndex, const EvtResonance2& res, double norm ) const
390{
391 // Initialise the resonance and daughter theta and phi angles
392 double thetaRes( 0.0 ), phiRes( 0.0 ), phiPrimeDaug( 0.0 ),
393 thetaPrimeDaug( 0.0 );
394 // Initialise beta rotation angle
395 double beta_res( 0.0 );
396
397 EvtVector4R res_atproton( 0.0, 0.0, 0.0, 0.0 ),
398 Lc_atproton( 0.0, 0.0, 0.0, 0.0 );
399
400 // Initialise Amplitude terms
401 EvtComplex term1( 0.0 ), term2( 0.0 ), term3( 0.0 ), term4( 0.0 );
402 // Normalised dynamical amplitude
403 EvtComplex resAmp( norm, 0.0 );
404
405 // Angles are not needed for the non-resonant amplitude
406 if ( resIndex != LcResLabel::NonReson ) {
407 resAmp = res.resAmpl() * norm;
408 // Resonance and daughter 4 momenta
409 EvtVector4R p4d1 = res.p4_d1();
410 EvtVector4R p4d2 = res.p4_d2();
411 EvtVector4R p4Res = p4d1 + p4d2;
412 EvtVector4R p4_d3 = res.p4_p() - p4Res;
413
414 double p4ResMag = p4Res.d3mag();
415
416 // 4-momenta for theta' and phi' angles
417 EvtVector4R yRes = -1.0 * p4_d3.cross( m_zprime );
418
419 EvtVector4R res_d1 = p4d1;
420 res_d1.applyBoostTo( p4Res, true );
421 double res_d1_Mag = res_d1.d3mag();
422
423 EvtVector4R res_d3 = -1.0 * p4_d3;
424 double res_d3_Mag = res_d3.d3mag();
425
426 thetaPrimeDaug = getACos( res_d1.dot( res_d3 ), res_d1_Mag * res_d3_Mag );
427
428 res_atproton = p4Res;
429 res_atproton.applyBoostTo( p4d1, true );
430 double res_atproton_mag = res_atproton.d3mag();
431
432 Lc_atproton = res.p4_p();
433 Lc_atproton.applyBoostTo( p4d1, true );
434 double Lc_atproton_mag = Lc_atproton.d3mag();
435
436 // Check that the momentum of the Lambda_c is not zero, as well as a valid zprime vector
437 if ( m_p4_Lambdac_Mag > 0.0 && m_zpMag > 0.0 ) {
438 thetaRes = getACos( -1.0 * p4Res.dot( m_zprime ), p4ResMag );
439 phiRes = getASin( -1.0 * p4Res.dot( m_p4_Lambda_c ),
440 sin( thetaRes ) * m_p4_Lambdac_Mag * p4ResMag );
441 phiPrimeDaug = getASin( res_d1.dot( yRes ), sin( thetaPrimeDaug ) *
442 res_d1_Mag *
443 yRes.d3mag() );
444
445 } else {
446 // Use randomised angles with flat probability distributions
447 thetaRes = EvtRandom::Flat( 0.0, EvtConst::pi );
448 phiRes = EvtRandom::Flat( 0.0, EvtConst::twoPi );
449 phiPrimeDaug = EvtRandom::Flat( 0.0, EvtConst::twoPi );
450 }
451
452 if ( res_atproton_mag > 0.0 && Lc_atproton_mag > 0.0 ) {
453 // Extra rotation to go to the proton helicity frame for the two resonances Delta++ and Lambda.
454 // No rotation is needed for K*. Use the momenta boosted to the proton restframe
455
456 beta_res = getACos( res_atproton.dot( Lc_atproton ),
457 res_atproton_mag * Lc_atproton_mag );
458
459 } else {
460 beta_res = EvtRandom::Flat( 0.0, EvtConst::pi );
461 }
462 }
463
464 // Find the spin-dependent amplitudes
465 if ( resIndex == LcResLabel::NonReson || resIndex == LcResLabel::Kstar ) {
466 term1 = resAmp * DecayAmp3( resIndex, 1, 1, thetaRes, phiRes,
467 thetaPrimeDaug, phiPrimeDaug );
468 term2 = resAmp * DecayAmp3( resIndex, 1, -1, thetaRes, phiRes,
469 thetaPrimeDaug, phiPrimeDaug );
470 term3 = resAmp * DecayAmp3( resIndex, -1, 1, thetaRes, phiRes,
471 thetaPrimeDaug, phiPrimeDaug );
472 term4 = resAmp * DecayAmp3( resIndex, -1, -1, thetaRes, phiRes,
473 thetaPrimeDaug, phiPrimeDaug );
474
475 } else {
476 double rotate_00 = EvtdFunction::d( 1, 1, 1, beta_res );
477 double rotate_10 = EvtdFunction::d( 1, -1, 1, beta_res );
478 double rotate_11 = EvtdFunction::d( 1, -1, -1, beta_res );
479 double rotate_01 = EvtdFunction::d( 1, 1, -1, beta_res );
480
481 // Delta and Lambda need to be rotated before summing over the proton helicity axis
482 EvtComplex termA = resAmp * DecayAmp3( resIndex, 1, 1, thetaRes, phiRes,
483 thetaPrimeDaug, phiPrimeDaug );
484 EvtComplex termB = resAmp * DecayAmp3( resIndex, 1, -1, thetaRes, phiRes,
485 thetaPrimeDaug, phiPrimeDaug );
486 EvtComplex termC = resAmp * DecayAmp3( resIndex, -1, 1, thetaRes, phiRes,
487 thetaPrimeDaug, phiPrimeDaug );
488 EvtComplex termD = resAmp * DecayAmp3( resIndex, -1, -1, thetaRes, phiRes,
489 thetaPrimeDaug, phiPrimeDaug );
490
491 term1 = rotate_00 * termA + rotate_10 * termB;
492 term2 = rotate_01 * termA + rotate_11 * termB;
493 term3 = rotate_00 * termC + rotate_10 * termD;
494 term4 = rotate_01 * termC + rotate_11 * termD;
495 }
496
497 // Return the spin amplitudes as a vector
498 std::vector<EvtComplex> ampVect;
499 ampVect.push_back( term1 );
500 ampVect.push_back( term2 );
501 ampVect.push_back( term3 );
502 ampVect.push_back( term4 );
503
504 return ampVect;
505}
506
507EvtComplex EvtLambdacPHH::DecayAmp3( EvtLambdacPHH::LcResLabel resonance, int m,
508 int mprime, double theta_res, double phi_res,
509 double theta_prime_daughter_res,
510 double phi_prime_daughter_res ) const
511{
512 // Find the amplitudes given in Tables 3 to 6 in the paper.
513 // Wigner d-functions use 2*spin, e.g. d(1/2, 1/2, 1/2) -> d(1, 1, 1)
514 EvtComplex term1( 0.0, 0.0 ), term2( 0.0, 0.0 );
515
516 if ( resonance == LcResLabel::NonReson ) {
517 // Non-resonant: table 6
518 if ( m == 1 && mprime == 1 ) {
519 term1 = m_Nplusplus *
520 EvtComplex( cos( m_phiNplusplus ), sin( m_phiNplusplus ) );
521
522 } else if ( m == 1 && mprime == -1 ) {
523 term1 = m_Nplusminus * EvtComplex( cos( m_phiNplusminus ),
524 sin( m_phiNplusminus ) );
525
526 } else if ( m == -1 && mprime == 1 ) {
527 term1 = m_Nminusplus * EvtComplex( cos( m_phiNminusplus ),
528 sin( m_phiNminusplus ) );
529
530 } else if ( m == -1 && mprime == -1 ) {
531 term1 = m_Nminusminus * EvtComplex( cos( m_phiNminusminus ),
532 sin( m_phiNminusminus ) );
533 }
534
535 } else if ( resonance == LcResLabel::Kstar ) {
536 // K*0(1-) resonance: table 3
537 if ( m == 1 && mprime == 1 ) {
538 term1 = fampl3( m_E1, m_phiE1, 1, 1, 1, theta_res, 2, 2, 0,
539 theta_prime_daughter_res, phi_prime_daughter_res );
540 term2 = fampl3( m_E2, m_phiE2, 1, 1, -1, theta_res, 2, 0, 0,
541 theta_prime_daughter_res, phi_res );
542
543 } else if ( m == 1 && mprime == -1 ) {
544 term1 = fampl3( m_E3, m_phiE3, 1, 1, 1, theta_res, 2, 0, 0,
545 theta_prime_daughter_res, 0.0 );
546 term2 = fampl3( m_E4, m_phiE4, 1, 1, -1, theta_res, 2, -2, 0,
547 theta_prime_daughter_res,
548 phi_res - phi_prime_daughter_res );
549
550 } else if ( m == -1 && mprime == 1 ) {
551 term1 = fampl3( m_E1, m_phiE1, 1, -1, 1, theta_res, 2, 2, 0,
552 theta_prime_daughter_res,
553 -( phi_res - phi_prime_daughter_res ) );
554 term2 = fampl3( m_E2, m_phiE2, 1, -1, -1, theta_res, 2, 0, 0,
555 theta_prime_daughter_res, 0.0 );
556
557 } else if ( m == -1 && mprime == -1 ) {
558 term1 = fampl3( m_E3, m_phiE3, 1, -1, 1, theta_res, 2, 0, 0,
559 theta_prime_daughter_res, -phi_res );
560 term2 = fampl3( m_E4, m_phiE4, 1, -1, -1, theta_res, 2, -2, 0,
561 theta_prime_daughter_res, -phi_prime_daughter_res );
562 }
563
564 } else if ( resonance == LcResLabel::Delta ) {
565 // Delta++(3/2+) resonance: table 4
566 if ( m == 1 && mprime == 1 ) {
567 term1 = fampl3( m_F1, m_phiF1, 1, 1, 1, theta_res, 3, 1, 1,
568 theta_prime_daughter_res, 0.0 );
569 term2 = fampl3( m_F2, m_phiF2, 1, 1, -1, theta_res, 3, -1, 1,
570 theta_prime_daughter_res,
571 phi_res - phi_prime_daughter_res );
572
573 } else if ( m == 1 && mprime == -1 ) {
574 term1 = fampl3( m_F1, m_phiF1, 1, 1, 1, theta_res, 3, 1, -1,
575 theta_prime_daughter_res, phi_prime_daughter_res );
576 term2 = fampl3( m_F2, m_phiF2, 1, 1, -1, theta_res, 3, -1, -1,
577 theta_prime_daughter_res, phi_res );
578
579 } else if ( m == -1 && mprime == 1 ) {
580 term1 = fampl3( m_F1, m_phiF1, 1, -1, 1, theta_res, 3, 1, 1,
581 theta_prime_daughter_res, -phi_res );
582 term2 = fampl3( m_F2, m_phiF2, 1, -1, -1, theta_res, 3, -1, 1,
583 theta_prime_daughter_res, -phi_prime_daughter_res );
584
585 } else if ( m == -1 && mprime == -1 ) {
586 term1 = fampl3( m_F1, m_phiF1, 1, -1, 1, theta_res, 3, 1, -1,
587 theta_prime_daughter_res,
588 -( phi_res - phi_prime_daughter_res ) );
589 term2 = fampl3( m_F2, m_phiF2, 1, -1, -1, theta_res, 3, -1, -1,
590 theta_prime_daughter_res, 0.0 );
591 }
592
593 } else if ( resonance == LcResLabel::Lambda ) {
594 // Lambda(1520)(3/2-) resonance: table 5
595 if ( m == 1 && mprime == 1 ) {
596 term1 = fampl3( m_H1, m_phiH1, 1, 1, 1, theta_res, 3, 1, 1,
597 theta_prime_daughter_res, 0.0 );
598 term2 = fampl3( m_H2, m_phiH2, 1, 1, -1, theta_res, 3, -1, 1,
599 theta_prime_daughter_res,
600 phi_res - phi_prime_daughter_res );
601
602 } else if ( m == 1 && mprime == -1 ) {
603 term1 = -1.0 * fampl3( m_H1, m_phiH1, 1, 1, 1, theta_res, 3, 1, -1,
604 theta_prime_daughter_res,
605 phi_prime_daughter_res );
606 term2 = -1.0 * fampl3( m_H2, m_phiH2, 1, 1, -1, theta_res, 3, -1,
607 -1, theta_prime_daughter_res, phi_res );
608
609 } else if ( m == -1 && mprime == 1 ) {
610 term1 = fampl3( m_H1, m_phiH1, 1, -1, 1, theta_res, 3, 1, 1,
611 theta_prime_daughter_res, -phi_res );
612 term2 = fampl3( m_H2, m_phiH2, 1, -1, -1, theta_res, 3, -1, 1,
613 theta_prime_daughter_res, -phi_prime_daughter_res );
614
615 } else if ( m == -1 && mprime == -1 ) {
616 term1 = -1.0 * fampl3( m_H1, m_phiH1, 1, -1, 1, theta_res, 3, 1, -1,
617 theta_prime_daughter_res,
618 -( phi_res - phi_prime_daughter_res ) );
619 term2 = -1.0 * fampl3( m_H2, m_phiH2, 1, -1, -1, theta_res, 3, -1,
620 -1, theta_prime_daughter_res, 0.0 );
621 }
622 }
623
624 EvtComplex Amplitude = term1 + term2;
625 return Amplitude;
626}
627
628EvtComplex EvtLambdacPHH::fampl3( double amplitude_res, double phi_res,
629 int spinMother, int m_spinMother,
630 int m_prime_spinMother, double theta_res,
631 float spin_res, float m_spin_res,
632 float m_prime_spin_res,
633 double theta_daughter_res,
634 double phi_prime_daughter_res ) const
635{
636 double dTerm1 = EvtdFunction::d( spinMother, m_spinMother,
637 m_prime_spinMother, theta_res );
638 double dTerm2 = EvtdFunction::d( spin_res, m_spin_res, m_prime_spin_res,
639 theta_daughter_res );
640
641 EvtComplex amp_phase1 = EvtComplex( cos( phi_res ), sin( phi_res ) );
642 EvtComplex amp_phase2 = EvtComplex( cos( phi_prime_daughter_res ),
643 sin( phi_prime_daughter_res ) );
644
645 EvtComplex partial_amp = amplitude_res * amp_phase1 * dTerm1 * amp_phase2 *
646 dTerm2;
647
648 return partial_amp;
649}
650
651double EvtLambdacPHH::getACos( double num, double denom ) const
652{
653 // Find inverse cosine, checking ratio is within +- 1
654 double angle( 0.0 ), ratio( 0.0 );
655 if ( fabs( denom ) > 0.0 ) {
656 ratio = num / denom;
657 }
658
659 if ( fabs( ratio ) <= 1.0 ) {
660 angle = acos( ratio );
661 }
662
663 return angle;
664}
665
666double EvtLambdacPHH::getASin( double num, double denom ) const
667{
668 // Find inverse sine, checking ratio is within +- 1
669 double angle( 0.0 ), ratio( 0.0 );
670 if ( fabs( denom ) > 0.0 ) {
671 ratio = num / denom;
672 }
673
674 if ( fabs( ratio ) <= 1.0 ) {
675 angle = asin( ratio );
676 }
677
678 return angle;
679}
abs2
double abs2(const EvtComplex &c)
Definition EvtComplex.hh:211
compareId
bool compareId(const std::pair< EvtId, int > &left, const std::pair< EvtId, int > &right)
Definition EvtLambdacPHH.cpp:96
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
EVTGEN_INFO
@ EVTGEN_INFO
Definition EvtReport.hh:52
EvtConst::pi
static const double pi
Definition EvtConst.hh:26
EvtConst::twoPi
static const double twoPi
Definition EvtConst.hh:27
EvtDecayAmp::vertex
void vertex(const EvtComplex &amp)
Definition EvtDecayAmp.hh:37
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cpp:547
EvtDecayBase::EvtDecayBase
EvtDecayBase()=default
EvtDecayBase::getNDaug
int getNDaug() const
Definition EvtDecayBase.hh:64
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cpp:534
EvtDecayBase::getNArg
int getNArg() const
Definition EvtDecayBase.hh:67
EvtDecayBase::getArg
double getArg(unsigned int j)
Definition EvtDecayBase.cpp:578
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cpp:295
EvtDecayBase::getParentId
EvtId getParentId() const
Definition EvtDecayBase.hh:60
EvtDecayBase::getDaug
EvtId getDaug(int i) const
Definition EvtDecayBase.hh:66
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cpp:516
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cpp:492
EvtDecayBase::getDaugs
const EvtId * getDaugs() const
Definition EvtDecayBase.hh:65
EvtGenKine::PhaseSpace
static double PhaseSpace(int ndaug, double mass[30], EvtVector4R p4[30], double mp)
Definition EvtGenKine.cpp:46
EvtId
Definition EvtId.hh:27
EvtLambdacPHH::m_phiNminusplus
double m_phiNminusplus
Definition EvtLambdacPHH.hh:89
EvtLambdacPHH::calcResAmpTerms
std::vector< EvtComplex > calcResAmpTerms(EvtLambdacPHH::LcResLabel resIndex, const EvtResonance2 &res, double norm) const
Definition EvtLambdacPHH.cpp:388
EvtLambdacPHH::m_p4_Lambda_c
EvtVector4R m_p4_Lambda_c
Definition EvtLambdacPHH.hh:99
EvtLambdacPHH::decay
void decay(EvtParticle *p) override
Definition EvtLambdacPHH.cpp:311
EvtLambdacPHH::m_phiNplusplus
double m_phiNplusplus
Definition EvtLambdacPHH.hh:89
EvtLambdacPHH::m_phiNplusminus
double m_phiNplusminus
Definition EvtLambdacPHH.hh:89
EvtLambdacPHH::clone
EvtDecayBase * clone() const override
Definition EvtLambdacPHH.cpp:91
EvtLambdacPHH::getACos
double getACos(double num, double denom) const
Definition EvtLambdacPHH.cpp:651
EvtLambdacPHH::initProbMax
void initProbMax() override
Definition EvtLambdacPHH.cpp:300
EvtLambdacPHH::getName
std::string getName() const override
Definition EvtLambdacPHH.cpp:86
EvtLambdacPHH::m_phiNminusminus
double m_phiNminusminus
Definition EvtLambdacPHH.hh:89
EvtLambdacPHH::m_p4_Lambdac_Mag
double m_p4_Lambdac_Mag
Definition EvtLambdacPHH.hh:100
EvtLambdacPHH::getASin
double getASin(double num, double denom) const
Definition EvtLambdacPHH.cpp:666
EvtLambdacPHH::calcNormalisations
void calcNormalisations()
Definition EvtLambdacPHH.cpp:161
EvtLambdacPHH::fampl3
EvtComplex fampl3(double amplitude_res, double phi_res, int spinMother, int m_spinMother, int m_prime_spinMother, double theta_res, float spin_res, float m_spin_res, float m_prime_spin_res, double theta_daughter_res, double phi_prime_daughter_res) const
Definition EvtLambdacPHH.cpp:628
EvtLambdacPHH::init
void init() override
Definition EvtLambdacPHH.cpp:112
EvtLambdacPHH::m_zprime
EvtVector4R m_zprime
Definition EvtLambdacPHH.hh:99
EvtLambdacPHH::DecayAmp3
EvtComplex DecayAmp3(EvtLambdacPHH::LcResLabel resonance, int m, int mprime, double theta_res, double phi_res, double theta_prime_daughter_res, double phi_prime_daughter_res) const
Definition EvtLambdacPHH.cpp:507
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.cpp:356
EvtPDL::getMass
static double getMass(EvtId i)
Definition EvtPDL.cpp:311
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.cpp:306
EvtPDL::chargeConj
static EvtId chargeConj(EvtId id)
Definition EvtPDL.cpp:201
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cpp:283
EvtParticle::initializePhaseSpace
double initializePhaseSpace(size_t numdaughter, const EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cpp:1100
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::getP4Lab
EvtVector4R getP4Lab() const
Definition EvtParticle.cpp:756
EvtParticle::getParent
EvtParticle * getParent() const
Definition EvtParticle.cpp:90
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cpp:95
EvtResonance2::p4_d1
const EvtVector4R & p4_d1() const
Definition EvtResonance2.hh:43
EvtResonance2::p4_p
const EvtVector4R & p4_p() const
Definition EvtResonance2.hh:42
EvtResonance2::resAmpl
EvtComplex resAmpl() const
Definition EvtResonance2.cpp:66
EvtResonance2::p4_d2
const EvtVector4R & p4_d2() const
Definition EvtResonance2.hh:44
EvtVector4R::dot
double dot(const EvtVector4R &v2) const
Definition EvtVector4R.cpp:205
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cpp:192
EvtVector4R::cross
EvtVector4R cross(const EvtVector4R &v2) const
Definition EvtVector4R.cpp:177
EvtVector4R::applyBoostTo
void applyBoostTo(const EvtVector4R &p4, bool inverse=false)
Definition EvtVector4R.cpp:111
EvtdFunction::d
static double d(int j, int m1, int m2, double theta)
Definition EvtdFunction.cpp:30
Generated by doxygen 1.16.1