diff --git a/EvtGenModels/EvtDToKpienu.hh b/EvtGenModels/EvtDToKpienu.hh
index d836c5d..a96a1c0 100644
--- a/EvtGenModels/EvtDToKpienu.hh
+++ b/EvtGenModels/EvtDToKpienu.hh
@@ -1,98 +1,96 @@
 #ifndef EVTDTOKPIENU_HH
 #define EVTDTOKPIENU_HH
 
 #include <array>
 #include "EvtGenBase/EvtDecayProb.hh"
 #include "EvtGenBase/EvtComplex.hh"
 #include "EvtGenBase/EvtVector4R.hh"
 
 class EvtParticle;
 
 class EvtDToKpienu : public EvtDecayProb {
 public:
-  EvtDToKpienu() {}
-  virtual ~EvtDToKpienu();
 
-  std::string getName();
-  EvtDecayBase* clone();
+  std::string getName() override;
+  EvtDecayBase* clone() override;
 
-  void init();
-  void initProbMax();
-  void decay( EvtParticle* p );
+  void init() override;
+  void initProbMax() override;
+  void decay( EvtParticle* p ) override;
 
 private:
   void KinVGen( const EvtVector4R& vp4_K, const EvtVector4R& vp4_Pi,
                 const EvtVector4R& vp4_Lep, const EvtVector4R& vp4_Nu,
                 const int charm, double& m2, double& q2, double& cosV, double& cosL,
-                double& chi );
+                double& chi ) const;
   double calPDF( const double m2, const double q2, const double cosV,
-                 const double cosL, const double chi );
+                 const double cosL, const double chi ) const;
   void ResonanceP( const double m, const double q, const double mV,
                    const double mA, const double V_0, const double A1_0,
                    const double A2_0, const double m0, const double width0,
                    const double rBW, double& amplitude, double& delta,
-                   EvtComplex& F11, EvtComplex& F21, EvtComplex& F31 );
+                   EvtComplex& F11, EvtComplex& F21, EvtComplex& F31 ) const;
   void NRS( const double m, const double q, const double rS, const double rS1,
             const double a_delta, const double b_delta, const double mA,
             const double m0, const double width0, double& amplitude,
-            double& delta, EvtComplex& F10 );
+            double& delta, EvtComplex& F10 ) const;
   void ResonanceD( const double m, const double q, const double mV, const double mA, const double TV_0,
                    const double T1_0, const double T2_0, const double m0, const double width0,
                    const double rBW, double& amplitude, double& delta,
-                   EvtComplex& F12, EvtComplex& F22, EvtComplex& F32 );
-  double getPStar( const double m, const double m1, const double m2 );
+                   EvtComplex& F12, EvtComplex& F22, EvtComplex& F32 ) const;
+  double getPStar( const double m, const double m1, const double m2 ) const;
   double getF1( const double m, const double m0, const double m_c1,
-                const double m_c2, const double rBW );
+                const double m_c2, const double rBW ) const;
   double getF2( const double m, const double m0, const double m_c1,
-                const double m_c2, const double rBW );
+                const double m_c2, const double rBW ) const;
   double getWidth0( const double m, const double m0, const double m_c1, const double m_c2,
-                    const double width0 );
+                    const double width0 ) const;
   double getWidth1( const double m, const double m0, const double m_c1, const double m_c2,
-                    const double width0, const double rBW );
+                    const double width0, const double rBW ) const;
   double getWidth2( const double m, const double m0, const double m_c1, const double m_c2,
-                    const double width0, const double rBW );
-  EvtComplex getCoef( const double rho, const double phi );
+                    const double width0, const double rBW ) const;
+  EvtComplex getCoef( const double rho, const double phi ) const;
 
   int nAmps;
   std::array<int,5> type;
 
   double rS;
   double rS1;
   double a_delta;
   double b_delta;
   double m0_1430_S;
   double width0_1430_S;
 
   double mV;
   double mA;
   double V_0;
   double A1_0;
   double A2_0;
   double m0;
   double width0;
   double rBW;
   double rho;
   double phi;
   double m0_1410;
   double width0_1410;
   double rho_1410;
   double phi_1410;
   double TV_0;
   double T1_0;
   double T2_0;
   double m0_1430;
   double width0_1430;
   double rho_1430;
   double phi_1430;
 
   double mD;
   double mPi;
   double mK;
   double Pi;
   double root2;
   double root2d3;
   double root1d2;
   double root3d2;
 };
 
 #endif
