include/EvtGenModels/EvtIntervalDecayAmp.hh

0001
0002 /***********************************************************************
0003 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
0004 *                                                                      *
0005 * This file is part of EvtGen.                                         *
0006 *                                                                      *
0007 * EvtGen is free software: you can redistribute it and/or modify       *
0008 * it under the terms of the GNU General Public License as published by *
0009 * the Free Software Foundation, either version 3 of the License, or    *
0010 * (at your option) any later version.                                  *
0011 *                                                                      *
0012 * EvtGen is distributed in the hope that it will be useful,            *
0013 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
0014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
0015 * GNU General Public License for more details.                         *
0016 *                                                                      *
0017 * You should have received a copy of the GNU General Public License    *
0018 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
0019 ***********************************************************************/
0020
0021 #ifndef EVT_INTERVAL_DECAY_AMP
0022 #define EVT_INTERVAL_DECAY_AMP
0023
0024 #define VERBOSE true
0025 #include "EvtGenBase/EvtAmpFactory.hh"
0026 #include "EvtGenBase/EvtAmpPdf.hh"
0027 #include "EvtGenBase/EvtCPUtil.hh"
0028 #include "EvtGenBase/EvtCyclic3.hh"
0029 #include "EvtGenBase/EvtDecayAmp.hh"
0030 #include "EvtGenBase/EvtMacros.hh"
0031 #include "EvtGenBase/EvtMultiChannelParser.hh"
0032 #include "EvtGenBase/EvtPDL.hh"
0033 #include "EvtGenBase/EvtParticle.hh"
0034 #include "EvtGenBase/EvtPdf.hh"
0035 #include "EvtGenBase/EvtReport.hh"
0036
0037 #include <iostream>
0038 #include <string>
0039 #include <vector>
0040
0041 // Decay model that uses the "amplitude on an interval"
0042 // templatization
0043
0044 template <class T>
0045 class EvtIntervalDecayAmp : public EvtDecayAmp {
0046   public:
0047     EvtIntervalDecayAmp() : _probMax( 0. ), _nScan( 0 ), _fact( 0 ) {}
0048
0049     EvtIntervalDecayAmp( const EvtIntervalDecayAmp<T>& other ) :
0050         _probMax( other._probMax ), _nScan( other._nScan ), COPY_PTR( _fact )
0051     {
0052     }
0053
0054     virtual ~EvtIntervalDecayAmp() { delete _fact; }
0055
0056     // Initialize model
0057
0058     void init() override
0059     {
0060         // Collect model parameters and parse them
0061
0062         vector<std::string> args;
0063         int i;
0064         for ( i = 0; i < getNArg(); i++ )
0065             args.push_back( getArgStr( i ) );
0066         EvtMultiChannelParser parser;
0067         parser.parse( args );
0068
0069         // Create factory and interval
0070
0071         if ( VERBOSE )
0072             EvtGenReport( EVTGEN_INFO, "EvtGen" )
0073                 << "Create factory and interval" << std::endl;
0074         _fact = createFactory( parser );
0075
0076         // Maximum PDF value over the Dalitz plot can be specified, or a scan
0077         // can be performed.
0078
0079         _probMax = parser.pdfMax();
0080         _nScan = parser.nScan();
0081         if ( VERBOSE )
0082             EvtGenReport( EVTGEN_INFO, "EvtGen" )
0083                 << "Pdf maximum " << _probMax << std::endl;
0084         if ( VERBOSE )
0085             EvtGenReport( EVTGEN_INFO, "EvtGen" )
0086                 << "Scan number " << _nScan << std::endl;
0087     }
0088
0089     void initProbMax() override
0090     {
0091         if ( 0 == _nScan ) {
0092             if ( _probMax > 0 )
0093                 setProbMax( _probMax );
0094             else
0095                 assert( 0 );
0096         } else {
0097             double factor = 1.2;    // increase maximum probability by 20%
0098             EvtAmpPdf<T> pdf( *_fact->getAmp() );
0099             EvtPdfSum<T>* pc = _fact->getPC();
0100             EvtPdfDiv<T> pdfdiv( pdf, *pc );
0101             printf( "Sampling %d points to find maximum\n", _nScan );
0102             EvtPdfMax<T> x = pdfdiv.findMax( *pc, _nScan );
0103             _probMax = factor * x.value();
0104             printf( "Found maximum %f\n", x.value() );
0105             printf( "Increase to   %f\n", _probMax );
0106             setProbMax( _probMax );
0107         }
0108     }
0109
0110     void decay( EvtParticle* p ) override
0111     {
0112         // Set things up in most general way
0113
0114         static EvtId B0 = EvtPDL::getId( "B0" );
0115         static EvtId B0B = EvtPDL::getId( "anti-B0" );
0116         double t;
0117         EvtId other_b;
0118         EvtComplex ampl( 0., 0. );
0119
0120         // Sample using pole-compensator pdf
0121
0122         EvtPdfSum<T>* pc = getPC();
0123         _x = pc->randomPoint();
0124
0125         if ( _fact->isCPModel() ) {
0126             // Time-dependent Dalitz plot changes
0127             // Dec 2005 (ddujmic@slac.stanford.edu)
0128
0129             EvtComplex A = _fact->getAmp()->evaluate( _x );
0130             EvtComplex Abar = _fact->getAmpConj()->evaluate( _x );
0131
0132             EvtCPUtil::getInstance()->OtherB( p, t, other_b );
0133
0134             double dm = _fact->dm();
0135             double mixAmpli = _fact->mixAmpli();
0136             double mixPhase = _fact->mixPhase();
0137             EvtComplex qoverp( cos( mixPhase ) * mixAmpli,
0138                                sin( mixPhase ) * mixAmpli );
0139             EvtComplex poverq( cos( mixPhase ) / mixAmpli,
0140                                -sin( mixPhase ) / mixAmpli );
0141
0142             if ( other_b == B0B )
0143                 ampl = A * cos( dm * t / ( 2 * EvtConst::c ) ) +
0144                        EvtComplex( 0., 1. ) * Abar *
0145                            sin( dm * t / ( 2 * EvtConst::c ) ) * qoverp;
0146             if ( other_b == B0 )
0147                 ampl = Abar * cos( dm * t / ( 2 * EvtConst::c ) ) +
0148                        EvtComplex( 0., 1. ) * A *
0149                            sin( dm * t / ( 2 * EvtConst::c ) ) * poverq;
0150
0151         } else {
0152             ampl = amplNonCP( _x );
0153         }
0154
0155         // Pole-compensate
0156
0157         double comp = sqrt( pc->evaluate( _x ) );
0158         assert( comp > 0 );
0159         vertex( ampl / comp );
0160
0161         // Now generate random angles, rotate and setup
0162         // the daughters
0163
0164         std::vector<EvtVector4R> v = initDaughters( _x );
0165
0166         size_t N = p->getNDaug();
0167         if ( v.size() != N ) {
0168             EvtGenReport( EVTGEN_INFO, "EvtGen" )
0169                 << "Number of daughters " << N << std::endl;
0170             EvtGenReport( EVTGEN_INFO, "EvtGen" )
0171                 << "Momentum vector size " << v.size() << std::endl;
0172             assert( 0 );
0173         }
0174
0175         for ( size_t i = 0; i < N; i++ ) {
0176             p->getDaug( i )->init( getDaugs()[i], v[i] );
0177         }
0178     }
0179
0180     virtual EvtAmpFactory<T>* createFactory(
0181         const EvtMultiChannelParser& parser ) = 0;
0182     virtual std::vector<EvtVector4R> initDaughters( const T& p ) const = 0;
0183
0184     // provide access to the decay point and to the amplitude of any decay point.
0185     // this is used by EvtBtoKD3P:
0186     const T& x() const { return _x; }
0187     EvtComplex amplNonCP( const T& x )
0188     {
0189         return _fact->getAmp()->evaluate( x );
0190     }
0191     EvtPdfSum<T>* getPC() { return _fact->getPC(); }
0192
0193   protected:
0194     double _probMax;    // Maximum probability
0195     int _nScan;         // Number of points for max prob DP scan
0196     T _x;               // Decay point
0197
0198     EvtAmpFactory<T>* _fact;    // factory
0199 };
0200
0201 #endif