diff --git a/EvtGenModels/EvtVtoSll.hh b/EvtGenModels/EvtVtoSll.hh --- a/EvtGenModels/EvtVtoSll.hh +++ b/EvtGenModels/EvtVtoSll.hh @@ -32,7 +32,14 @@ void initProbMax() override; void init() override; - void decay( EvtParticle* p ) override; + void decay( EvtParticle* parent ) override; + + protected: + void calcAmp( const EvtParticle& parent, EvtAmp& amp ) const; + + private: + static constexpr double m_poleSize{ 0.00005 }; + bool m_electronMode{ false }; }; #endif diff --git a/History.md b/History.md --- a/History.md +++ b/History.md @@ -11,6 +11,10 @@ === ## R02-0X-00 +22 Aug 2023 Fernando Abudinen +* D97: Bugfix probmax issue and introduced pole compensation for VTOSLL model + Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation. + 22 Aug 2023 Tom Latham * CMake updates - Update default C++ standard to 17 diff --git a/src/EvtGenModels/EvtVtoSll.cpp b/src/EvtGenModels/EvtVtoSll.cpp --- a/src/EvtGenModels/EvtVtoSll.cpp +++ b/src/EvtGenModels/EvtVtoSll.cpp @@ -22,12 +22,14 @@ #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtGenKine.hh" +#include "EvtGenBase/EvtIdSet.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtVector4C.hh" +#include "EvtGenBase/EvtVectorParticle.hh" #include #include @@ -54,58 +56,185 @@ checkSpinDaughter( 0, EvtSpinType::SCALAR ); checkSpinDaughter( 1, EvtSpinType::DIRAC ); checkSpinDaughter( 2, EvtSpinType::DIRAC ); + + // Work out whether we have electron mode + const EvtIdSet leptons{ "e-", "e+" }; + if ( leptons.contains( getDaug( 1 ) ) || leptons.contains( getDaug( 2 ) ) ) { + m_electronMode = true; + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " EvtVtoSll has dielectron final state" << std::endl; + } } void EvtVtoSll::initProbMax() { - setProbMax( 10.0 ); + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " EvtVtoSll is finding maximum probability ... " << std::endl; + + double theProbMax = 0; + double theProbMax_q2 = 0; + double theProbMax_theta = 0; + + EvtVectorParticle parent{}; + parent.noLifeTime(); + parent.init( getParentId(), + EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) ); + parent.setDiagonalSpinDensity(); + + EvtAmp amp; + EvtId daughters[3] = { getDaug( 0 ), getDaug( 1 ), getDaug( 2 ) }; + amp.init( getParentId(), 3, daughters ); + parent.makeDaughters( 3, daughters ); + EvtParticle* scalar = parent.getDaug( 0 ); + EvtParticle* lep1 = parent.getDaug( 1 ); + EvtParticle* lep2 = parent.getDaug( 2 ); + scalar->noLifeTime(); + lep1->noLifeTime(); + lep2->noLifeTime(); + + EvtSpinDensity rho; + rho.setDiag( parent.getSpinStates() ); + + const double M0 = EvtPDL::getMass( getParentId() ); + const double mL = EvtPDL::getMass( getDaug( 0 ) ); + const double m1 = EvtPDL::getMass( getDaug( 1 ) ); + const double m2 = EvtPDL::getMass( getDaug( 2 ) ); + + const double m12Sum = m1 + m2; + const double m12Del = m1 - m2; + + const double q2min = ( m1 + m2 ) * ( m1 + m2 ); + const double q2max = ( M0 - mL ) * ( M0 - mL ); + const double mSqSum = M0 * M0 + mL * mL; + + EvtVector4R p4scalar, p4lep1, p4lep2, boost; + + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " EvtVtoSll is probing whole phase space ..." << std::endl; + + double prob = 0; + const int nsteps = 5000; + for ( int i = 0; i <= nsteps; i++ ) { + const double q2 = q2min + i * ( q2max - q2min ) / nsteps; + const double eScalar = ( mSqSum - q2 ) / ( 2 * M0 ); + const double pstar = i == 0 ? 0 + : sqrt( q2 - m12Sum * m12Sum ) * + sqrt( q2 - m12Del * m12Del ) / + ( 2 * sqrt( q2 ) ); + + boost.set( M0 - eScalar, 0, 0, +sqrt( eScalar * eScalar - mL * mL ) ); + if ( i != nsteps ) { + p4scalar.set( eScalar, 0, 0, -sqrt( eScalar * eScalar - mL * mL ) ); + } else { + p4scalar.set( mL, 0, 0, 0 ); + } + + const double pstarSq = pstar * pstar; + const double E1star = sqrt( pstarSq + m1 * m1 ); + const double E2star = sqrt( pstarSq + m2 * m2 ); + + for ( int j = 0; j <= 45; j++ ) { + const double theta = j * EvtConst::pi / 45; + + const double pstarT = pstar * sin( theta ); + const double pstarZ = pstar * cos( theta ); + + p4lep1.set( E1star, 0, pstarT, pstarZ ); + p4lep2.set( E2star, 0, -pstarT, -pstarZ ); + + if ( i != nsteps ) // At maximal q2 we are already in correct frame as scalar child and W/Zvirtual are at rest + { + p4lep1 = boostTo( p4lep1, boost ); + p4lep2 = boostTo( p4lep2, boost ); + } + scalar->init( getDaug( 0 ), p4scalar ); + lep1->init( getDaug( 1 ), p4lep1 ); + lep2->init( getDaug( 2 ), p4lep2 ); + calcAmp( parent, amp ); + prob = rho.normalizedProb( amp.getSpinDensity() ); + // In case of electron mode add pole + if ( m_electronMode ) { + prob /= 1.0 + m_poleSize / ( q2 * q2 ); + } + + if ( prob > theProbMax ) { + theProbMax = prob; + theProbMax_q2 = q2; + theProbMax_theta = theta; + } + } + } + + theProbMax *= 1.01; + + setProbMax( theProbMax ); + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " EvtVtoSll set up maximum probability to " << theProbMax + << " found at q2 = " << theProbMax_q2 << " (" + << nsteps * ( theProbMax_q2 - q2min ) / ( q2max - q2min ) + << " %) and theta = " << theProbMax_theta * 180 / EvtConst::pi + << std::endl; } -void EvtVtoSll::decay( EvtParticle* p ) +void EvtVtoSll::decay( EvtParticle* parent ) { - p->initializePhaseSpace( getNDaug(), getDaugs() ); - - EvtParticle *l1, *l2; - l1 = p->getDaug( 1 ); - l2 = p->getDaug( 2 ); - - EvtVector4C l11, l12, l21, l22; - l11 = EvtLeptonVCurrent( l1->spParent( 0 ), l2->spParent( 0 ) ); - l12 = EvtLeptonVCurrent( l1->spParent( 0 ), l2->spParent( 1 ) ); - l21 = EvtLeptonVCurrent( l1->spParent( 1 ), l2->spParent( 0 ) ); - l22 = EvtLeptonVCurrent( l1->spParent( 1 ), l2->spParent( 1 ) ); - - EvtVector4C eps0 = p->eps( 0 ); - EvtVector4C eps1 = p->eps( 1 ); - EvtVector4C eps2 = p->eps( 2 ); - - EvtVector4R P = p->getP4Restframe(); - EvtVector4R k = l1->getP4() + l2->getP4(); - double k2 = k * k; - - EvtTensor4C T( dual( EvtGenFunctions::directProd( P, ( 1.0 / k2 ) * k ) ) ); + // Phase space initialization depends on what leptons are + if ( m_electronMode ) { + setWeight( parent->initializePhaseSpace( getNDaug(), getDaugs(), false, + m_poleSize, 1, 2 ) ); + } else { + parent->initializePhaseSpace( getNDaug(), getDaugs() ); + } + + calcAmp( *parent, _amp2 ); +} - double M2 = p->mass(); +void EvtVtoSll::calcAmp( const EvtParticle& parent, EvtAmp& amp ) const +{ + const EvtParticle& l1 = *parent.getDaug( 1 ); + const EvtParticle& l2 = *parent.getDaug( 2 ); + + const EvtVector4C l11 = EvtLeptonVCurrent( l1.spParent( 0 ), + l2.spParent( 0 ) ); + const EvtVector4C l12 = EvtLeptonVCurrent( l1.spParent( 0 ), + l2.spParent( 1 ) ); + const EvtVector4C l21 = EvtLeptonVCurrent( l1.spParent( 1 ), + l2.spParent( 0 ) ); + const EvtVector4C l22 = EvtLeptonVCurrent( l1.spParent( 1 ), + l2.spParent( 1 ) ); + + const EvtVector4C eps0 = parent.eps( 0 ); + const EvtVector4C eps1 = parent.eps( 1 ); + const EvtVector4C eps2 = parent.eps( 2 ); + + const EvtVector4R P = parent.getP4Restframe(); + const EvtVector4R k = l1.getP4() + l2.getP4(); + const double k2 = k * k; + + const EvtTensor4C T( + dual( EvtGenFunctions::directProd( P, ( 1.0 / k2 ) * k ) ) ); + + double M2 = parent.mass(); M2 *= M2; - double m2 = l1->mass(); + double m2 = l1.mass(); m2 *= m2; - double norm = 1.0 / sqrt( 2 * M2 + 4 * m2 - 4 * m2 * m2 / M2 ); + const double norm = 1.0 / sqrt( 2 * M2 + 4 * m2 - 4 * m2 * m2 / M2 ); - vertex( 0, 0, 0, norm * ( eps0 * T.cont2( l11 ) ) ); - vertex( 0, 0, 1, norm * ( eps0 * T.cont2( l12 ) ) ); - vertex( 0, 1, 0, norm * ( eps0 * T.cont2( l21 ) ) ); - vertex( 0, 1, 1, norm * ( eps0 * T.cont2( l22 ) ) ); + amp.vertex( 0, 0, 0, norm * ( eps0 * T.cont2( l11 ) ) ); + amp.vertex( 0, 0, 1, norm * ( eps0 * T.cont2( l12 ) ) ); + amp.vertex( 0, 1, 0, norm * ( eps0 * T.cont2( l21 ) ) ); + amp.vertex( 0, 1, 1, norm * ( eps0 * T.cont2( l22 ) ) ); - vertex( 1, 0, 0, norm * ( eps1 * T.cont2( l11 ) ) ); - vertex( 1, 0, 1, norm * ( eps1 * T.cont2( l12 ) ) ); - vertex( 1, 1, 0, norm * ( eps1 * T.cont2( l21 ) ) ); - vertex( 1, 1, 1, norm * ( eps1 * T.cont2( l22 ) ) ); + amp.vertex( 1, 0, 0, norm * ( eps1 * T.cont2( l11 ) ) ); + amp.vertex( 1, 0, 1, norm * ( eps1 * T.cont2( l12 ) ) ); + amp.vertex( 1, 1, 0, norm * ( eps1 * T.cont2( l21 ) ) ); + amp.vertex( 1, 1, 1, norm * ( eps1 * T.cont2( l22 ) ) ); - vertex( 2, 0, 0, norm * ( eps2 * T.cont2( l11 ) ) ); - vertex( 2, 0, 1, norm * ( eps2 * T.cont2( l12 ) ) ); - vertex( 2, 1, 0, norm * ( eps2 * T.cont2( l21 ) ) ); - vertex( 2, 1, 1, norm * ( eps2 * T.cont2( l22 ) ) ); + amp.vertex( 2, 0, 0, norm * ( eps2 * T.cont2( l11 ) ) ); + amp.vertex( 2, 0, 1, norm * ( eps2 * T.cont2( l12 ) ) ); + amp.vertex( 2, 1, 0, norm * ( eps2 * T.cont2( l21 ) ) ); + amp.vertex( 2, 1, 1, norm * ( eps2 * T.cont2( l22 ) ) ); return; } diff --git a/test/jsonFiles/VTOSLL__Jpsi_pi0ee.json b/test/jsonFiles/VTOSLL__Jpsi_pi0ee.json new file mode 100644 --- /dev/null +++ b/test/jsonFiles/VTOSLL__Jpsi_pi0ee.json @@ -0,0 +1,26 @@ +{ + "parent" : "J/psi", + "daughters" : ["pi0", "e+", "e-"], + "models" : ["VTOSLL"], + "parameters" : [[]], + "extras" : ["noPhotos"], + "events" : 10000, + "histograms" : [ + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 5.0}, + {"variable" : "mass", "title" : "m(pi0 e+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.14, "xmax" : 3.1}, + {"variable" : "mass", "title" : "m(pi0 e-)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.14, "xmax" : 3.1}, + {"variable" : "mass", "title" : "m(e+e-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.1}, + {"variable" : "pSumSq", "title" : "qSq(e+ e-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 5.0}, + {"variable" : "p", "title" : "p(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.6}, + {"variable" : "pSq", "title" : "pSq(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.6}, + {"variable" : "pz", "title" : "pz(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6}, + {"variable" : "cosHel", "title" : "cosHel(pi0,e+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosHel", "title" : "cosHel(e+,e-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(e+)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(e-)", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "phi(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0} + ], + "outfile" : "VTOSLL__Jpsi_pi0ee.root", + "reference" : "RefVTOSLL__Jpsi_pi0ee.root" +} diff --git a/test/jsonFiles/VTOSLL__Jpsi_pi0mumu.json b/test/jsonFiles/VTOSLL__Jpsi_pi0mumu.json new file mode 100644 --- /dev/null +++ b/test/jsonFiles/VTOSLL__Jpsi_pi0mumu.json @@ -0,0 +1,26 @@ +{ + "parent" : "J/psi", + "daughters" : ["pi0", "mu+", "mu-"], + "models" : ["VTOSLL"], + "parameters" : [[]], + "extras" : ["noPhotos"], + "events" : 10000, + "histograms" : [ + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 5.0}, + {"variable" : "mass", "title" : "m(pi0 mu+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.14, "xmax" : 3.1}, + {"variable" : "mass", "title" : "m(pi0 mu-)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.14, "xmax" : 3.1}, + {"variable" : "mass", "title" : "m(mu+mu-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.1}, + {"variable" : "pSumSq", "title" : "qSq(mu+ mu-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 5.0}, + {"variable" : "p", "title" : "p(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.6}, + {"variable" : "pSq", "title" : "pSq(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.6}, + {"variable" : "pz", "title" : "pz(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6}, + {"variable" : "cosHel", "title" : "cosHel(pi0,mu+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosHel", "title" : "cosHel(mu+,mu-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(mu+)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta(mu-)", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "phi(pi0)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0} + ], + "outfile" : "VTOSLL__Jpsi_pi0mumu.root", + "reference" : "RefVTOSLL__Jpsi_pi0mumu.root" +}