diff --git a/src/EvtGenModels/EvtDToKpienu.cpp b/src/EvtGenModels/EvtDToKpienu.cpp
index 06ba5dd..3874d5b 100644
--- a/src/EvtGenModels/EvtDToKpienu.cpp
+++ b/src/EvtGenModels/EvtDToKpienu.cpp
@@ -1,498 +1,519 @@
+//--------------------------------------------------------------------------
+//
+// Environment:
+//      This software is part of the EvtGen package developed jointly
+//      for the BaBar and CLEO collaborations. If you use all or part
+//      of it, please give an appropriate acknowledgement.
+//
+// Copyright Information: See EvtGen/COPYRIGHT
+//      Copyright (C) 2000      Caltech, UCSB
+//
+// Module: EvtDToKpienu.cpp
+//
+// Description: Decay model for D- --> Kpienu decays based on PWA analysis
+//              in Phys. Rev. D94 (2016) 0032001 
+//
+// Modification history:
+//
+// Liaoyuan Dong                     Aug 28 2019    
+//
+//------------------------------------------------------------------------
+
 #include "EvtGenModels/EvtDToKpienu.hh"
 
 #include "EvtGenBase/EvtComplex.hh"
 #include "EvtGenBase/EvtDecayTable.hh"
 #include "EvtGenBase/EvtGenKine.hh"
 #include "EvtGenBase/EvtPDL.hh"
 #include "EvtGenBase/EvtParticle.hh"
 #include "EvtGenBase/EvtPatches.hh"
 #include "EvtGenBase/EvtReport.hh"
 
