Page MenuHomeHEPForge
Files F7876973

D97.diff
No OneTemporary
Actions

D97.diff
View Options

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 <iostream>
#include <stdlib.h>
@@ -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"
+}

File Metadata

Mime Type
text/plain
Expires
Mon, Nov 18, 2:46 PM (12 h, 24 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804807
Default Alt Text
D97.diff (14 KB)

Event Timeline