Changeset View
Changeset View
Standalone View
Standalone View
EvtGenModels/EvtIntervalDecayAmp.hh
| Show All 38 Lines | |||||
| #include <vector> | #include <vector> | ||||
| // Decay model that uses the "amplitude on an interval" | // Decay model that uses the "amplitude on an interval" | ||||
| // templatization | // templatization | ||||
| template <class T> | template <class T> | ||||
| class EvtIntervalDecayAmp : public EvtDecayAmp { | class EvtIntervalDecayAmp : public EvtDecayAmp { | ||||
| public: | public: | ||||
| EvtIntervalDecayAmp() : _probMax( 0. ), _nScan( 0 ), _fact( 0 ) {} | EvtIntervalDecayAmp() : m_probMax( 0. ), m_nScan( 0 ), m_fact( 0 ) {} | ||||
| EvtIntervalDecayAmp( const EvtIntervalDecayAmp<T>& other ) : | EvtIntervalDecayAmp( const EvtIntervalDecayAmp<T>& other ) : | ||||
| _probMax( other._probMax ), _nScan( other._nScan ), COPY_PTR( _fact ) | m_probMax( other.m_probMax ), m_nScan( other.m_nScan ), COPY_PTR( m_fact ) | ||||
| { | { | ||||
| } | } | ||||
| virtual ~EvtIntervalDecayAmp() { delete _fact; } | virtual ~EvtIntervalDecayAmp() { delete m_fact; } | ||||
| // Initialize model | // Initialize model | ||||
| void init() override | void init() override | ||||
| { | { | ||||
| // Collect model parameters and parse them | // Collect model parameters and parse them | ||||
| vector<std::string> args; | vector<std::string> args; | ||||
| int i; | int i; | ||||
| for ( i = 0; i < getNArg(); i++ ) | for ( i = 0; i < getNArg(); i++ ) | ||||
| args.push_back( getArgStr( i ) ); | args.push_back( getArgStr( i ) ); | ||||
| EvtMultiChannelParser parser; | EvtMultiChannelParser parser; | ||||
| parser.parse( args ); | parser.parse( args ); | ||||
| // Create factory and interval | // Create factory and interval | ||||
| if ( VERBOSE ) | if ( VERBOSE ) | ||||
| EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
| << "Create factory and interval" << std::endl; | << "Create factory and interval" << std::endl; | ||||
| _fact = createFactory( parser ); | m_fact = createFactory( parser ); | ||||
| // Maximum PDF value over the Dalitz plot can be specified, or a scan | // Maximum PDF value over the Dalitz plot can be specified, or a scan | ||||
| // can be performed. | // can be performed. | ||||
| _probMax = parser.pdfMax(); | m_probMax = parser.pdfMax(); | ||||
| _nScan = parser.nScan(); | m_nScan = parser.nScan(); | ||||
| if ( VERBOSE ) | if ( VERBOSE ) | ||||
| EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
| << "Pdf maximum " << _probMax << std::endl; | << "Pdf maximum " << m_probMax << std::endl; | ||||
| if ( VERBOSE ) | if ( VERBOSE ) | ||||
| EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
| << "Scan number " << _nScan << std::endl; | << "Scan number " << m_nScan << std::endl; | ||||
| } | } | ||||
| void initProbMax() override | void initProbMax() override | ||||
| { | { | ||||
| if ( 0 == _nScan ) { | if ( 0 == m_nScan ) { | ||||
| if ( _probMax > 0 ) | if ( m_probMax > 0 ) | ||||
| setProbMax( _probMax ); | setProbMax( m_probMax ); | ||||
| else | else | ||||
| assert( 0 ); | assert( 0 ); | ||||
| } else { | } else { | ||||
| double factor = 1.2; // increase maximum probability by 20% | double factor = 1.2; // increase maximum probability by 20% | ||||
| EvtAmpPdf<T> pdf( *_fact->getAmp() ); | EvtAmpPdf<T> pdf( *m_fact->getAmp() ); | ||||
| EvtPdfSum<T>* pc = _fact->getPC(); | EvtPdfSum<T>* pc = m_fact->getPC(); | ||||
| EvtPdfDiv<T> pdfdiv( pdf, *pc ); | EvtPdfDiv<T> pdfdiv( pdf, *pc ); | ||||
| printf( "Sampling %d points to find maximum\n", _nScan ); | printf( "Sampling %d points to find maximum\n", m_nScan ); | ||||
| EvtPdfMax<T> x = pdfdiv.findMax( *pc, _nScan ); | EvtPdfMax<T> x = pdfdiv.findMax( *pc, m_nScan ); | ||||
| _probMax = factor * x.value(); | m_probMax = factor * x.value(); | ||||
| printf( "Found maximum %f\n", x.value() ); | printf( "Found maximum %f\n", x.value() ); | ||||
| printf( "Increase to %f\n", _probMax ); | printf( "Increase to %f\n", m_probMax ); | ||||
| setProbMax( _probMax ); | setProbMax( m_probMax ); | ||||
| } | } | ||||
| } | } | ||||
| void decay( EvtParticle* p ) override | void decay( EvtParticle* p ) override | ||||
| { | { | ||||
| // Set things up in most general way | // Set things up in most general way | ||||
| static EvtId B0 = EvtPDL::getId( "B0" ); | static EvtId B0 = EvtPDL::getId( "B0" ); | ||||
| static EvtId B0B = EvtPDL::getId( "anti-B0" ); | static EvtId B0B = EvtPDL::getId( "anti-B0" ); | ||||
| double t; | double t; | ||||
| EvtId other_b; | EvtId other_b; | ||||
| EvtComplex ampl( 0., 0. ); | EvtComplex ampl( 0., 0. ); | ||||
| // Sample using pole-compensator pdf | // Sample using pole-compensator pdf | ||||
| EvtPdfSum<T>* pc = getPC(); | EvtPdfSum<T>* pc = getPC(); | ||||
| _x = pc->randomPoint(); | m_x = pc->randomPoint(); | ||||
| if ( _fact->isCPModel() ) { | if ( m_fact->isCPModel() ) { | ||||
| // Time-dependent Dalitz plot changes | // Time-dependent Dalitz plot changes | ||||
| // Dec 2005 (ddujmic@slac.stanford.edu) | // Dec 2005 (ddujmic@slac.stanford.edu) | ||||
| EvtComplex A = _fact->getAmp()->evaluate( _x ); | EvtComplex A = m_fact->getAmp()->evaluate( m_x ); | ||||
| EvtComplex Abar = _fact->getAmpConj()->evaluate( _x ); | EvtComplex Abar = m_fact->getAmpConj()->evaluate( m_x ); | ||||
| EvtCPUtil::getInstance()->OtherB( p, t, other_b ); | EvtCPUtil::getInstance()->OtherB( p, t, other_b ); | ||||
| double dm = _fact->dm(); | double dm = m_fact->dm(); | ||||
| double mixAmpli = _fact->mixAmpli(); | double mixAmpli = m_fact->mixAmpli(); | ||||
| double mixPhase = _fact->mixPhase(); | double mixPhase = m_fact->mixPhase(); | ||||
| EvtComplex qoverp( cos( mixPhase ) * mixAmpli, | EvtComplex qoverp( cos( mixPhase ) * mixAmpli, | ||||
| sin( mixPhase ) * mixAmpli ); | sin( mixPhase ) * mixAmpli ); | ||||
| EvtComplex poverq( cos( mixPhase ) / mixAmpli, | EvtComplex poverq( cos( mixPhase ) / mixAmpli, | ||||
| -sin( mixPhase ) / mixAmpli ); | -sin( mixPhase ) / mixAmpli ); | ||||
| if ( other_b == B0B ) | if ( other_b == B0B ) | ||||
| ampl = A * cos( dm * t / ( 2 * EvtConst::c ) ) + | ampl = A * cos( dm * t / ( 2 * EvtConst::c ) ) + | ||||
| EvtComplex( 0., 1. ) * Abar * | EvtComplex( 0., 1. ) * Abar * | ||||
| sin( dm * t / ( 2 * EvtConst::c ) ) * qoverp; | sin( dm * t / ( 2 * EvtConst::c ) ) * qoverp; | ||||
| if ( other_b == B0 ) | if ( other_b == B0 ) | ||||
| ampl = Abar * cos( dm * t / ( 2 * EvtConst::c ) ) + | ampl = Abar * cos( dm * t / ( 2 * EvtConst::c ) ) + | ||||
| EvtComplex( 0., 1. ) * A * | EvtComplex( 0., 1. ) * A * | ||||
| sin( dm * t / ( 2 * EvtConst::c ) ) * poverq; | sin( dm * t / ( 2 * EvtConst::c ) ) * poverq; | ||||
| } else { | } else { | ||||
| ampl = amplNonCP( _x ); | ampl = amplNonCP( m_x ); | ||||
| } | } | ||||
| // Pole-compensate | // Pole-compensate | ||||
| double comp = sqrt( pc->evaluate( _x ) ); | double comp = sqrt( pc->evaluate( m_x ) ); | ||||
| assert( comp > 0 ); | assert( comp > 0 ); | ||||
| vertex( ampl / comp ); | vertex( ampl / comp ); | ||||
| // Now generate random angles, rotate and setup | // Now generate random angles, rotate and setup | ||||
| // the daughters | // the daughters | ||||
| std::vector<EvtVector4R> v = initDaughters( _x ); | std::vector<EvtVector4R> v = initDaughters( m_x ); | ||||
| size_t N = p->getNDaug(); | size_t N = p->getNDaug(); | ||||
| if ( v.size() != N ) { | if ( v.size() != N ) { | ||||
| EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
| << "Number of daughters " << N << std::endl; | << "Number of daughters " << N << std::endl; | ||||
| EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
| << "Momentum vector size " << v.size() << std::endl; | << "Momentum vector size " << v.size() << std::endl; | ||||
| assert( 0 ); | assert( 0 ); | ||||
| } | } | ||||
| for ( size_t i = 0; i < N; i++ ) { | for ( size_t i = 0; i < N; i++ ) { | ||||
| p->getDaug( i )->init( getDaugs()[i], v[i] ); | p->getDaug( i )->init( getDaugs()[i], v[i] ); | ||||
| } | } | ||||
| } | } | ||||
| virtual EvtAmpFactory<T>* createFactory( | virtual EvtAmpFactory<T>* createFactory( | ||||
| const EvtMultiChannelParser& parser ) = 0; | const EvtMultiChannelParser& parser ) = 0; | ||||
| virtual std::vector<EvtVector4R> initDaughters( const T& p ) const = 0; | virtual std::vector<EvtVector4R> initDaughters( const T& p ) const = 0; | ||||
| // provide access to the decay point and to the amplitude of any decay point. | // provide access to the decay point and to the amplitude of any decay point. | ||||
| // this is used by EvtBtoKD3P: | // this is used by EvtBtoKD3P: | ||||
| const T& x() const { return _x; } | const T& x() const { return m_x; } | ||||
| EvtComplex amplNonCP( const T& x ) | EvtComplex amplNonCP( const T& x ) | ||||
| { | { | ||||
| return _fact->getAmp()->evaluate( x ); | return m_fact->getAmp()->evaluate( x ); | ||||
| } | } | ||||
| EvtPdfSum<T>* getPC() { return _fact->getPC(); } | EvtPdfSum<T>* getPC() { return m_fact->getPC(); } | ||||
| protected: | protected: | ||||
| double _probMax; // Maximum probability | double m_probMax; // Maximum probability | ||||
| int _nScan; // Number of points for max prob DP scan | int m_nScan; // Number of points for max prob DP scan | ||||
| T _x; // Decay point | T m_x; // Decay point | ||||
| EvtAmpFactory<T>* _fact; // factory | EvtAmpFactory<T>* m_fact; // factory | ||||
| }; | }; | ||||
| #endif | #endif | ||||