-#include <stdlib.h>
-
-EvtDToKpienu::~EvtDToKpienu()
-{
-}
 
 std::string EvtDToKpienu::getName()
 {
   return "DToKpienu";
 }
 
 EvtDecayBase* EvtDToKpienu::clone()
 {
   return new EvtDToKpienu;
 }
 
 void EvtDToKpienu::init()
 {
   checkNArg( 0 );
   checkNDaug( 4 );
   checkSpinParent( EvtSpinType::SCALAR );
   checkSpinDaughter( 0, EvtSpinType::SCALAR );
   checkSpinDaughter( 1, EvtSpinType::SCALAR );
 
-  std::cout << "EvtDToKpienu ==> Initialization !" << std::endl;
+  EvtGenReport(EVTGEN_INFO,"EvtGen") << "EvtDToKpienu ==> Initialization !" << std::endl;
   nAmps = 2;
 
   rS            = -11.57;    // S-wave
   rS1           = 0.08;
   a_delta       = 1.94;
   b_delta       = -0.81;
   m0_1430_S     = 1.425;
   width0_1430_S = 0.270;
   type[0]       = 0;
 
   mV   = 1.81;
   mA   = 2.61;
   V_0  = 1.411;
   A1_0 = 1;
   A2_0 = 0.788;
 
   m0      = 0.8946;    // P-wave K*
   width0  = 0.04642;
   rBW     = 3.07;
   rho     = 1;
   phi     = 0;
   type[1] = 1;
 
   m0_1410     = 1.414;    // P-wave K*(1410)
   width0_1410 = 0.232;
   rho_1410    = 0.1;
   phi_1410    = 0.;
   type[2]     = 2;
 
   TV_0        = 1;    // D-wave K*2(1430)
   T1_0        = 1;
   T2_0        = 1;
   m0_1430     = 1.4324;
   width0_1430 = 0.109;
   rho_1430    = 15;
   phi_1430    = 0;
   type[3]     = 3;
 
   mD      = 1.86962;
   mPi     = 0.13957;
   mK      = 0.49368;
   Pi      = atan2( 0.0, -1.0 );
   root2   = sqrt( 2. );
   root2d3 = sqrt( 2. / 3 );
   root1d2 = sqrt( 0.5 );
   root3d2 = sqrt( 1.5 );
 }
 
 void EvtDToKpienu::initProbMax()
 {
-  setProbMax( 22750.0 );
-}
-
-void EvtDToKpienu::decay( EvtParticle* p )
-{
   /*
+   * This piece of code could in principle be used to calculate maximum
+   * probablity on fly. But as it uses high number of points and model
+   * deals with single final state, we keep hardcoded number for now rather
+   * than adapting code to work here.
+
      double maxprob = 0.0;
      for(int ir=0;ir<=60000000;ir++){
         p->initializePhaseSpace(getNDaug(),getDaugs());
         EvtVector4R _K  = p->getDaug(0)->getP4();
         EvtVector4R _pi = p->getDaug(1)->getP4();
         EvtVector4R _e  = p->getDaug(2)->getP4();
         EvtVector4R _nu = p->getDaug(3)->getP4();
         int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
         int charm;
         if(pid == -321) charm = 1;
         else charm = -1;
         double m2, q2, cosV, cosL, chi;
         KinVGen(_K, _pi, _e, _nu, charm, m2, q2, cosV, cosL, chi);
         double _prob = calPDF(m2, q2, cosV, cosL, chi);
         if(_prob>maxprob) {
             maxprob=_prob;
-            std::cout << "Max PDF = " << ir << " charm= " << charm << " prob= "
+            EvtGenReport(EVTGEN_INFO,"EvtGen") << "Max PDF = " << ir << " charm= " << charm << " prob= "
      << _prob << std::endl;
         }
      }
-     std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
+     EvtGenReport(EVTGEN_INFO,"EvtGen") << "Max!!!!!!!!!!! " << maxprob<< std::endl;
   */
+  setProbMax( 22750.0 );
+}
+
+void EvtDToKpienu::decay( EvtParticle* p )
+{
 
   p->initializePhaseSpace( getNDaug(), getDaugs() );
   EvtVector4R K  = p->getDaug( 0 )->getP4();
   EvtVector4R pi = p->getDaug( 1 )->getP4();
   EvtVector4R e  = p->getDaug( 2 )->getP4();
   EvtVector4R nu = p->getDaug( 3 )->getP4();
 
   int pid = EvtPDL::getStdHep( p->getDaug( 0 )->getId() );
   int charm;
   if ( pid == -321 ) {
     charm = 1;
   } else {
     charm = -1;
   }
   double m2, q2, cosV, cosL, chi;
   KinVGen( K, pi, e, nu, charm, m2, q2, cosV, cosL, chi );
   double prob = calPDF( m2, q2, cosV, cosL, chi );
   setProb( prob );
   return;
 }
 
 void EvtDToKpienu::KinVGen( const EvtVector4R& vp4_K, const EvtVector4R& vp4_Pi,
                             const EvtVector4R& vp4_Lep,
                             const EvtVector4R& vp4_Nu, const int charm, double& m2,
                             double& q2, double& cosV, double& cosL,
-                            double& chi )
+                            double& chi ) const
 {
   EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
   EvtVector4R vp4_W   = vp4_Lep + vp4_Nu;
 
   m2 = vp4_KPi.mass2();
   q2 = vp4_W.mass2();
 
   EvtVector4R boost;
   boost.set( vp4_KPi.get( 0 ), -vp4_KPi.get( 1 ), -vp4_KPi.get( 2 ),
              -vp4_KPi.get( 3 ) );
   EvtVector4R vp4_Kp = boostTo( vp4_K, boost );
   cosV = vp4_Kp.dot( vp4_KPi ) / ( vp4_Kp.d3mag() * vp4_KPi.d3mag() );
 
   boost.set( vp4_W.get( 0 ), -vp4_W.get( 1 ), -vp4_W.get( 2 ),
              -vp4_W.get( 3 ) );
   EvtVector4R vp4_Lepp = boostTo( vp4_Lep, boost );
   cosL = vp4_Lepp.dot( vp4_W ) / ( vp4_Lepp.d3mag() * vp4_W.d3mag() );
 
   EvtVector4R V = vp4_KPi / vp4_KPi.d3mag();
   EvtVector4R C = vp4_Kp.cross( V );
   C /= C.d3mag();
   EvtVector4R D = vp4_Lepp.cross( V );
   D /= D.d3mag();
   double sinx = C.cross( V ).dot( D );
   double cosx = C.dot( D );
   chi         = sinx > 0 ? acos( cosx ) : -acos( cosx );
   if ( charm == -1 )
     chi = -chi;
 }
 
 double EvtDToKpienu::calPDF( const double m2, const double q2, const double cosV, const double cosL,
-                             const double chi )
+                             const double chi ) const
 {
   double m = sqrt( m2 );
   double q = sqrt( q2 );
 
   // begin to calculate form factor
   EvtComplex F10( 0.0, 0.0 );
   EvtComplex F11( 0.0, 0.0 );
   EvtComplex F21( 0.0, 0.0 );
   EvtComplex F31( 0.0, 0.0 );
   EvtComplex F12( 0.0, 0.0 );
   EvtComplex F22( 0.0, 0.0 );
   EvtComplex F32( 0.0, 0.0 );
 
   EvtComplex f10( 0.0, 0.0 );
   EvtComplex f11( 0.0, 0.0 );
   EvtComplex f21( 0.0, 0.0 );
   EvtComplex f31( 0.0, 0.0 );
   EvtComplex f12( 0.0, 0.0 );
   EvtComplex f22( 0.0, 0.0 );
   EvtComplex f32( 0.0, 0.0 );
   EvtComplex coef( 0.0, 0.0 );
   double amplitude_temp, delta_temp;
 
   for ( int index = 0; index < nAmps; index++ ) {
     switch ( type[index] ) {
       case 0:    // calculate form factor of S wave
       {
         NRS( m, q, rS, rS1, a_delta, b_delta, mA, m0_1430_S, width0_1430_S,
              amplitude_temp, delta_temp, f10 );
         F10              = F10 + f10;
         break;
       }
       case 1:    // calculate form factor of P wave (K*)
       {
         ResonanceP( m, q, mV, mA, V_0, A1_0, A2_0, m0, width0, rBW,
                     amplitude_temp, delta_temp, f11, f21, f31 );
         coef             = getCoef( rho, phi );
         F11              = F11 + coef * f11;
         F21              = F21 + coef * f21;
         F31              = F31 + coef * f31;
         break;
       }
       case 2:    // calculate form factor of P wave (K*(1410))
       {
         ResonanceP( m, q, mV, mA, V_0, A1_0, A2_0, m0_1410, width0_1410, rBW,
                     amplitude_temp, delta_temp, f11, f21, f31 );
         coef             = getCoef( rho_1410, phi_1410 );
         F11              = F11 + coef * f11;
         F21              = F21 + coef * f21;
         F31              = F31 + coef * f31;
         break;
       }
       case 3:    // calculate form factor of D wave
       {
         ResonanceD( m, q, mV, mA, TV_0, T1_0, T2_0, m0_1430, width0_1430, rBW,
                     amplitude_temp, delta_temp, f12, f22, f32 );
         coef             = getCoef( rho_1430, phi_1430 );
         F12              = F12 + coef * f12;
         F22              = F22 + coef * f22;
         F32              = F32 + coef * f32;
         break;
       }
       default: {
-        std::cout << "No such form factor type!!!" << std::endl;
+        EvtGenReport(EVTGEN_ERROR,"EvtGen") << "No such form factor type!!!" << std::endl;
         break;
       }
     }
   }
 
   // begin to calculate pdf value
   double I, I1, I2, I3, I4, I5, I6, I7, I8, I9;
 
   double cosV2 = cosV * cosV;
   double sinV  = sqrt( 1.0 - cosV2 );
   double sinV2 = sinV * sinV;
 
   EvtComplex F1 = F10 + F11 * cosV + F12 * ( 1.5 * cosV2 - 0.5 );
   EvtComplex F2 = F21 * root1d2 + F22 * cosV * root3d2;
   EvtComplex F3 = F31 * root1d2 + F32 * cosV * root3d2;
 
   I1 = 0.25 * ( abs2( F1 ) + 1.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
   I2 = -0.25 * ( abs2( F1 ) - 0.5 * sinV2 * ( abs2( F2 ) + abs2( F3 ) ) );
   I3 = -0.25 * ( abs2( F2 ) - abs2( F3 ) ) * sinV2;
   I4 = real( conj( F1 ) * F2 ) * sinV * 0.5;
   I5 = real( conj( F1 ) * F3 ) * sinV;
   I6 = real( conj( F2 ) * F3 ) * sinV2;
   I7 = imag( conj( F2 ) * F1 ) * sinV;
   I8 = imag( conj( F3 ) * F1 ) * sinV * 0.5;
   I9 = imag( conj( F3 ) * F2 ) * sinV2 * ( -0.5 );
 
   double coschi  = cos( chi );
   double sinchi  = sin( chi );
   double sin2chi = 2.0 * sinchi * coschi;
   double cos2chi = 1.0 - 2.0 * sinchi * sinchi;
 
   double sinL  = sqrt( 1. - cosL * cosL );
   double sinL2 = sinL * sinL;
   double sin2L = 2.0 * sinL * cosL;
   double cos2L = 1.0 - 2.0 * sinL2;
 
   I = I1 + I2 * cos2L + I3 * sinL2 * cos2chi + I4 * sin2L * coschi +
       I5 * sinL * coschi + I6 * cosL + I7 * sinL * sinchi +
       I8 * sin2L * sinchi + I9 * sinL2 * sin2chi;
   return I;
 }
 
 void EvtDToKpienu::ResonanceP( const double m, const double q, const double mV,
                                const double mA, const double V_0,
                                const double A1_0, const double A2_0,
                                const double m0, const double width0,
                                const double rBW, double& amplitude,
                                double& delta, EvtComplex& F11, EvtComplex& F21,
-                               EvtComplex& F31 )
+                               EvtComplex& F31 ) const
 {
   double pKPi   = getPStar( mD, m, q );
   double mD2    = mD * mD;
   double m2     = m * m;
   double m02    = m0 * m0;
   double q2     = q * q;
   double mV2    = mV * mV;
   double mA2    = mA * mA;
   double summDm = mD + m;
   double V      = V_0 / ( 1.0 - q2 / ( mV2 ) );
   double A1     = A1_0 / ( 1.0 - q2 / ( mA2 ) );
   double A2     = A2_0 / ( 1.0 - q2 / ( mA2 ) );
   double A      = summDm * A1;
   double B      = 2.0 * mD * pKPi / summDm * V;
 
   // construct the helicity form factor
   double H0 = 0.5 / ( m * q ) *
               ( ( mD2 - m2 - q2 ) * summDm * A1 -
                 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
   double Hp = A - B;
   double Hm = A + B;
 
   // calculate alpha
   double B_Kstar = 2. / 3.;    // B_Kstar = Br(Kstar(892)->k pi)
   double pStar0  = getPStar( m0, mPi, mK );
   double alpha   = sqrt( 3. * Pi * B_Kstar / ( pStar0 * width0 ) );
 
   // construct amplitudes of (non)resonance
   double F     = getF1( m, m0, mPi, mK, rBW );
   double width = getWidth1( m, m0, mPi, mK, width0, rBW );
 
   EvtComplex C( m0 * width0 * F, 0.0 );
   double AA      = m02 - m2;
   double BB      = -m0 * width;
   EvtComplex amp = C / EvtComplex( AA, BB );
   amplitude      = abs( amp );
   delta          = atan2( imag( amp ), real( amp ) );
 
   double alpham2 = alpha * 2.0;
   F11            = amp * alpham2 * q * H0 * root2;
   F21            = amp * alpham2 * q * ( Hp + Hm );
   F31            = amp * alpham2 * q * ( Hp - Hm );
 }
 
 void EvtDToKpienu::NRS( const double m, const double q, const double rS,
                         const double rS1, const double a_delta,
                         const double b_delta, const double mA, const double m0,
                         const double width0, double& amplitude, double& delta,
-                        EvtComplex& F10 )
+                        EvtComplex& F10 ) const
 {
   static const double tmp = ( mK + mPi ) * ( mK + mPi );
 
   double m2            = m * m;
   double q2            = q * q;
   double mA2           = mA * mA;
   double pKPi          = getPStar( mD, m, q );
   double m_K0_1430     = m0;
   double width_K0_1430 = width0;
   double m2_K0_1430    = m_K0_1430 * m_K0_1430;
   double width         = getWidth0( m, m_K0_1430, mPi, mK, width_K0_1430 );
 
   // calculate modul of the amplitude
   double x, Pm;
   if ( m < m_K0_1430 ) {
     x  = sqrt( m2 / tmp - 1.0 );
     Pm = 1.0 + rS1 * x;
   } else {
     x  = sqrt( m2_K0_1430 / tmp - 1.0 );
     Pm = 1.0 + rS1 * x;
     Pm *= m_K0_1430 * width_K0_1430 /
           sqrt( ( m2_K0_1430 - m2 ) * ( m2_K0_1430 - m2 ) +
                 m2_K0_1430 * width * width );
   }
 
   // calculate phase of the amplitude
   double pStar = getPStar( m, mPi, mK );
   double delta_bg =
     atan( 2. * a_delta * pStar / ( 2. + a_delta * b_delta * pStar * pStar ) );
   delta_bg = ( delta_bg > 0 ) ? delta_bg : ( delta_bg + Pi );
 
   double delta_K0_1430 = atan( m_K0_1430 * width / ( m2_K0_1430 - m2 ) );
   delta_K0_1430 =
     ( delta_K0_1430 > 0 ) ? delta_K0_1430 : ( delta_K0_1430 + Pi );
   delta = delta_bg + delta_K0_1430;
 
   EvtComplex ci( cos( delta ), sin( delta ) );
   EvtComplex amp = ci * rS * Pm;
   amplitude      = rS * Pm;
 
   F10 = amp * pKPi * mD / ( 1. - q2 / mA2 );
 }
 
 void EvtDToKpienu::ResonanceD( const double m, const double q, const double mV,
                                const double mA, const double TV_0,
                                const double T1_0, const double T2_0,
                                const double m0, const double width0,
                                const double rBW, double& amplitude,
                                double& delta, EvtComplex& F12, EvtComplex& F22,
-                               EvtComplex& F32 )
+                               EvtComplex& F32 ) const
 {
   double pKPi   = getPStar( mD, m, q );
   double mD2    = mD * mD;
   double m2     = m * m;
   double m02    = m0 * m0;
   double q2     = q * q;
   double mV2    = mV * mV;
   double mA2    = mA * mA;
   double summDm = mD + m;
   double TV     = TV_0 / ( 1.0 - q2 / ( mV2 ) );
   double T1     = T1_0 / ( 1.0 - q2 / ( mA2 ) );
   double T2     = T2_0 / ( 1.0 - q2 / ( mA2 ) );
 
   // construct amplitudes of (non)resonance
   double F     = getF2( m, m0, mPi, mK, rBW );
   double width = getWidth2( m, m0, mPi, mK, width0, rBW );
   EvtComplex C( m0 * width0 * F, 0.0 );
   double AA = m02 - m2;
   double BB = -m0 * width;
 
   EvtComplex amp = C / EvtComplex( AA, BB );
   amplitude      = abs( amp );
   delta          = atan2( imag( amp ), real( amp ) );
 
   F12 = amp * mD * pKPi / 3. *
         ( ( mD2 - m2 - q2 ) * summDm * T1 - mD2 * pKPi * pKPi / summDm * T2 );
   F22 = amp * root2d3 * mD * m * q * pKPi * summDm * T1;
   F32 = amp * root2d3 * 2. * mD2 * m * q * pKPi * pKPi / summDm * TV;
 }
 
 double EvtDToKpienu::getPStar( const double m, const double m1,
-                               const double m2 )
+                               const double m2 ) const
 {
   double s  = m * m;
   double s1 = m1 * m1;
   double s2 = m2 * m2;
   double x  = s + s1 - s2;
   double t  = 0.25 * x * x / s - s1;
   double p;
   if ( t > 0.0 ) {
     p = sqrt( t );
   } else {
-    std::cout << " Hello, pstar is less than 0.0" << std::endl;
+    EvtGenReport(EVTGEN_ERROR,"EvtGen") << " Hello, pstar is less than 0.0" << std::endl;
     p = 0.04;
   }
   return p;
 }
 
 double EvtDToKpienu::getF1( const double m, const double m0, const double m_c1, const double m_c2,
-                            const double rBW )
+                            const double rBW ) const
 {
   double pStar   = getPStar( m, m_c1, m_c2 );
   double pStar0  = getPStar( m0, m_c1, m_c2 );
   double rBW2    = rBW * rBW;
   double pStar2  = pStar * pStar;
   double pStar02 = pStar0 * pStar0;
   double B       = 1. / sqrt( 1. + rBW2 * pStar2 );
   double B0      = 1. / sqrt( 1. + rBW2 * pStar02 );
   double F       = pStar / pStar0 * B / B0;
   return F;
 }
 
 double EvtDToKpienu::getF2( const double m, const double m0, const double m_c1,
-                            const double m_c2, const double rBW )
+                            const double m_c2, const double rBW ) const
 {
   double pStar   = getPStar( m, m_c1, m_c2 );
   double pStar0  = getPStar( m0, m_c1, m_c2 );
   double rBW2    = rBW * rBW;
   double pStar2  = pStar * pStar;
   double pStar02 = pStar0 * pStar0;
   double B       = 1. / sqrt( ( rBW2 * pStar2 - 3. ) * ( rBW2 * pStar2 - 3. ) +
                         9. * rBW2 * pStar2 );
   double B0 = 1. / sqrt( ( rBW2 * pStar02 - 3. ) * ( rBW2 * pStar02 - 3. ) +
                          9. * rBW2 * pStar02 );
   double F  = pStar2 / pStar02 * B / B0;
   return F;
 }
 
 double EvtDToKpienu::getWidth0( const double m, const double m0,
                                 const double m_c1, const double m_c2,
-                                const double width0 )
+                                const double width0 ) const
 {
   double pStar  = getPStar( m, m_c1, m_c2 );
   double pStar0 = getPStar( m0, m_c1, m_c2 );
   double width  = width0 * pStar / pStar0 * m0 / m;
   return width;
 }
 
 double EvtDToKpienu::getWidth1( const double m, const double m0,
                                 const double m_c1, const double m_c2,
-                                const double width0, const double rBW )
+                                const double width0, const double rBW ) const
 {
   double pStar  = getPStar( m, m_c1, m_c2 );
   double pStar0 = getPStar( m0, m_c1, m_c2 );
   double F      = getF1( m, m0, m_c1, m_c2, rBW );
   double width  = width0 * pStar / pStar0 * m0 / m * F * F;
   return width;
 }
 
 double EvtDToKpienu::getWidth2( const double m, const double m0,
                                 const double m_c1, const double m_c2,
-                                const double width0, const double rBW )
+                                const double width0, const double rBW ) const
 {
   double pStar  = getPStar( m, m_c1, m_c2 );
   double pStar0 = getPStar( m0, m_c1, m_c2 );
   double F      = getF2( m, m0, m_c1, m_c2, rBW );
   double width  = width0 * pStar / pStar0 * m0 / m * F * F;
   return width;
 }
 
-EvtComplex EvtDToKpienu::getCoef( const double rho, const double phi )
+EvtComplex EvtDToKpienu::getCoef( const double rho, const double phi ) const
 {
   EvtComplex coef( rho * cos( phi ), rho * sin( phi ) );
   return coef;
 }