QCMCFilter.cxx

Go to the documentation of this file.

// QCMCFilter-00-00-02 Jan. 2013 Dan Ambrose
// Based on QCMCReweightProc program by Werner Sun
// This version has a corrected Y input from original version
#include "QCMCFilter.h"
#include "D0To2pip2pim.h"
#include "D0ToKSpipi.h"
#include "D0ToKSpipipi0.h"
#include "D0Topipipi0.h"
#include "D0Topippim2pi0.h"
#include "Dalitz.h"
// #include "D0ToKSKK.h"
#include "D0ToKLpipi.h"
#include "D0ToKSLKK.h"
// #include "D0ToKLpipi2018.h"
// #include "D0ToKSpipi2018.h"
#include "D0ToKKpipi.h"
#include "D0ToKpipipiLHCb.h"
#include "D0TopiKpipi.h"
// #include "models/dcs_LHCb_BESIIIIMP.h"
// #include "models/cf_LHCb_BESIIIIMP.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/RegistryEntry.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "TMath.h"
 
#include "HepPDT/ParticleData.hh"
#include "HepPDT/ParticleDataTable.hh"
#include <cmath>
// #include "PartPropSvc/PartPropSvc.h"
#include "GaudiKernel/IPartPropSvc.h"
 
// Monte-Carlo data
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
 
///////////////////
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
// #include "EvtGenBase/EvtComplex.hh"
// #include "EvtGenBase/EvtVector4R.hh"
 
#include "EvtRecEvent/EvtRecEvent.h"
 
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/TwoVector.h"
using CLHEP::Hep2Vector;
using CLHEP::Hep3Vector;
using CLHEP::HepLorentzVector;
using CLHEP::HepVector;
 
#include <complex>
#include <vector>
 
//
// constants, enums and typedefs
//
// Partilcle masses
const double xmpi0  = 0.134976; // BOSS 6.4.1 pdt.table value
const double xmeta  = 0.54784;  // PDG08 = 0.547853, BOSS 6.4.1 pdt.table = 0.54784
const double xmkaon = 0.49368;
const double xmpion = 0.13957;
const double xmk0   = 0.49761;
const double xmrho  = 0.77549;
const double xmd0   = 1.86484; // PDG08
 
const double PI = 3.1415926; // pi
 
// Antiparticles have negative ID.
 
static const int kPsi3770ID      = 30443;
static const int kD0ID           = 421;
static const int kD0BarID        = -421;
static const int kDpID           = 411;
static const int kDmID           = -411;
static const int kGammaID        = 22;
static const int kGammaFSRID     = -22;
static const int kPiPlusID       = 211;
static const int kPiMinusID      = -211;
static const int kPi0ID          = 111;
static const int kEtaID          = 221;
static const int kEtaPrimeID     = 331;
static const int kF0600ID        = 9000221;
static const int kF0ID           = 9010221;
static const int kFPrime0ID      = 10221;
static const int kF0m1500ID      = 50221;
static const int kF0m1710ID      = 10331;
static const int kF2ID           = 225;
static const int kA00ID          = 10111;
static const int kA0PlusID       = 10211;
static const int kA0MinusID      = -10211;
static const int kRhoPlusID      = 213;
static const int kRhoMinusID     = -213;
static const int kRho0ID         = 113;
static const int kRho2SPlusID    = 30213;
static const int kRho2SMinusID   = -30213;
static const int kRho2S0ID       = 30113;
static const int kA1PlusID       = 20213;
static const int kA1MinusID      = -20213;
static const int kA10ID          = 20113;
static const int kOmegaID        = 223;
static const int kPhiID          = 333;
static const int kKPlusID        = 321;
static const int kKMinusID       = -321;
static const int kK0SID          = 310;
static const int kK0LID          = 130;
static const int kK0ID           = 311;
static const int kK0BarID        = -311;
static const int kKStarPlusID    = 323;
static const int kKStarMinusID   = -323;
static const int kKStar0ID       = 313;
static const int kKStar0BarID    = -313;
static const int kKDPStar0ID     = 30313;
static const int kKDPStar0BarID  = -30313;
static const int kK0Star0ID      = 10311;
static const int kK0Star0BarID   = -10311;
static const int kK0StarPlusID   = 10321;
static const int kK0StarMinusID  = -10321;
static const int kK1PlusID       = 10323;
static const int kK1MinusID      = -10323;
static const int kK10ID          = 10313;
static const int kK10BarID       = -10313;
static const int kK1PrimePlusID  = 20323;
static const int kK1PrimeMinusID = -20323;
static const int kK1Prime0ID     = 20313;
static const int kK1Prime0BarID  = -20313;
static const int kK2StarPlusID   = 325;
static const int kK2StarMinusID  = -325;
static const int kK2Star0ID      = 315;
static const int kK2Star0BarID   = -315;
static const int kEMinusID       = 11;
static const int kEPlusID        = -11;
static const int kMuMinusID      = 13;
static const int kMuPlusID       = -13;
static const int kNuEID          = 12;
static const int kNuEBarID       = -12;
static const int kNuMuID         = 14;
static const int kNuMuBarID      = -14;
 
// D0 decay modes
static const int kFlavoredCF_0    = 0;
static const int kFlavoredCFBar_0 = 1;
static const int kFlavoredCF_1    = 2;
static const int kFlavoredCFBar_1 = 3;
static const int kFlavoredCF_3    = 4;
static const int kFlavoredCFBar_3 = 5;
static const int kFlavoredCS      = 6;
static const int kFlavoredCSBar   = 7;
static const int kSLPlus          = 8;
static const int kSLMinus         = 9;
static const int kCPPlus          = 10;
static const int kCPMinus         = 11;
static const int kDalitz          = 12;
static const int k2Pip2Pim        = 13;
static const int kPipPim2Pi0      = 14;
static const int kK0PiPiPi0       = 15;
static const int kKKPiPi          = 16;
static const int kK0KK            = 17;
static const int kPipPimPi0       = 18;
static const int kNDecayTypes     = 19;
 
// Varibales to keep track of the Events
int          m_nUnknownEvents     = 0;
int          m_nUnknownDecays     = 0;
int          m_nD0D0barEvents     = 0;
int          m_nD0bar             = 0;
int          m_nDpDmEvents        = 0;
int          m_nD0D0barDiscarded  = 0;
int          m_nDpDmDiscarded     = 0;
int          m_nCPPlus            = 0;
int          m_nCPMinus           = 0;
int          m_nFlavoredCFD0_0    = 0;
int          m_nFlavoredCFD0_1    = 0;
int          m_nFlavoredCFD0_3    = 0;
int          m_nFlavoredCSD0      = 0;
int          m_nFlavoredDCSD0_0   = 0;
int          m_nFlavoredDCSD0_1   = 0;
int          m_nFlavoredDCSD0_3   = 0;
int          m_nSL                = 0;
int          m_nDalitz            = 0;
int          m_n2Pip2Pim          = 0;
int          m_nPipPim2Pi0        = 0;
int          m_nK0PiPiPi0         = 0;
int          m_nKKPiPi            = 0;
int          m_nK0KK              = 0;
int          m_nPiPiPi0           = 0;
double       m_dalitzNumer1       = 0;
double       m_dalitzNumer2       = 0;
double       m_dalitzDenom        = 0;
double       dalitzNumer1_fil     = 0;
double       dalitzNumer2_fil     = 0;
double       dalitzDenom_fil      = 0;
double       m_2Pip2PimNumer1     = 0;
double       m_2Pip2PimNumer2     = 0;
double       m_2Pip2PimDenom      = 0;
double       Pip2Pim2Numer1_fil   = 0;
double       Pip2Pim2Numer2_fil   = 0;
double       Pip2Pim2Denom_fil    = 0;
double       m_PipPim2Pi0Numer1   = 0;
double       m_PipPim2Pi0Numer2   = 0;
double       m_PipPim2Pi0Denom    = 0;
double       PipPim2Pi0Numer1_fil = 0;
double       PipPim2Pi0Numer2_fil = 0;
double       PipPim2Pi0Denom_fil  = 0;
double       m_KSPiPiPi0Numer1    = 0;
double       m_KSPiPiPi0Numer2    = 0;
double       m_KSPiPiPi0Denom     = 0;
double       KSPiPiPi0Numer1_fil  = 0;
double       KSPiPiPi0Numer2_fil  = 0;
double       KSPiPiPi0Denom_fil   = 0;
HepSymMatrix m_weights( 20, 0 );
double       m_rwsCF     = 0.;
double       m_rwsCF_0   = 0.;
double       m_rwsCF_1   = 0.;
double       m_rwsCF_3   = 0.;
double       m_rwsCS     = 0.;
double       m_deltaCF   = 0.;
double       m_deltaCF_0 = 0.;
double       m_deltaCF_1 = 0.;
double       m_deltaCF_3 = 0.;
double       m_deltaCS   = 0.;
HepMatrix    m_modeCounter( 21, 21, 0 );
HepMatrix    m_keptModeCounter( 21, 21, 0 );
 
int m_rho_flag;
 
/////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( QCMCFilter )
 
QCMCFilter::QCMCFilter( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  // Declare the properties
  declareProperty( "Debug", m_debug = true );
  declareProperty( "x", m_x = 0.00407 );
  declareProperty( "y", m_y = 0.00647 ); // For values outside of (-.11 , .06) measured y value
                                         // may vary from input y
 
  declareProperty( "rForCFModes",
                   m_rCF = 0.0621 ); // should be kept near this. based on MC for d0->kpi and
                                     // d0bar->kpi 0.0621 = sqrt( .00015/.0389)
  declareProperty( "zForCFModes", m_zCF = 2.0 );
  declareProperty( "wSignForCFModes", m_wCFSign = true );
  declareProperty( "rForCFMode_0",
                   m_rCF_0 = 0.0587 ); // should be kept near this. based on MC for d0->kpi and
                                       // d0bar->kpi 0.0621 = sqrt( .00015/.0389)
  declareProperty( "RForCFMode_0", m_RCF_0 = 1. ); // should be kept near this.
  declareProperty( "zForCFMode_0", m_zCF_0 = 1.9585 );
  declareProperty( "dForCFMode_0", m_dCF_0 = 191.7 * 3.14159 / 180. );
  declareProperty( "wSignForCFMode_0", m_wCFSign_0 = true );
  declareProperty( "rForCFModes_1",
                   m_rCF_1 = 0.044 ); // should be kept near this. based on MC for d0->kpi and
                                      // d0bar->kpi 0.0621 = sqrt( .00015/.0389)
  declareProperty( "RForCFMode_1", m_RCF_1 = 0.79 ); // should be kept near this.
  declareProperty( "zForCFModes_1", m_zCF_1 = 1.9225 );
  declareProperty( "dForCFMode_1", m_dCF_1 = 196. * 3.14159 / 180. );
  declareProperty( "wSignForCFModes_1", m_wCFSign_1 = true );
  declareProperty( "rForCFModes_3",
                   m_rCF_3 = 0.0549 ); // should be kept near this. based on MC for d0->kpi and
                                       // d0bar->kpi 0.0621 = sqrt( .00015/.0389)
  declareProperty( "RForCFMode_3", m_RCF_3 = 0.43 ); // should be kept near this.
  declareProperty( "zForCFModes_3", m_zCF_3 = 1.2313 );
  declareProperty( "dForCFMode_3", m_dCF_3 = 161. * 3.14159 / 180. );
  declareProperty( "wSignForCFModes_3", m_wCFSign_3 = true );
 
  declareProperty( "FCP_PiPiPi0", m_FCP_PiPiPi0 = 0.973 );
 
  declareProperty( "rForCSModes", m_rCS = 1.0 );
  declareProperty( "zForCSModes", m_zCS = -2.0 );
  declareProperty( "wSignForCSModes", m_wCSSign = true );
  declareProperty( "DplusDminusWeight", m_dplusDminusWeight = 0.5 );
  declareProperty( "dalitzModel", m_dalitzModel = 0 ); // 0 CLEO-c, 1 Babar, 2 Belle
  declareProperty( "UseNewWeights", m_useNewWeights = true );
  declareProperty( "InvertSelection", m_invertSelection = false );
}
 
//***************************************************************
 
StatusCode QCMCFilter::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  // Get the Particle Properties Service
  IPartPropSvc*     p_PartPropSvc;
  static const bool CREATEIFNOTTHERE( true );
  StatusCode        PartPropStatus = service( "PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE );
  if ( !PartPropStatus.isSuccess() || 0 == p_PartPropSvc )
  {
    log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
    return PartPropStatus;
  }
  m_particleTable = p_PartPropSvc->PDT();
 
  // Calculates parameters here based on Parameter Input by User
  //  (e.g. run expensive algorithms that are based on parameters
  //   specified by user at run-time)
 
  // m_y = (m_y - 0.008)/0.292;
  // corrects the y input so that we get the desired measurable y from CP:Semi-leptonic
  // doubletag method
  //  This equation is a result of solving equations for y dependence.
  // If parent MC Branching ratio's this equation will need resolved. January 2013 -DA
  double x     = m_x;
  double y     = m_y;
  double rCF_0 = m_rCF_0;
  double dCF_0 = m_dCF_0;
  double RCF_0 = m_RCF_0;
  // double zCF_0 = m_zCF_0 ;
  double zCF_0  = 2 * cos( dCF_0 );
  double rCF2_0 = rCF_0 * rCF_0;
  double zCF2_0 = zCF_0 * zCF_0;
  double rCF_1  = m_rCF_1;
  double dCF_1  = m_dCF_1;
  double RCF_1  = m_RCF_1;
  // double zCF_1 = m_zCF_1 ;
  double zCF_1  = 2 * cos( dCF_1 );
  double rCF2_1 = rCF_1 * rCF_1;
  double zCF2_1 = zCF_1 * zCF_1;
  double rCF_3  = m_rCF_3;
  double dCF_3  = m_dCF_3;
  double RCF_3  = m_RCF_3;
  // double zCF_3 = m_zCF_3 ;
  double zCF_3  = 2 * cos( dCF_3 );
  double rCF2_3 = rCF_3 * rCF_3;
  double zCF2_3 = zCF_3 * zCF_3;
  double rCS    = m_rCS;
  double zCS    = m_zCS;
  double rCS2   = rCS * rCS;
  double zCS2   = zCS * zCS;
 
  double rmix         = ( x * x + y * y ) / 2.;
  double wCF_0        = 2 * cos( dCF_0 );
  double wCF_1        = 2 * cos( dCF_1 );
  double wCF_3        = 2 * cos( dCF_3 );
  double wCS          = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 );
  double vCF_0CSPlus  = ( zCF_0 * zCS + wCF_0 * wCS ) / 2.;
  double vCF_1CSPlus  = ( zCF_1 * zCS + wCF_1 * wCS ) / 2.;
  double vCF_3CSPlus  = ( zCF_3 * zCS + wCF_3 * wCS ) / 2.;
  double vCF_0CSMinus = ( zCF_0 * zCS - wCF_0 * wCS ) / 2.;
  double vCF_1CSMinus = ( zCF_1 * zCS - wCF_1 * wCS ) / 2.;
  double vCF_3CSMinus = ( zCF_3 * zCS - wCF_3 * wCS ) / 2.;
 
  m_deltaCF_0 = acos( zCF_0 / 2. ) * ( m_wCFSign_0 ? 1. : -1. );
  m_deltaCF_1 = acos( zCF_1 / 2. ) * ( m_wCFSign_1 ? 1. : -1. );
  m_deltaCF_3 = acos( zCF_3 / 2. ) * ( m_wCFSign_3 ? 1. : -1. );
  m_deltaCS   = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. );
 
  double rCF        = m_rCF;
  double zCF        = m_zCF;
  double rCF2       = rCF * rCF;
  double zCF2       = zCF * zCF;
  double wCF        = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 );
  double vCFCSPlus  = ( zCF * zCS + wCF * wCS ) / 2.;
  double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2.;
  m_deltaCF         = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. );
 
  double bCF_0 = 1. + 0.5 * rCF_0 * ( zCF_0 * y + wCF_0 * x ) +
                 0.5 * ( y * y - x * x ); // second order y term added -DA
  double bCF_1 = 1. + 0.5 * rCF_1 * ( zCF_1 * y + wCF_1 * x ) +
                 0.5 * ( y * y - x * x ); // second order y term added -DA
  double bCF_3 = 1. + 0.5 * rCF_3 * ( zCF_3 * y + wCF_3 * x ) +
                 0.5 * ( y * y - x * x ); // second order y term added -DA
  double bCS      = 1. + 0.5 * rCS * ( zCS * y + wCS * x ) + 0.5 * ( y * y - x * x );
  double bCPPlus  = 1. - y;
  double bCPMinus = 1. + y;
 
  double bCPPlus_PiPiPi0 = 1. - ( 2. * m_FCP_PiPiPi0 - 1. ) * y;
 
  if ( !m_useNewWeights )
  {
    bCF_0     = 1.;
    bCF_1     = 1.;
    bCF_3     = 1.;
    m_rwsCF_0 = rCF2_0;
    m_rwsCF_1 = rCF2_1;
    m_rwsCF_3 = rCF2_3;
    bCS       = 1.;
    m_rwsCS   = rCS2;
    bCPPlus   = 1.;
    bCPMinus  = 1.;
  }
 
  m_rwsCF_0 = ( rCF2_0 + 0.5 * rCF_0 * ( zCF_0 * y - wCF_0 * x ) + rmix ) /
              bCF_0; // though division by b not mentioned in memo it is correct
  m_rwsCF_1 = ( rCF2_1 + 0.5 * rCF_1 * ( zCF_1 * y - wCF_1 * x ) + rmix ) /
              bCF_1; // though division by b not mentioned in memo it is correct
  m_rwsCF_3 = ( rCF2_3 + 0.5 * rCF_3 * ( zCF_3 * y - wCF_3 * x ) + rmix ) /
              bCF_3; // though division by b not mentioned in memo it is correct
  m_rwsCS = ( rCS2 + 0.5 * rCS * ( zCS * y - wCS * x ) + rmix ) / bCS;
 
  double bfCF_0    = ( 1. - RCF_0 * rCF_0 * ( 0.5 * zCF_0 * y + 0.5 * wCF_0 * x ) +
                    0.5 * ( -1 * x * x + y * y ) ); // B(D0->F)=AF*AF*brCF_0
  double bfCF_1    = ( 1. - RCF_1 * rCF_1 * ( 0.5 * zCF_1 * y + 0.5 * wCF_1 * x ) +
                    0.5 * ( -1 * x * x + y * y ) );
  double bfCF_3    = ( 1. - RCF_3 * rCF_3 * ( 0.5 * zCF_3 * y + 0.5 * wCF_3 * x ) +
                    0.5 * ( -1 * x * x + y * y ) );
  double bfCFbar_0 = ( rCF2_0 - RCF_0 * rCF_0 * ( 0.5 * zCF_0 * y - 0.5 * wCF_0 * x ) +
                       0.5 * ( x * x + y * y ) ); // B(D0->Fbar)=AF*AF*brCFbar_0
  double bfCFbar_1 = ( rCF2_1 - RCF_1 * rCF_1 * ( 0.5 * zCF_1 * y - 0.5 * wCF_1 * x ) +
                       0.5 * ( x * x + y * y ) );
  double bfCFbar_3 = ( rCF2_3 - RCF_3 * rCF_3 * ( 0.5 * zCF_3 * y - 0.5 * wCF_3 * x ) +
                       0.5 * ( x * x + y * y ) );
 
  // Calculate reweighting factors, normalized to the largest factor.
  // Note:A fake Dalitz weight is added in to just to make initial predictions
  // These Dalitz weights are not used in actual code
 
  // FlavoredCF
 
  // First calculate rate factors.
  //    m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] = ( ( 1 -  RCF_0 * RCF_0 ) / ( 1 -  RCF_0 *
  //(1./rCF_0) * (0.5*zCF_0*y - 0.5*wCF_0*x) + rmix/(2*rCF2_0) ) ) ;        rmix *
  //m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;    m_weights[ kFlavoredCF_0 ][
  // kFlavoredCFBar_0 ] = ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( brCF_0 *
  // brCF_0 );
  //        1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
 
  // Nov 2007: correct for r2 -> RWS and BR != A^2
  //    m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
  //    m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
  //        ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
  //    m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
  //        ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
  //    m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
  //        ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
  m_weights[kFlavoredCF_0][kFlavoredCF_0] =
      ( ( 1 - RCF_0 * RCF_0 ) /
        ( 1 - RCF_0 * ( 1. / rCF_0 ) * ( 0.5 * zCF_0 * y - 0.5 * wCF_0 * x ) +
          rmix / ( 2 * rCF2_0 ) ) );
  m_weights[kFlavoredCF_0][kFlavoredCFBar_0] =
      ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( bfCF_0 * bfCF_0 );
  m_weights[kFlavoredCF_0][kFlavoredCF_1] =
      ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) -
        2. * ( rCF2_0 / rCF2_1 ) * RCF_0 * RCF_1 * cos( dCF_0 - dCF_1 ) ) /
      ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) -
        RCF_1 * ( 1. / rCF_1 ) * ( 0.5 * zCF_1 * y - 0.5 * wCF_1 * x ) -
        ( RCF_0 * rCF_1 / ( rCF_0 * rCF_0 ) ) * ( 0.5 * zCF_0 * y - 0.5 * wCF_0 * x ) +
        rmix / ( rCF_0 * rCF_0 ) );
  cout << ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) -
            RCF_1 * ( 1. / rCF_1 ) * ( 0.5 * zCF_1 * y - 0.5 * wCF_1 * x ) -
            ( RCF_0 * rCF_1 / ( rCF_0 * rCF_0 ) ) * ( 0.5 * zCF_0 * y - 0.5 * wCF_0 * x ) +
            rmix / ( rCF_0 * rCF_0 ) )
       << endl;
  m_weights[kFlavoredCF_0][kFlavoredCFBar_1] =
      ( 1. + rCF2_0 * rCF2_1 - 2. * RCF_0 * RCF_1 * rCF_0 * rCF_1 ) / ( bfCF_0 * bfCF_1 );
  m_weights[kFlavoredCF_0][kFlavoredCF_3] =
      ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) -
        2. * ( rCF2_0 / rCF2_3 ) * RCF_0 * RCF_3 * cos( dCF_0 - dCF_3 ) ) /
      ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) -
        RCF_3 * ( 1. / rCF_3 ) * ( 0.5 * zCF_3 * y - 0.5 * wCF_3 * x ) -
        ( RCF_0 / ( rCF_0 * rCF_0 ) ) * ( 0.5 * zCF_0 * y - 0.5 * wCF_0 * x ) +
        rmix / ( rCF_0 * rCF_0 ) );
  m_weights[kFlavoredCF_0][kFlavoredCFBar_3] =
      ( 1. + rCF2_0 * rCF2_3 - 2. * RCF_0 * RCF_3 * rCF_0 * rCF_3 ) / ( bfCF_0 * bfCF_3 );
  m_weights[kFlavoredCF_0][kFlavoredCSBar] = 1. + rCF2_0 * rCS2 - rCF_0 * rCS * vCFCSMinus;
  m_weights[kFlavoredCF_0][kFlavoredCS]    = rCF2_0 + rCS2 - rCF_0 * rCS * vCF_0CSPlus +
                                          rmix * m_weights[kFlavoredCF_0][kFlavoredCSBar];
 
  m_weights[kFlavoredCF_0][kSLPlus]  = 1. / bfCF_0; // inexist
  m_weights[kFlavoredCF_0][kSLMinus] = 1. / bfCF_0; //
 
  // m_weights[ kFlavoredCF_0 ][ kCPPlus ] =
  //    ( 1. + rCF2_0 + rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPPlus ;
  // m_weights[ kFlavoredCF_0 ][ kCPMinus ] =
  //    ( 1. + rCF2_0 - rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPMinus ;
  m_weights[kFlavoredCF_0][kCPPlus] =
      ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos( dCF_0 ) ) /
      ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos( dCF_0 ) * y + y * y ) * bCPPlus );
  m_weights[kFlavoredCF_0][kCPMinus] =
      ( 1. + rCF2_0 + 2. * rCF_0 * RCF_0 * cos( dCF_0 ) ) /
      ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos( dCF_0 ) * y + y * y ) * bCPMinus );
  m_weights[kFlavoredCF_0][kDalitz]     = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][k2Pip2Pim]   = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][kPipPim2Pi0] = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][kK0PiPiPi0]  = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][kKKPiPi]     = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][kK0KK]       = 1. / bfCF_0;
  m_weights[kFlavoredCF_0][kPipPimPi0] =
      ( 1. + rCF2_0 - 2. * ( 2. * m_FCP_PiPiPi0 - 1. ) * rCF_0 * RCF_0 * cos( dCF_0 ) ) /
      ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos( dCF_0 ) * y + y * y ) * bCPPlus_PiPiPi0 );
 
  // FlavoredCFBar
  m_weights[kFlavoredCFBar_0][kFlavoredCFBar_0] = m_weights[kFlavoredCF_0][kFlavoredCF_0];
  m_weights[kFlavoredCFBar_0][kFlavoredCF_1]    = m_weights[kFlavoredCF_0][kFlavoredCFBar_1];
  m_weights[kFlavoredCFBar_0][kFlavoredCFBar_1] = m_weights[kFlavoredCF_0][kFlavoredCF_1];
  m_weights[kFlavoredCFBar_0][kFlavoredCF_3]    = m_weights[kFlavoredCF_0][kFlavoredCFBar_3];
  m_weights[kFlavoredCFBar_0][kFlavoredCFBar_3] = m_weights[kFlavoredCF_0][kFlavoredCF_3];
  m_weights[kFlavoredCFBar_0][kFlavoredCS]      = m_weights[kFlavoredCF_0][kFlavoredCSBar];
  m_weights[kFlavoredCFBar_0][kFlavoredCSBar]   = m_weights[kFlavoredCF_0][kFlavoredCS];
  m_weights[kFlavoredCFBar_0][kSLPlus]          = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kSLMinus]         = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kCPPlus]          = m_weights[kFlavoredCF_0][kCPPlus];
  m_weights[kFlavoredCFBar_0][kCPMinus]         = m_weights[kFlavoredCF_0][kCPMinus];
  m_weights[kFlavoredCFBar_0][kDalitz]          = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][k2Pip2Pim]        = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kPipPim2Pi0]      = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kK0PiPiPi0]       = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kKKPiPi]          = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kK0KK]            = 1. / bCF_0;
  m_weights[kFlavoredCFBar_0][kPipPimPi0]       = m_weights[kFlavoredCF_0][kPipPimPi0];
 
  m_weights[kFlavoredCF_1][kFlavoredCF_1] =
      ( ( 1 - RCF_1 * RCF_1 ) /
        ( 1 - RCF_1 * ( 1. / rCF_1 ) * ( 0.5 * zCF_1 * y - 0.5 * wCF_1 * x ) +
          rmix / ( 2 * rCF2_1 ) ) );
  m_weights[kFlavoredCF_1][kFlavoredCFBar_1] =
      ( 1. + rCF2_1 * rCF2_1 - 2. * RCF_1 * RCF_1 * rCF2_1 ) / ( bfCF_1 * bfCF_1 );
  m_weights[kFlavoredCF_1][kFlavoredCF_3] =
      ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) -
        2. * ( rCF2_1 / rCF2_3 ) * RCF_1 * RCF_3 * cos( dCF_1 - dCF_3 ) ) /
      ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) -
        RCF_3 * ( 1. / rCF_3 ) * ( 0.5 * zCF_3 * y - 0.5 * wCF_3 * x ) -
        ( RCF_1 / ( rCF_1 * rCF_1 ) ) * ( 0.5 * zCF_1 * y - 0.5 * wCF_1 * x ) +
        rmix / ( rCF_1 * rCF_1 ) );
  ;
  m_weights[kFlavoredCF_1][kFlavoredCFBar_3] =
      ( 1. + rCF2_1 * rCF2_3 - 2. * RCF_1 * RCF_3 * rCF_1 * rCF_3 ) / ( bfCF_1 * bfCF_3 );
  ;
  m_weights[kFlavoredCF_1][kFlavoredCSBar] = 1. + rCF2_1 * rCS2 - rCF_1 * rCS * vCFCSMinus;
  m_weights[kFlavoredCF_1][kFlavoredCS]    = rCF2_1 + rCS2 - rCF_1 * rCS * vCF_1CSPlus +
                                          rmix * m_weights[kFlavoredCF_1][kFlavoredCSBar];
 
  m_weights[kFlavoredCF_1][kSLPlus]  = 1. / bfCF_1;
  m_weights[kFlavoredCF_1][kSLMinus] = 1. / bfCF_1;
 
  //    m_weights[ kFlavoredCF_1 ][ kCPPlus ] =
  //        ( 1. + rCF2_1 + rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPPlus ;
  //    m_weights[ kFlavoredCF_1 ][ kCPMinus ] =
  //        ( 1. + rCF2_1 - rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPMinus ;
  m_weights[kFlavoredCF_1][kCPPlus] =
      ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos( dCF_1 ) ) /
      ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos( dCF_1 ) * y + y * y ) * bCPPlus );
  m_weights[kFlavoredCF_1][kCPMinus] =
      ( 1. + rCF2_1 + 2. * rCF_1 * RCF_1 * cos( dCF_1 ) ) /
      ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos( dCF_1 ) * y + y * y ) * bCPMinus );
  m_weights[kFlavoredCF_1][kDalitz]     = 1. / bCF_1;
  m_weights[kFlavoredCF_1][k2Pip2Pim]   = 1. / bCF_1;
  m_weights[kFlavoredCF_1][kPipPim2Pi0] = 1. / bCF_1;
  m_weights[kFlavoredCF_1][kK0PiPiPi0]  = 1. / bCF_1;
  m_weights[kFlavoredCF_1][kKKPiPi]     = 1. / bCF_1;
  m_weights[kFlavoredCF_1][kK0KK]       = 1. / bCF_1;
  m_weights[kFlavoredCF_1][kPipPimPi0] =
      ( 1. + rCF2_1 - 2. * ( 2. * m_FCP_PiPiPi0 - 1. ) * rCF_1 * RCF_1 * cos( dCF_1 ) ) /
      ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos( dCF_1 ) * y + y * y ) * bCPPlus_PiPiPi0 );
 
  // FlavoredCFBar
  m_weights[kFlavoredCFBar_1][kFlavoredCFBar_1] = m_weights[kFlavoredCF_1][kFlavoredCF_1];
  m_weights[kFlavoredCFBar_1][kFlavoredCF_3]    = m_weights[kFlavoredCF_1][kFlavoredCFBar_3];
  m_weights[kFlavoredCFBar_1][kFlavoredCFBar_3] = m_weights[kFlavoredCF_1][kFlavoredCF_3];
  m_weights[kFlavoredCFBar_1][kFlavoredCS]      = m_weights[kFlavoredCF_1][kFlavoredCSBar];
  m_weights[kFlavoredCFBar_1][kFlavoredCSBar]   = m_weights[kFlavoredCF_1][kFlavoredCS];
  m_weights[kFlavoredCFBar_1][kSLPlus]          = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kSLMinus]         = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kCPPlus]          = m_weights[kFlavoredCF_1][kCPPlus];
  m_weights[kFlavoredCFBar_1][kCPMinus]         = m_weights[kFlavoredCF_1][kCPMinus];
  m_weights[kFlavoredCFBar_1][kDalitz]          = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][k2Pip2Pim]        = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kPipPim2Pi0]      = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kK0PiPiPi0]       = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kKKPiPi]          = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kK0KK]            = 1. / bCF_1;
  m_weights[kFlavoredCFBar_1][kPipPimPi0]       = m_weights[kFlavoredCF_1][kPipPimPi0];
 
  m_weights[kFlavoredCF_3][kFlavoredCF_3] =
      ( ( 1 - RCF_3 * RCF_3 ) /
        ( 1 - RCF_3 * ( 1. / rCF_3 ) * ( 0.5 * zCF_3 * y - 0.5 * wCF_3 * x ) +
          rmix / ( 2 * rCF2_3 ) ) );
  m_weights[kFlavoredCF_3][kFlavoredCFBar_3] =
      ( 1. + rCF2_3 * rCF2_3 - 2. * RCF_3 * RCF_3 * rCF2_3 ) / ( bfCF_3 * bfCF_3 );
  m_weights[kFlavoredCF_3][kFlavoredCSBar] = 1. + rCF2_3 * rCS2 - rCF_3 * rCS * vCFCSMinus;
  m_weights[kFlavoredCF_3][kFlavoredCS]    = rCF2_3 + rCS2 - rCF_3 * rCS * vCF_3CSPlus +
                                          rmix * m_weights[kFlavoredCF_3][kFlavoredCSBar];
 
  m_weights[kFlavoredCF_3][kSLPlus]  = 1. / bfCF_3;
  m_weights[kFlavoredCF_3][kSLMinus] = 1. / bfCF_3;
 
  // m_weights[ kFlavoredCF_3 ][ kCPPlus ] =
  //    ( 1. + rCF2_3 + rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPPlus ;
  // m_weights[ kFlavoredCF_3 ][ kCPMinus ] =
  //    ( 1. + rCF2_3 - rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPMinus ;
  m_weights[kFlavoredCF_3][kCPPlus] =
      ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos( dCF_3 ) ) /
      ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos( dCF_3 ) * y + y * y ) * bCPPlus );
  m_weights[kFlavoredCF_3][kCPMinus] =
      ( 1. + rCF2_3 + 2. * rCF_3 * RCF_3 * cos( dCF_3 ) ) /
      ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos( dCF_3 ) * y + y * y ) * bCPMinus );
  m_weights[kFlavoredCF_3][kDalitz]     = 1. / bCF_3;
  m_weights[kFlavoredCF_3][k2Pip2Pim]   = 1. / bCF_3;
  m_weights[kFlavoredCF_3][kPipPim2Pi0] = 1. / bCF_3;
  m_weights[kFlavoredCF_3][kK0PiPiPi0]  = 1. / bCF_3;
  m_weights[kFlavoredCF_3][kKKPiPi]     = 1. / bCF_3;
  m_weights[kFlavoredCF_3][kK0KK]       = 1. / bCF_3;
  m_weights[kFlavoredCF_3][kPipPimPi0] =
      ( 1. + rCF2_3 - 2. * ( 2. * m_FCP_PiPiPi0 - 1. ) * rCF_3 * RCF_3 * cos( dCF_3 ) ) /
      ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos( dCF_3 ) * y + y * y ) * bCPPlus_PiPiPi0 );
 
  // FlavoredCFBar
  m_weights[kFlavoredCFBar_3][kFlavoredCFBar_3] = m_weights[kFlavoredCF_3][kFlavoredCF_3];
  m_weights[kFlavoredCFBar_3][kFlavoredCS]      = m_weights[kFlavoredCF_3][kFlavoredCSBar];
  m_weights[kFlavoredCFBar_3][kFlavoredCSBar]   = m_weights[kFlavoredCF_3][kFlavoredCS];
  m_weights[kFlavoredCFBar_3][kSLPlus]          = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kSLMinus]         = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kCPPlus]          = m_weights[kFlavoredCF_3][kCPPlus];
  m_weights[kFlavoredCFBar_3][kCPMinus]         = m_weights[kFlavoredCF_3][kCPMinus];
  m_weights[kFlavoredCFBar_3][kDalitz]          = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][k2Pip2Pim]        = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kPipPim2Pi0]      = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kK0PiPiPi0]       = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kKKPiPi]          = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kK0KK]            = 1. / bCF_3;
  m_weights[kFlavoredCFBar_3][kPipPimPi0]       = m_weights[kFlavoredCF_3][kPipPimPi0];
 
  // FlavoredCS
 
  // First calculate rate factors.
  m_weights[kFlavoredCS][kFlavoredCSBar] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2;
  m_weights[kFlavoredCS][kFlavoredCS]    = rmix * m_weights[kFlavoredCS][kFlavoredCSBar];
 
  // Nov 2007: correct for r2 -> RWS and BR != A^2
  m_weights[kFlavoredCS][kFlavoredCS] /= m_rwsCS * bCS * bCS;
  m_weights[kFlavoredCS][kFlavoredCSBar] /= ( 1. + m_rwsCS * m_rwsCS ) * bCS * bCS;
 
  m_weights[kFlavoredCS][kSLPlus]  = 1. / bCS;
  m_weights[kFlavoredCS][kSLMinus] = 1. / bCS;
 
  m_weights[kFlavoredCS][kCPPlus] =
      ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus;
  m_weights[kFlavoredCS][kCPMinus] =
      ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus;
 
  m_weights[kFlavoredCS][kDalitz]     = 1. / bCS;
  m_weights[kFlavoredCS][k2Pip2Pim]   = 1. / bCS;
  m_weights[kFlavoredCS][kPipPim2Pi0] = 1. / bCS;
  m_weights[kFlavoredCS][kK0PiPiPi0]  = 1. / bCS;
  m_weights[kFlavoredCS][kKKPiPi]     = 1. / bCS;
  m_weights[kFlavoredCS][kK0KK]       = 1. / bCS;
  m_weights[kFlavoredCS][kPipPimPi0] =
      ( 1. + rCS2 + ( 2. * m_FCP_PiPiPi0 - 1. ) * rCS * zCS ) / ( 1. + m_rwsCS ) / bCS /
      bCPPlus_PiPiPi0;
 
  // FlavoredCSBar
 
  m_weights[kFlavoredCSBar][kFlavoredCSBar] = m_weights[kFlavoredCS][kFlavoredCS];
  m_weights[kFlavoredCSBar][kSLPlus]        = 1. / bCS;
  m_weights[kFlavoredCSBar][kSLMinus]       = 1. / bCS;
  m_weights[kFlavoredCSBar][kCPPlus]        = m_weights[kFlavoredCS][kCPPlus];
  m_weights[kFlavoredCSBar][kCPMinus]       = m_weights[kFlavoredCS][kCPMinus];
  m_weights[kFlavoredCSBar][kDalitz]        = 1. / bCS;
  m_weights[kFlavoredCSBar][k2Pip2Pim]      = 1. / bCS;
  m_weights[kFlavoredCSBar][kPipPim2Pi0]    = 1. / bCS;
  m_weights[kFlavoredCSBar][kK0PiPiPi0]     = 1. / bCS;
  m_weights[kFlavoredCSBar][kKKPiPi]        = 1. / bCS;
  m_weights[kFlavoredCSBar][kK0KK]          = 1. / bCS;
  m_weights[kFlavoredCSBar][kPipPimPi0]     = m_weights[kFlavoredCS][kPipPimPi0];
 
  // SLPlus
 
  // SLPlus/SLPlus should have rate factor = Rmix, but none of these are
  // generated in uncorrelated MC.
  m_weights[kSLPlus][kSLPlus]     = 0.;
  m_weights[kSLPlus][kSLMinus]    = 1.;
  m_weights[kSLPlus][kCPPlus]     = 1. / bCPPlus;
  m_weights[kSLPlus][kCPMinus]    = 1. / bCPMinus;
  m_weights[kSLPlus][kDalitz]     = 1.;
  m_weights[kSLPlus][k2Pip2Pim]   = 1.;
  m_weights[kSLPlus][kPipPim2Pi0] = 1.;
  m_weights[kSLPlus][kK0PiPiPi0]  = 1.;
  m_weights[kSLPlus][kKKPiPi]     = 1.;
  m_weights[kSLPlus][kK0KK]       = 1.;
  m_weights[kSLPlus][kPipPimPi0]  = 1. / bCPPlus_PiPiPi0;
 
  // SLMinus
 
  m_weights[kSLMinus][kSLMinus]    = 0.;
  m_weights[kSLMinus][kCPPlus]     = 1. / bCPPlus;
  m_weights[kSLMinus][kCPMinus]    = 1. / bCPMinus;
  m_weights[kSLMinus][kDalitz]     = 1.;
  m_weights[kSLMinus][k2Pip2Pim]   = 1.;
  m_weights[kSLMinus][kPipPim2Pi0] = 1.;
  m_weights[kSLMinus][kK0PiPiPi0]  = 1.;
  m_weights[kSLMinus][kKKPiPi]     = 1.;
  m_weights[kSLMinus][kK0KK]       = 1.;
  m_weights[kSLMinus][kPipPimPi0]  = 1. / bCPPlus_PiPiPi0;
 
  // CPPlus
 
  m_weights[kCPPlus][kCPPlus]     = 0.;
  m_weights[kCPPlus][kCPMinus]    = 2. / bCPPlus / bCPMinus;
  m_weights[kCPPlus][kDalitz]     = 1. / bCPPlus;
  m_weights[kCPPlus][k2Pip2Pim]   = 1. / bCPPlus;
  m_weights[kCPPlus][kPipPim2Pi0] = 1. / bCPPlus;
  m_weights[kCPPlus][kK0PiPiPi0]  = 1. / bCPPlus;
  m_weights[kCPPlus][kKKPiPi]     = 1. / bCPPlus;
  m_weights[kCPPlus][kK0KK]       = 1. / bCPPlus;
  m_weights[kCPPlus][kPipPimPi0] =
      ( 1 - ( 2. * m_FCP_PiPiPi0 - 1. ) ) / bCPPlus / bCPPlus_PiPiPi0;
 
  // CPMinus
 
  m_weights[kCPMinus][kCPMinus]    = 0.;
  m_weights[kCPMinus][kDalitz]     = 1. / bCPMinus;
  m_weights[kCPMinus][k2Pip2Pim]   = 1. / bCPMinus;
  m_weights[kCPMinus][kPipPim2Pi0] = 1. / bCPMinus;
  m_weights[kCPMinus][kK0PiPiPi0]  = 1. / bCPMinus;
  m_weights[kCPMinus][kKKPiPi]     = 1. / bCPMinus;
  m_weights[kCPMinus][kK0KK]       = 1. / bCPMinus;
  m_weights[kCPMinus][kPipPimPi0] =
      2. * ( 1 - ( 2. * m_FCP_PiPiPi0 - 1. ) ) / bCPMinus / bCPPlus_PiPiPi0;
 
  // Dalitz estimate
  m_weights[kDalitz][kDalitz]     = 1;
  m_weights[kDalitz][k2Pip2Pim]   = 1;
  m_weights[kDalitz][kPipPim2Pi0] = 1;
  m_weights[kDalitz][kK0PiPiPi0]  = 1;
  m_weights[kDalitz][kKKPiPi]     = 1;
  m_weights[kDalitz][kK0KK]       = 1;
  m_weights[kDalitz][kPipPimPi0]  = 1;
 
  m_weights[k2Pip2Pim][k2Pip2Pim]   = 1;
  m_weights[k2Pip2Pim][kPipPim2Pi0] = 1;
  m_weights[k2Pip2Pim][kK0PiPiPi0]  = 1;
  m_weights[k2Pip2Pim][kKKPiPi]     = 1;
  m_weights[k2Pip2Pim][kK0KK]       = 1;
  m_weights[k2Pip2Pim][kPipPimPi0]  = 1;
 
  m_weights[kPipPim2Pi0][kPipPim2Pi0] = 1;
  m_weights[kPipPim2Pi0][kK0PiPiPi0]  = 1;
  m_weights[kPipPim2Pi0][kKKPiPi]     = 1;
  m_weights[kPipPim2Pi0][kK0KK]       = 1;
  m_weights[kPipPim2Pi0][kPipPimPi0]  = 1;
 
  m_weights[kK0PiPiPi0][kK0PiPiPi0] = 1;
  m_weights[kK0PiPiPi0][kKKPiPi]    = 1;
  m_weights[kK0PiPiPi0][kK0KK]      = 1;
  m_weights[kK0PiPiPi0][kPipPimPi0] = 1;
 
  m_weights[kKKPiPi][kKKPiPi]    = 1;
  m_weights[kKKPiPi][kK0KK]      = 1;
  m_weights[kKKPiPi][kPipPimPi0] = 1;
 
  m_weights[kK0KK][kK0KK]      = 1;
  m_weights[kK0KK][kPipPimPi0] = 1;
 
  m_weights[kPipPimPi0][kPipPimPi0] =
      ( 1 - ( 2. * m_FCP_PiPiPi0 - 1. ) * ( 2. * m_FCP_PiPiPi0 - 1. ) ) / bCPPlus_PiPiPi0 /
      bCPPlus_PiPiPi0;
 
  log << MSG::INFO << "Weight matrix" << m_weights << endmsg;
 
  // Boss 6.6.2 MC Branching fractions
  HepVector fractionsD0( kNDecayTypes + 1, 0 );
  fractionsD0[kFlavoredCF_0]    = 0.305;
  fractionsD0[kFlavoredCFBar_0] = 0.0076;
  fractionsD0[kFlavoredCF_1]    = 0.144;
  fractionsD0[kFlavoredCFBar_1] = 0.00031;
  fractionsD0[kFlavoredCF_3]    = 0.0822;
  fractionsD0[kFlavoredCFBar_3] = 0.00027;
  fractionsD0[kFlavoredCS]      = 0.031;
  fractionsD0[kFlavoredCSBar]   = 0.031;
  fractionsD0[kSLPlus]          = 0.125;
  fractionsD0[kSLMinus]         = 0.;
  fractionsD0[kCPPlus]          = 0.113;
  fractionsD0[kCPMinus]         = 0.102;
  fractionsD0[kDalitz]          = 0.05855;
  fractionsD0[k2Pip2Pim]        = 0.00720;
  fractionsD0[kPipPim2Pi0]      = 0.0102;
  fractionsD0[kK0PiPiPi0]       = 0.105;
  fractionsD0[kKKPiPi]          = 0.00273;
  fractionsD0[kK0KK]            = 0.00902;
  fractionsD0[kPipPimPi0]       = 0.0149;
  HepVector scales              = m_weights * fractionsD0;
  log << MSG::INFO << "Integrated scale factors " << scales << endmsg;
 
  HepVector fractionsD0bar( kNDecayTypes + 1, 0 );
  fractionsD0bar[kFlavoredCF_0]    = 0.0076;
  fractionsD0bar[kFlavoredCFBar_0] = 0.305;
  fractionsD0bar[kFlavoredCF_1]    = 0.00031;
  fractionsD0bar[kFlavoredCFBar_1] = 0.144;
  fractionsD0bar[kFlavoredCF_3]    = 0.00027;
  fractionsD0bar[kFlavoredCFBar_3] = 0.0822;
  fractionsD0bar[kFlavoredCS]      = 0.031;
  fractionsD0bar[kFlavoredCSBar]   = 0.031;
  fractionsD0bar[kSLPlus]          = 0.;
  fractionsD0bar[kSLMinus]         = 0.125;
  fractionsD0bar[kCPPlus]          = 0.113;
  fractionsD0bar[kCPMinus]         = 0.102;
  fractionsD0bar[kDalitz]          = 0.056;
  fractionsD0bar[k2Pip2Pim]        = 0.00720;
  fractionsD0bar[kPipPim2Pi0]      = 0.0102;
  fractionsD0bar[kK0PiPiPi0]       = 0.104;
  fractionsD0bar[kKKPiPi]          = 0.00273;
  fractionsD0bar[kK0KK]            = 0.00902;
  fractionsD0bar[kPipPimPi0]       = 0.0149;
  HepVector weight_norm            = fractionsD0bar.T() * scales;
  log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg;
 
  // Original MC
  // Branching fraction predictions after weight
  double brCF_0 =
      ( fractionsD0bar[kFlavoredCFBar_0] * scales[kFlavoredCFBar_0] ) / weight_norm[0];
  double brCF_1 =
      ( fractionsD0bar[kFlavoredCFBar_1] * scales[kFlavoredCFBar_1] ) / weight_norm[0];
  double brCF_3 =
      ( fractionsD0bar[kFlavoredCFBar_3] * scales[kFlavoredCFBar_3] ) / weight_norm[0];
  double brCS = ( fractionsD0bar[kFlavoredCS] * scales[kFlavoredCS] +
                  fractionsD0bar[kFlavoredCSBar] * scales[kFlavoredCSBar] ) /
                weight_norm[0];
  double brDCS_0  = ( fractionsD0bar[kFlavoredCF_0] * scales[kFlavoredCF_0] ) / weight_norm[0];
  double brDCS_1  = ( fractionsD0bar[kFlavoredCF_1] * scales[kFlavoredCF_1] ) / weight_norm[0];
  double brDCS_3  = ( fractionsD0bar[kFlavoredCF_3] * scales[kFlavoredCF_3] ) / weight_norm[0];
  double brCPPlus = ( fractionsD0bar[kCPPlus] * scales[kCPPlus] ) / weight_norm[0];
  double brCPMinus = ( fractionsD0bar[kCPMinus] * scales[kCPMinus] ) / weight_norm[0];
 
  // y=0.0 dalitz example used for prediction
  double dalitz_y_num = -0.000177536;
  double dalitz_y_den = -0.0150846;
 
  double numFactCF_0 = -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x );
  double numFactCF_1 = -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x );
  double numFactCF_3 = -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x );
  double numFactCS   = -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x );
  double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0;
  double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1;
  double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3;
  double denFactCS   = 0.5 * rCS2 * zCS2;
 
  double num_pre = brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 +
                   brCF_3 * numFactCF_3 + 0.5 * brCS * numFactCS + dalitz_y_num;
  double den_pre = 1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 -
                   brCF_3 * denFactCF_3 - 0.5 * brCS * denFactCS + dalitz_y_den;
  double y_pre = num_pre / den_pre;
  double num   = fractionsD0bar[kCPPlus] - fractionsD0bar[kCPMinus] +
               fractionsD0bar[kFlavoredCFBar_0] * numFactCF_0 +
               fractionsD0bar[kFlavoredCFBar_1] * numFactCF_1 +
               fractionsD0bar[kFlavoredCFBar_3] * numFactCF_3 +
               fractionsD0bar[kFlavoredCS] * numFactCS + dalitz_y_num;
  double den = 1. - fractionsD0bar[kCPPlus] - fractionsD0bar[kCPMinus] -
               fractionsD0bar[kFlavoredCFBar_0] * denFactCF_0 -
               fractionsD0bar[kFlavoredCFBar_1] * denFactCF_1 -
               fractionsD0bar[kFlavoredCFBar_3] * denFactCF_3 -
               fractionsD0bar[kFlavoredCS] * denFactCS + dalitz_y_den;
  double y_prematrix = num / den;
  double y_test1     = 0.25 * ( ( ( scales[kCPMinus] * m_weights[kCPPlus][kSLPlus] ) /
                              ( scales[kCPPlus] * m_weights[kCPMinus][kSLPlus] ) ) -
                            ( ( scales[kCPPlus] * m_weights[kCPMinus][kSLPlus] ) /
                              ( scales[kCPMinus] * m_weights[kCPPlus][kSLPlus] ) ) );
 
  log << MSG::INFO << "An Input  Y of " << m_y << endmsg;
  log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg;
  log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg;
  log << MSG::INFO << "The CP-SL double tag method predicts y of " << y_test1 << endmsg;
 
  // Renormalize
  m_largestWeight = 0.;
  for ( int i = 0; i < kNDecayTypes; ++i )
  {
    for ( int j = 0; j < kNDecayTypes; ++j )
    {
      if ( m_weights[i][j] < 0 )
        m_weights[i][j] = 0; // gets rid of negative weights in early calculations
      if ( m_weights[i][j] > m_largestWeight ) { m_largestWeight = m_weights[i][j]; }
    }
  }
  m_weights /= m_largestWeight;
 
  log << MSG::INFO << "Renormalized weight matrix " << m_weights << endmsg;
 
  // Set up weight parameters
  Dalitz   t( 8 );
  double   D0Sum[10], D0barSum[10]; // last index is for sum over bins
  TComplex D0D0barSum[10];
  int      points[10];
 
  for ( int i = 0; i < 10; ++i )
  {
    D0Sum[i]      = 0.;
    D0barSum[i]   = 0.;
    D0D0barSum[i] = TComplex( 0., 0. );
    points[i]     = 0;
  }
 
  double m2min    = 0.;
  double m2max    = 3.;
  int    nsteps   = 1000;
  double stepsize = ( m2max - m2min ) / (double)nsteps;
 
  for ( int i = 0; i < nsteps; ++i )
  {
    double m2i = m2min + ( (double)i + 0.5 ) * stepsize;
 
    for ( int j = 0; j < nsteps; ++j )
    {
      double m2j = m2min + ( (double)j + 0.5 ) * stepsize;
 
      if ( t.Point_on_DP( m2i, m2j ) )
      {
        double m2k = 1.86450 * 1.86450 + 0.497648 * 0.497648 + 0.139570 * 0.139570 +
                     0.139570 * 0.139570 - m2i - m2j;
 
        TComplex D0, D0bar;
        if ( m_dalitzModel == 0 )
        { // Cleo model
          D0    = t.CLEO_Amplitude( m2i, m2j, m2k );
          D0bar = t.CLEO_Amplitude( m2j, m2i, m2k );
        }
        if ( m_dalitzModel == 1 )
        { // Babar model
          D0    = t.Babar_Amplitude( m2i, m2j, m2k );
          D0bar = t.Babar_Amplitude( m2j, m2i, m2k );
        }
        if ( m_dalitzModel == 2 )
        { // Belle model
          D0    = t.Amplitude( m2i, m2j, m2k );
          D0bar = t.Amplitude( m2j, m2i, m2k );
        }
 
        if ( j <= i )
        { // lower half only
          int bin = t.getBin( m2i, m2j, m2k );
          if ( bin == -1 ) bin = 8;
 
          ++points[bin];
          D0Sum[bin] += norm( D0 );
          D0barSum[bin] += norm( D0bar );
          D0D0barSum[bin] += D0 * conj( D0bar );
 
          ++points[9];
          D0Sum[9] += norm( D0 );
          D0barSum[9] += norm( D0bar );
          D0D0barSum[9] += D0 * conj( D0bar );
        }
      }
    }
  }
 
  for ( int i = 0; i < 10; ++i )
  {
    double r2 = D0barSum[i] / D0Sum[i];
    double c  = real( D0D0barSum[i] ) / sqrt( D0barSum[i] ) / sqrt( D0Sum[i] );
    double s  = imag( D0D0barSum[i] ) / sqrt( D0barSum[i] ) / sqrt( D0Sum[i] );
 
    cout << "BIN " << i << ": " << points[i] << " pts " << r2 << " " << c << " " << s << " "
         << c * c + s * s << endmsg;
  }
 
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
//********************************************************************
StatusCode QCMCFilter::execute() {
 
  // interface to event data service
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode   sc         = svcLocator->service( "EventDataSvc", m_evtSvc );
  if ( sc.isFailure() ) std::cout << "Could not access EventDataSvc!" << std::endl;
 
  // setFilterPassed(false);
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
 
  if ( !eventHeader )
  {
    cout << " eventHeader " << endl;
    return StatusCode::FAILURE;
  }
  long runNo = eventHeader->runNumber();
  long event = eventHeader->eventNumber();
  log << MSG::DEBUG << "run, evtnum = " << runNo << " , " << event << endmsg;
 
  unsigned int QC_status = 0;
  // eventHeader->setEventTag(0);
 
  // Get McParticle collection
  SmartDataPtr<Event::McParticleCol> mcParticles( eventSvc(), EventModel::MC::McParticleCol );
 
  // Check if event is Psi(3770)
  int psipp_daughter = 0;
  int d0_count       = 0;
  int chd_count      = 0;
  for ( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end();
        it++ )
  {
    int pdg_code = ( *it )->particleProperty();
    if ( ( *it )->primaryParticle() ) continue;
    Event::McParticle it2        = ( *it )->mother();
    int               mother_pdg = 0;
    // finds if the particle is daughter of Psi(3770)
    mother_pdg = it2.particleProperty();
    if ( mother_pdg == kPsi3770ID )
    {
      psipp_daughter++; // Found daughter
      if ( abs( pdg_code ) == kD0ID ) d0_count++;
      if ( abs( pdg_code ) == kDpID ) chd_count++;
      if ( pdg_code == kD0BarID ) m_nD0bar++;
    }
  }
 
  if ( psipp_daughter == 0 ) // didn't find Psi(3770)
  {
    ++m_nUnknownEvents;
 
    if ( m_invertSelection )
    {
      // eventHeader->setEventTag(0);
      QC_status = 0;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      log << MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg;
      return StatusCode::SUCCESS; // discard event
    }
    else
    {
      // -------- Write to root file
      // setFilterPassed(true);
      // eventHeader->setEventTag(1);
      QC_status = 1;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
 
      return StatusCode::SUCCESS; // kept event
    }
  }
 
  log << MSG::DEBUG << "Found psi(3770) -->" << endmsg;
 
  // Check that top particle has only two daughters + possible FSR photon
  if ( psipp_daughter > 3 )
  {
    ++m_nUnknownEvents;
 
    if ( m_invertSelection )
    {
      // eventHeader->setEventTag(0);
      QC_status = 0;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      log << MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg;
      return StatusCode::SUCCESS; // discard event
    }
    else
    {
      // -------- Write to root file
      // eventHeader->setEventTag(1);
      QC_status = 1;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      // setFilterPassed(true);
 
      return StatusCode::SUCCESS; // kept event
    }
  }
 
  // Check if D+D- state
  if ( chd_count == 2 )
  { // D+D- event
    ++m_nDpDmEvents;
 
    // D+D- must be rewieghted otherwise too many D+D- events relative to D0D0bar.
    double random = RandFlat::shoot();
 
    if ( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
         ( random > m_dplusDminusWeight && m_invertSelection ) )
    {
      // -------- Write to root file
      // setFilterPassed(true);
      // eventHeader->setEventTag(1);
      QC_status = 1;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
 
      return StatusCode::SUCCESS; // event kept
    }
    else
    {
      // eventHeader->setEventTag(0);
      QC_status = 0;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      log << MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg;
      ++m_nDpDmDiscarded;
      return StatusCode::SUCCESS; // event discarded
    }
  }
 
  // Check if D0D0Bar
  if ( d0_count != 2 )
  {
    if ( !m_invertSelection )
    {
      // eventHeader->setEventTag(0);
      QC_status = 0;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      log << MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg;
      return StatusCode::SUCCESS; // discard event
    }
    else
    {
      // -------- Write to root file
      // eventHeader->setEventTag(1);
      QC_status = 1;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      // setFilterPassed(true);
 
      return StatusCode::SUCCESS; // kept event
    }
  }
 
  log << MSG::INFO << "D0-D0bar event." << endmsg;
  m_nD0D0barEvents++;
 
  m_rho_flag = 0;
 
  // Get D0-D0bar modes and info
  vector<double> d0_decay    = findD0Decay( 1 );
  vector<double> d0bar_decay = findD0Decay( -1 );
 
  int d0mode    = d0_decay[0];
  int d0barmode = d0bar_decay[0];
  if ( m_debug ) cout << "D0 " << d0mode << endl;
  if ( m_debug ) cout << "D0bar " << d0barmode << endl;
  m_modeCounter[d0mode][d0barmode]++;
 
  // if any event unidentified remove
  if ( d0_decay[0] == 20 || d0bar_decay[0] == 20 )
  {
    if ( !m_invertSelection )
    {
      // eventHeader->setEventTag(0);
      QC_status = 0;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      log << MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg;
      return StatusCode::SUCCESS; // discard event
    }
    else
    {
      // -------- Write to root file
      // eventHeader->setEventTag(1);
      QC_status = 1;
      eventHeader->setEventTag( QC_status );
 
      log << MSG::INFO << "QC_status " << QC_status << endmsg;
      // setFilterPassed(true);
 
      return StatusCode::SUCCESS; // kept event
    }
  }
 
  // Loop over D's
  double r2D0       = d0_decay[1];
  double deltaD0    = d0_decay[2];
  double RD0        = d0_decay[3];
  double FCP        = d0_decay[4];
  double r2D0bar    = d0bar_decay[1];
  double deltaD0bar = d0bar_decay[2];
  double RD0bar     = d0bar_decay[3];
  double FCPbar     = d0bar_decay[4];
 
  // Weight based on DDbar combination
  double weight;
  if ( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz || d0_decay[0] == k2Pip2Pim ||
       d0bar_decay[0] == k2Pip2Pim || d0_decay[0] == kPipPim2Pi0 ||
       d0bar_decay[0] == kPipPim2Pi0 || d0_decay[0] == kK0PiPiPi0 ||
       d0bar_decay[0] == kK0PiPiPi0 || d0bar_decay[0] == kPipPimPi0 ||
       d0_decay[0] == kKKPiPi || d0bar_decay[0] == kKKPiPi || d0_decay[0] == kK0KK ||
       d0bar_decay[0] == kK0KK || d0bar_decay[0] == kFlavoredCF_3 ||
       d0bar_decay[0] == kFlavoredCFBar_3 || d0_decay[0] == kPipPimPi0 )
  {
    double r2prod = r2D0 * r2D0bar;
    double x      = m_x;
    double y      = m_y;
 
    // double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
    // double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) /
    // bD0 ; double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar )
    // * x ) ; double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin(
    // deltaD0bar ) * x ) ) / bD0bar ;
 
    // weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
    // weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
    double bD0 = 1. - ( 2 * FCP - 1. ) * RD0 * sqrt( r2D0 ) *
                          ( cos( deltaD0 ) * y + sin( deltaD0 ) * x );
    double rwsD0 = ( r2D0 - ( 2 * FCP - 1. ) * RD0 * sqrt( r2D0 ) *
                                ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) /
                   bD0;
    double bD0bar = 1. - ( 2 * FCPbar - 1. ) * RD0bar * sqrt( r2D0bar ) *
                             ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x );
    double rwsD0bar = ( r2D0bar - ( 2 * FCPbar - 1. ) * RD0bar * sqrt( r2D0bar ) *
                                      ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) /
                      bD0bar;
 
    log << MSG::INFO << "bD0 == " << bD0 << endmsg;
    log << MSG::INFO << "bD0bar == " << bD0bar << endmsg;
    log << MSG::INFO << "rwsD0 == " << rwsD0 << endmsg;
    log << MSG::INFO << "rwsD0bar == " << rwsD0bar << endmsg;
    weight = 1. + r2prod -
             2. * RD0 * ( 2 * FCP - 1. ) * ( 2 * FCPbar - 1. ) * sqrt( r2prod ) *
                 cos( deltaD0 + deltaD0bar );
    log << MSG::INFO << "weight == " << weight << endmsg;
    weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar;
    log << MSG::INFO << "weight == " << weight << endmsg;
    weight /= m_largestWeight;
    log << MSG::INFO << "weight == " << weight << endmsg;
 
    // get dalitz graph information before cut
    int    k   = -10;
    double m2i = 0;
    double m2j = 0;
    double m2k = 0;
    if ( d0_decay[0] == kDalitz )
    {
      k   = d0bar_decay[0];
      m2i = d0_decay[4];
      m2j = d0_decay[5];
      m2k = d0_decay[7];
    }
    if ( d0bar_decay[0] == kDalitz )
    {
      k   = d0_decay[0];
      m2i = d0bar_decay[4];
      m2j = d0bar_decay[5];
      m2k = d0bar_decay[7];
    }
  }
  else { weight = m_weights[d0_decay[0]][d0bar_decay[0]]; }
 
  double random = RandFlat::shoot();
  log << MSG::DEBUG << "type:  " << d0_decay[0] << " " << d0bar_decay[0] << ", weight "
      << weight << " <? random " << random << endmsg;
 
  if ( ( random < weight && !m_invertSelection ) || ( random > weight && m_invertSelection ) )
  {
    // -------- Write to root file
    QC_status = 1;
    eventHeader->setEventTag( QC_status );
 
    log << MSG::INFO << "QC_status " << QC_status << endmsg;
    // eventHeader->setEventTag(1);
    // setFilterPassed(true);
 
    // get dalitz graph information after filter
    /*int k = -10;
            double m2i= 0;
            double m2j= 0;
            double m2k= 0;
            double cosd = cos( deltaD0 ) ;
            double sind = sin( deltaD0 ) ;
            if( d0_decay[0] == kDalitz) {
            k = d0bar_decay[0];
            m2i= d0_decay[4];
            m2j= d0_decay[5];
            m2k= d0_decay[7];
            if ( m2j - m2i > 0. ) {  //avoids doublecounting the branching fraction
            dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
            dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
            dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
            }
            if( d0bar_decay[0] == kDalitz) {
            k = d0_decay[0];
            m2i= d0bar_decay[4];
            m2j= d0bar_decay[5];
            m2k= d0bar_decay[7];
            cosd = cos( deltaD0bar ) ;
            sind = sin( deltaD0bar ) ;
            if( m2j - m2i < 0. ) { //avoids doublecounting the branching fraction
            dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
            dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
            dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
            }
 
            m_keptModeCounter[d0mode][d0barmode]++;*/
    return StatusCode::SUCCESS; // event kept
  }
  else
  {
    QC_status = 0;
    eventHeader->setEventTag( QC_status );
 
    log << MSG::INFO << "QC_status " << QC_status << endmsg;
    // eventHeader->setEventTag(0);
    log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg;
 
    ++m_nD0D0barDiscarded;
    return StatusCode::SUCCESS; // event discarded
  }
}
 
//*********************************************************************
StatusCode QCMCFilter::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  log << MSG::INFO << "Number of unknown events: " << m_nUnknownEvents << endmsg;
  log << MSG::INFO << "Number of D+D- events seen: " << m_nDpDmEvents << endmsg;
  log << MSG::INFO << "Number of D+D- events discarded: " << m_nDpDmDiscarded << endmsg;
  if ( m_nDpDmEvents > 0 )
  {
    log << MSG::INFO
        << "Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents )
        << endmsg;
  }
 
  log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar
      << endmsg;
  log << MSG::INFO << "Number of D0D0bar events discarded: " << m_nD0D0barDiscarded << endmsg;
  /*if( m_nD0D0barEvents > 0 &&  m_nCPPlus > 0 && m_nCPMinus > 0)
          {
          log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double(
m_nD0D0barEvents ) << endl << "Fraction unknown decays: " <<    0.5 * double( m_nUnknownDecays
) / double( m_nD0D0barEvents ) <<endmsg ;
 
          double nD0 = 2. * m_nD0D0barEvents ;
          double brSL = m_nSL / nD0 ;
          double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
          double brCF_0 = m_nFlavoredCFD0_0 / nD0 ;
          double dBrCF_0 = sqrt( brCF_0 * ( 1. - brCF_0 ) / nD0 ) ;
          double brCF_1 = m_nFlavoredCFD0_1 / nD0 ;
          double dBrCF_1 = sqrt( brCF_1 * ( 1. - brCF_1 ) / nD0 ) ;
          double brCF_3 = m_nFlavoredCFD0_3 / nD0 ;
          double dBrCF_3 = sqrt( brCF_3 * ( 1. - brCF_3 ) / nD0 ) ;
          double brCS = m_nFlavoredCSD0 / nD0 ;
          double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
          double brDCS_0 = m_nFlavoredDCSD0_0 / nD0 ;
          double dBrDCS_0 = sqrt( brDCS_0 * ( 1. - brDCS_0 ) / nD0 ) ;
          double brDCS_1 = m_nFlavoredDCSD0_1 / nD0 ;
          double dBrDCS_1 = sqrt( brDCS_1 * ( 1. - brDCS_1 ) / nD0 ) ;
          double brDCS_3 = m_nFlavoredDCSD0_3 / nD0 ;
          double dBrDCS_3 = sqrt( brDCS_3 * ( 1. - brDCS_3 ) / nD0 ) ;
          double brCPPlus = m_nCPPlus / nD0 ;
          double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
          double brCPMinus = m_nCPMinus / nD0 ;
          double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
          double brDalitz = m_nDalitz / nD0 ;
          double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
          double fdBrDalitz = dBrDalitz / brDalitz ;
          double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
          double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
          double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
          double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
          double dalitzDenomNorm = m_dalitzDenom / nD0 ;
          double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
 
          double br2Pip2Pim = m_n2Pip2Pim / nD0 ;
          double dBr2Pip2Pim = sqrt( br2Pip2Pim * ( 1. - br2Pip2Pim ) / nD0 ) ;
          double fdBr2Pip2Pim= dBr2Pip2Pim / br2Pip2Pim ;
          double Pip2Pim2Numer1Norm = m_2Pip2PimNumer1 / nD0 ;
          double d2Pip2PimNumer1Norm = Pip2Pim2Numer1Norm * fdBr2Pip2Pim ;
          double Pip2Pim2Numer2Norm = m_2Pip2PimNumer2 / nD0 ;
          double d2Pip2PimNumer2Norm = Pip2Pim2Numer2Norm * fdBr2Pip2Pim ;
          double Pip2Pim2DenomNorm = m_2Pip2PimDenom / nD0 ;
          double d2Pip2PimDenomNorm = Pip2Pim2DenomNorm * fdBr2Pip2Pim ;
 
          double brPipPim2Pi0 = m_nPipPim2Pi0 / nD0 ;
          double dBrPipPim2Pi0 = sqrt( brPipPim2Pi0 * ( 1. - brPipPim2Pi0 ) / nD0 ) ;
          double fdBrPipPim2Pi0= dBrPipPim2Pi0 / brPipPim2Pi0 ;
          double Pip2Pim2Numer1Norm = m_PipPim2Pi0Numer1 / nD0 ;
          double dPipPim2Pi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrPipPim2Pi0 ;
          double Pip2Pim2Numer2Norm = m_PipPim2Pi0Numer2 / nD0 ;
          double dPipPim2Pi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrPipPim2Pi0 ;
          double Pip2Pim2DenomNorm = m_PipPim2Pi0Denom / nD0 ;
          double dPipPim2Pi0DenomNorm = Pip2Pim2DenomNorm * fdBrPipPim2Pi0 ;
 
          double brKSPiPiPi0 = m_nKSPiPiPi0 / nD0 ;
          double dBrKSPiPiPi0 = sqrt( brKSPiPiPi0 * ( 1. - brKSPiPiPi0 ) / nD0 ) ;
          double fdBrKSPiPiPi0= dBrKSPiPiPi0 / brKSPiPiPi0 ;
          double Pip2Pim2Numer1Norm = m_KSPiPiPi0Numer1 / nD0 ;
          double dKSPiPiPi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrKSPiPiPi0 ;
          double Pip2Pim2Numer2Norm = m_KSPiPiPi0Numer2 / nD0 ;
          double dKSPiPiPi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrKSPiPiPi0 ;
          double Pip2Pim2DenomNorm = m_KSPiPiPi0Denom / nD0 ;
          double dKSPiPiPi0DenomNorm = Pip2Pim2DenomNorm * fdBrKSPiPiPi0 ;
 
  //after the filter was applied
  double fil_SL = 0, fil_FlavoredCFD0_0 = 0, fil_FlavoredCFD0_1 = 0, fil_FlavoredCFD0_3 = 0,
fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0_0 =0, fil_FlavoredDCSD0_1 =0, fil_FlavoredDCSD0_3 =0,
fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0, fil_2Pip2Pim = 0; for( int i = 0 ; i < 14 ;
++i )
  {
  fil_SL += m_keptModeCounter[i][9] + m_keptModeCounter[9][i] + m_keptModeCounter[i][8] +
m_keptModeCounter[8][i] ; fil_FlavoredCFD0_0 += m_keptModeCounter[i][1] +
m_keptModeCounter[0][i]; fil_FlavoredCFD0_1 += m_keptModeCounter[i][3] +
m_keptModeCounter[2][i]; fil_FlavoredCFD0_3 += m_keptModeCounter[i][5] +
m_keptModeCounter[4][i]; fil_FlavoredCSD0 += m_keptModeCounter[i][7] + m_keptModeCounter[7][i]
+ m_keptModeCounter[i][6] + m_keptModeCounter[6][i]; fil_FlavoredDCSD0_0 +=
m_keptModeCounter[i][0] + m_keptModeCounter[1][i]; fil_FlavoredDCSD0_1 +=
m_keptModeCounter[i][2] + m_keptModeCounter[3][i]; fil_FlavoredDCSD0_3 +=
m_keptModeCounter[i][4] + m_keptModeCounter[5][i]; fil_CPPlus += m_keptModeCounter[i][10] +
m_keptModeCounter[10][i]; fil_CPMinus += m_keptModeCounter[i][11] + m_keptModeCounter[11][i];
  fil_Dalitz += m_keptModeCounter[i][12] + m_keptModeCounter[12][i];
  fil_2Pip2Pim += m_keptModeCounter[i][13] + m_keptModeCounter[13][i];
  fil_PipPim2Pi0 += m_keptModeCounter[i][14] + m_keptModeCounter[14][i];
  fil_KSPiPiPi0 += m_keptModeCounter[i][15] + m_keptModeCounter[15][i];
}
double nD0_fil = 2. * ( m_nD0D0barEvents -  m_nD0D0barDiscarded ) ;
double brSL_fil = fil_SL / nD0_fil ;
double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
double brCF_0_fil = fil_FlavoredCFD0_0 / nD0_fil ;
double dBrCF_0_fil = sqrt( brCF_0_fil * ( 1. - brCF_0_fil ) / nD0_fil ) ;
double brCF_1_fil = fil_FlavoredCFD0_1 / nD0_fil ;
double dBrCF_1_fil = sqrt( brCF_1_fil * ( 1. - brCF_1_fil ) / nD0_fil ) ;
double brCF_3_fil = fil_FlavoredCFD0_3 / nD0_fil ;
double dBrCF_3_fil = sqrt( brCF_3_fil * ( 1. - brCF_3_fil ) / nD0_fil ) ;
double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
double brDCS_0_fil = fil_FlavoredDCSD0_0 / nD0_fil ;
double brDCS_1_fil = fil_FlavoredDCSD0_1 / nD0_fil ;
double brDCS_3_fil = fil_FlavoredDCSD0_3 / nD0_fil ;
double dBrDCS_0_fil = sqrt( brDCS_0_fil * ( 1. - brDCS_0_fil ) / nD0_fil ) ;
double dBrDCS_1_fil = sqrt( brDCS_1_fil * ( 1. - brDCS_1_fil ) / nD0_fil ) ;
double dBrDCS_3_fil = sqrt( brDCS_3_fil * ( 1. - brDCS_3_fil ) / nD0_fil ) ;
double brCPPlus_fil = fil_CPPlus / nD0_fil ;
double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
double brCPMinus_fil = fil_CPMinus / nD0_fil ;
double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
double brDalitz_fil = fil_Dalitz / nD0_fil ;
double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
 
double br2Pip2Pim_fil = fil_2Pip2Pim / nD0_fil ;
double dBr2Pip2Pim_fil = sqrt( br2Pip2Pim_fil * ( 1. - br2Pip2Pim_fil ) / nD0_fil ) ;
double fdBr2Pip2Pim_fil = dBr2Pip2Pim_fil / br2Pip2Pim_fil ;
double Pip2Pim2Numer1Norm_fil = Pip2Pim2Numer1_fil / nD0_fil ;
double d2Pip2PimNumer1Norm_fil = Pip2Pim2Numer1Norm_fil * fdBr2Pip2Pim_fil ;
double Pip2Pim2Numer2Norm_fil = Pip2Pim2Numer2_fil / nD0_fil ;
double d2Pip2PimNumer2Norm_fil = Pip2Pim2Numer2Norm_fil * fdBr2Pip2Pim_fil ;
double Pip2Pim2DenomNorm_fil = Pip2Pim2Denom_fil / nD0_fil ;
double d2Pip2PimDenomNorm_fil = Pip2Pim2DenomNorm_fil * fdBr2Pip2Pim_fil ;
 
log << MSG::INFO <<
"Original number of SL decays: " << m_nSL <<"  ==> Br(SL) = " << brSL << " +- " << dBrSL <<
endl << "Filtered number of SL decays: " << fil_SL <<"  ==> Br(SL) = " << brSL_fil << " +- " <<
dBrSL_fil << endl << "Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: "
<<m_nFlavoredCFD0_0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0_0 << endl << "  ==>
Br(CF) = " << brCF_0 << " +- " << dBrCF_0 << endl << "Original number of D0->CF1/DCS1 or
D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_1  << "/" << m_nFlavoredDCSD0_1 << endl << "  ==>
Br(CF) = " << brCF_1 << " +- " << dBrCF_1 << endl << "Original number of D0->CF1/DCS1 or
D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_3  << "/" << m_nFlavoredDCSD0_3 << endl << "  ==>
Br(CF) = " << brCF_3 << " +- " << dBrCF_3 << endl << "  ==> Br(CS) = " << brCS << " +- " <<
dBrCS << endl << "  ==> Br(DCS) = " << brDCS_0 << " +- " << dBrDCS_0 << endl << "Filtered
number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0_0 << "/" <<
fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0_0 << endl << "  ==> Br(DCS) = " << brDCS_1 << " +-
" << dBrDCS_1 << endl << "Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: "
<<fil_FlavoredCFD0_1 << "/" << fil_FlavoredDCSD0_1 << endl << "  ==> Br(DCS) = " << brDCS_3 <<
" +- " << dBrDCS_3 << endl << "Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: "
<<fil_FlavoredCFD0_3 << "/" << fil_FlavoredDCSD0_3 << endl << "  ==> Br(CF) = " << brCF_0_fil
<< " +- " << dBrCF_0_fil << endl << "  ==> Br(CF) = " << brCF_1_fil << " +- " << dBrCF_1_fil <<
endl << "  ==> Br(CF) = " << brCF_3_fil << " +- " << dBrCF_3_fil << endl << "  ==> Br(CS) = "
<< brCS_fil << " +- " << dBrCS_fil << endl << "  ==> Br(DCS) = " << brDCS_0_fil << " +- " <<
dBrDCS_0_fil << endl << "  ==> Br(DCS) = " << brDCS_1_fil << " +- " << dBrDCS_1_fil << endl <<
"  ==> Br(DCS) = " << brDCS_3_fil << " +- " << dBrDCS_3_fil << endl <<
 
"Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl <<
"  ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl <<
"  ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl <<
"Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl <<
"  ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl <<
"  ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl <<
 
"Original number of Dalitz: " << m_nDalitz << endl <<
"  ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
"  y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
"  numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
"  denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
"Filtered number of Dalitz: " << fil_Dalitz << endl <<
"  ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
"  y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
"  numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
"  denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg <<
 
"Original number of 2Pip2Pim: " << m_n2Pip2Pim << endl <<
"  ==> Br(2Pip2Pim) = " << br2Pip2Pim << " +- " << dBr2Pip2Pim <<endl <<
"  y contrib. numer1 " << Pip2Pim2Numer1Norm << " +- " << d2Pip2PimNumer1Norm << endl <<
"  numer2 " << Pip2Pim2Numer2Norm << " +- " << d2Pip2PimNumer2Norm << endl <<
"  denom " << Pip2Pim2DenomNorm << " +- " << d2Pip2PimDenomNorm << endmsg <<
"Filtered number of 2Pip2Pim: " << fil_2Pip2Pim << endl <<
"  ==> Br(2Pip2Pim) = " << br2Pip2Pim_fil << " +- " << dBr2Pip2Pim_fil <<endl <<
"  y contrib. numer1 " << Pip2Pim2Numer1Norm_fil << " +- " << d2Pip2PimNumer1Norm_fil << endl
<< "  numer2 " << Pip2Pim2Numer2Norm_fil << " +- " << d2Pip2PimNumer2Norm_fil << endl << "
denom " << Pip2Pim2DenomNorm_fil << " +- " << d2Pip2PimDenomNorm_fil << endmsg ;
 
 
HepMatrix    differencetoWeight(17,17,0);
for( int i = 0 ; i < 16 ; ++i )  {
  for( int j = 0 ; j < 16 ; ++j )  {
          if (m_modeCounter[i][j] != 0) {
                  differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] -
m_weights[i][j]; } } }
 
                  log << MSG::INFO << "Weight matrix" << m_weights <<  endmsg ;
                  log << MSG::INFO << "D0 Modes before filter" << m_modeCounter <<  endmsg ;
                  log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter <<  endmsg ;
                  log << MSG::INFO << "Compare percentage difference to weights " <<
differencetoWeight <<  endmsg ;
 
                  //Calculating and comparing Y before/after the filter
                  // To calculate y, we want half the CS BF
                  brCS /= 2. ;
                  dBrCS /= 2. ;
                  brCS_fil /= 2 ;
                  dBrCS_fil /= 2 ;
 
                  double y, dy, y_fil, dy_fil ;
                  double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
 
                  //Original MC
                  double rCF_0 = m_rCF_0 ;
                  double zCF_0 = m_zCF_0 ;
                  double rCF2_0 = rCF_0 * rCF_0 ;
                  double zCF2_0 = zCF_0 * zCF_0 ;
                  double rCF_1 = m_rCF_1 ;
                  double zCF_1 = m_zCF_1 ;
                  double rCF2_1 = rCF_1 * rCF_1 ;
                  double zCF2_1 = zCF_1 * zCF_1 ;
                  double rCF_3 = m_rCF_3 ;
                  double zCF_3 = m_zCF_3 ;
                  double rCF2_3 = rCF_3 * rCF_3 ;
                  double zCF2_3 = zCF_3 * zCF_3 ;
                  double rCS = m_rCS ;
                  double zCS = m_zCS ;
                  double rCS2 = rCS * rCS ;
                  double zCS2 = zCS * zCS ;
 
                  double wCF_0 = ( m_wCFSign_0 ? 1. : -1. ) * sqrt( 4. - zCF2_0 ) ;
                  double wCF_1 = ( m_wCFSign_1 ? 1. : -1. ) * sqrt( 4. - zCF2_1 ) ;
                  double wCF_3 = ( m_wCFSign_3 ? 1. : -1. ) * sqrt( 4. - zCF2_3 ) ;
                  double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
 
                  double numFactCF_0 =
                  -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x ) ;
                  double numFactCF_1 =
                  -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x ) ;
                  double numFactCF_3 =
                  -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x ) ;
                  double numFactCS =
                  -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
                  double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0 ;
                  double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1 ;
                  double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3 ;
                  double denFactCS = 0.5 * rCS2 * zCS2 ;
                  double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ;
                  double Pip2Pim2NumerNorm = Pip2Pim2Numer1Norm + Pip2Pim2Numer2Norm ;
 
                  double num =
                  brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 + brCF_3 *
numFactCF_3 + brCS * numFactCS + dalitzNumerNorm + Pip2Pim2NumerNorm; double den =
                  1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 -
brCF_3 * denFactCF_3 - brCS * denFactCS + dalitzDenomNorm + Pip2Pim2DenomNorm; y = num / den ;
                  dy = sqrt(
                                  ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
                                  ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
                                  ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 *
                                  ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 +
                                  ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 *
                                  ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 +
                                  ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 *
                                  ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 +
                                  ( numFactCS * den + denFactCS * num ) * dBrCS *
                                  ( numFactCS * den + denFactCS * num ) * dBrCS +
                                  ( dalitzNumerNorm * den + dalitzDenomNorm * num ) *
fdBrDalitz * ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz + (
Pip2Pim2NumerNorm * den + Pip2Pim2DenomNorm * num ) * fdBrDalitz * ( Pip2Pim2NumerNorm * den +
Pip2Pim2DenomNorm * num ) * fdBrDalitz ) / ( den * den ) ;
 
                  double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] +
m_modeCounter[4][6] +  m_modeCounter[5][6] ; double n_cpminussl = m_modeCounter[7][4] +
m_modeCounter[7][5] + m_modeCounter[4][7] +  m_modeCounter[5][7] ; double a_cpsl = ( n_cpplussl
* m_nCPMinus ) / ( n_cpminussl * m_nCPPlus  ) ; double b_cpsl = ( n_cpminussl  * m_nCPPlus ) /
( n_cpplussl * m_nCPMinus  ) ; y_cpsl = 0.25 * ( a_cpsl  - b_cpsl ) ; dy_cpsl = 0.25 * ( a_cpsl
+ b_cpsl )  * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) )
;
 
 
                  //Filtered MC
                  double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
                  double Pip2Pim2NumerNorm_fil = Pip2Pim2Numer1Norm_fil +
Pip2Pim2Numer2Norm_fil ; double num_fil = brCPPlus_fil - brCPMinus_fil + brCF_0_fil *
numFactCF_0 + brCF_1_fil * numFactCF_1 + brCF_3_fil * numFactCF_3 + brCS_fil * numFactCS +
                  dalitzNumerNorm_fil + Pip2Pim2NumerNorm_fil;
                  double den_fil =
                  1. - brCPPlus_fil - brCPMinus_fil - brCF_0_fil * denFactCF_0 - brCF_1_fil *
denFactCF_1 - brCF_3_fil * denFactCF_3 - brCS_fil * denFactCS + dalitzDenomNorm_fil +
Pip2Pim2DenomNorm_fil; y_fil = num_fil / den_fil ; dy_fil = sqrt( ( num_fil + den_fil ) * (
num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil + ( num_fil - den_fil ) * ( num_fil -
den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil + ( numFactCF_0 * den_fil + denFactCF_0 * num_fil )
* dBrCF_0_fil * ( numFactCF_0 * den_fil + denFactCF_0 * num_fil ) * dBrCF_0_fil + ( numFactCF_1
* den_fil + denFactCF_1 * num_fil ) * dBrCF_1_fil * ( numFactCF_1 * den_fil + denFactCF_1 *
num_fil ) * dBrCF_1_fil + ( numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil * (
numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil + ( numFactCS * den_fil +
denFactCS * num_fil ) * dBrCS_fil * ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
                                  ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil *
num_fil ) * fdBrDalitz_fil * ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil )
* fdBrDalitz_fil + ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) *
fdBr2Pip2Pim_fil * ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) *
fdBr2Pip2Pim_fil ) / ( den_fil * den_fil ) ;
 
                  double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] +
m_keptModeCounter[4][6] +  m_keptModeCounter[5][6] ; double fil_cpminussl =
m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] +
m_keptModeCounter[5][7] ; a_cpsl = ( fil_cpplussl  * fil_CPMinus ) / ( fil_cpminussl *
fil_CPPlus  ) ; b_cpsl = ( fil_cpminussl  * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus  ) ;
                  y_fil_cpsl = 0.25 * ( a_cpsl  - b_cpsl ) ;
                  dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + (
1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ;
 
 
 
                  log << MSG::INFO <<"BR Method : Parent MC    ==> y = " << y << " +- " << dy
<< endmsg ; log << MSG::INFO <<"BR Method : Filtered MC  ==> y = " << y_fil << " +- " << dy_fil
<< endmsg ; log << MSG::INFO <<"CP-SL doubletag : Parent MC   ==> y = " << y_cpsl << " +- " <<
dy_cpsl <<endl<<"Previous line should be equivalent with 0 as parent MC not correalated"<<
endmsg ; log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " <<
dy_fil_cpsl << endmsg ;
 
                  }*/
 
  return StatusCode::SUCCESS;
}
 
///////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
vector<double> QCMCFilter::findD0Decay( int charm ) {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "  In findD0Decay " << endmsg;
  vector<double> d0_info( 9 );
  for ( int i = 0; i < 9; i++ ) d0_info[i] = 0;
  double decay_mode = 20; // d0_info[0] = decay mode;
  double r2D0       = -1; // d0_info[1] = r2D0;
  double deltaD0    = -1; // d0_info[2] = deltaD0;
  double RD0        = 1;  // d0_info[3] = RD0;
  double FCP        = 1;  // d0_info[4] = FCP;
 
  double num[26];
  for ( int i = 0; i < 26; i++ ) num[i] = 0;
  // Number of each partile type
  int kPiPlus        = 0; // num[0] = Pi+, a0+
  int kPiMinus       = 1; // num[1] = Pi-, a0-
  int kPi0           = 2; // num[2] = Pi0
  int kEta           = 3; // num[3] = Eta
  int kEtaPrime      = 4; // num[4] = Eta Prime
  int kNeutralScalar = 5; // num[5] = Neutral Scalar : f0, f'0, f0(1500), a0,  (f_2 can be
                          // included here because it is only used in f_2 pi0 decay)
  int kUFVPlus     = 6;   // num[6] = unflavored positive (axial) vectors: rho+, a1+, rho(2S)+
  int kUFVMinus    = 7;   // num[7] = unflavored negative (axial) vectors: rho-, a1-, rho(2S)-
  int kRho0        = 8;   // num[8] = Rho0, Rho(2S)0
  int kOmega       = 9;   // num[9] = Omega
  int kPhi         = 10;  // num[10] = Phi
  int kKPlus       = 11;  // num[11] = K+
  int kKMinus      = 12;  // num[12] = K-
  int kK0S         = 13;  // num[13] = K short
  int kK0L         = 14;  // num[14] = K long
  int kFVPlus      = 15;  // num[15] = flavored positive (axial) vectors: K*+, K_1+, K_2*+
  int kFVMinus     = 16;  // num[16] = flavored negative (axial) vectors: K*-, K_1-, K_2*-
  int kKStar0      = 17;  // num[17] = K*0
  int kKStar0Bar   = 18;  // num[18] = anti-K*0
  int kK10         = 19;  // num[19] = K1_0
  int kK10Bar      = 20;  // num[20] = anti-K1_0
  int kLeptonPlus  = 21;  // num[21] =  LeptonPlus: e+, mu+
  int kLeptonMinus = 22;  // num[22] =  LeptonMinus: e-, mu-
  int kNeutrino    = 23;  // num[23] =  Neutrino: nu, nubar
  int kGamma       = 24;  // num[24] =  gamma
  int kUnknown     = 25;  // num[25] =  Not identified particle
  // int kGammaFSR      GammaFSR is ignored as we are interested in the PreFSR quantum #
 
  // Get McParticle collection
  SmartDataPtr<Event::McParticleCol> mcParticles( eventSvc(), EventModel::MC::McParticleCol );
  int                                d0_pdg = charm == 1 ? kD0ID : kD0BarID;
 
  HepLorentzVector         piPlus, piMinus, k0;
  vector<HepLorentzVector> piPlus_4pi, piMinus_4pi, pi0_p4;
  vector<HepLorentzVector> kPlus_kkpipi, kMinus_kkpipi;
  piPlus_4pi.clear();
  piMinus_4pi.clear(), pi0_p4.clear();
  kPlus_kkpipi.clear();
  kMinus_kkpipi.clear();
 
  for ( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end();
        it++ )
  {
    int pdg_code = ( *it )->particleProperty();
    if ( ( *it )->primaryParticle() ) continue;
    Event::McParticle it2        = ( *it )->mother();
    int               mother_pdg = 0;
 
    // finds if the particle is related to the d0 we need
    mother_pdg = it2.particleProperty();
 
    // special cases
    if ( pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID ||
         pdg_code == kK0BarID || pdg_code == kKDPStar0ID || pdg_code == kKDPStar0BarID ||
         pdg_code == kK10ID || pdg_code == kK10BarID || pdg_code == kK0Star0ID ||
         pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID )
      continue; // don't count special intermediate states
 
    if ( mother_pdg == kK0ID || mother_pdg == kK0BarID )
    { // Looks another step up if mother was k0 or k0bar
      it2        = it2.mother();
      mother_pdg = it2.particleProperty();
    }
 
    if ( mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID )
    { // Looks another step up if mother was k*0 or k*0bar
      // if code is found to be K*0 -> K+Pi- we save  it as a KStar0 and K*0Bar -> K-Pi+ as
      // KStar0Bar if is it a neutral decay K*0(Bar)->K0(bar) Pi0 or Gamma.  We save them as
      // K0(bar) and Pi0 or Gamma.
      if ( pdg_code == kPiPlusID || pdg_code == kPiMinusID ) continue;
      if ( mother_pdg == kKStar0ID && pdg_code == kKPlusID ) pdg_code = kKStar0ID;
      if ( mother_pdg == kKStar0BarID && pdg_code == kKMinusID ) pdg_code = kKStar0BarID;
      it2        = it2.mother();
      mother_pdg = it2.particleProperty();
    }
 
    if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID )
    { // Looks another step up if mother was K_0
      it2        = it2.mother();
      mother_pdg = it2.particleProperty();
    }
 
    if ( mother_pdg == kK10ID || mother_pdg == kK10BarID )
    { // Looks another step up if mother was k_10 or k_10bar
      it2        = it2.mother();
      mother_pdg = it2.particleProperty();
    }
 
    // Identifies what particles are the daughters
    int nFSR( 0 );
    if ( mother_pdg == d0_pdg ) // Found daughter
    {
      // if(pdg_code ==  kGammaFSRID)continue;
      // if(m_debug)cout<<"pdgid = "<<pdg_code<<endl;
      if ( pdg_code == kPiPlusID || pdg_code == kA0PlusID ) num[0]++;
      else if ( pdg_code == kPiMinusID || pdg_code == kA0MinusID ) num[1]++;
      else if ( pdg_code == kPi0ID ) num[2]++;
      else if ( pdg_code == kEtaID ) num[3]++;
      else if ( pdg_code == kEtaPrimeID ) num[4]++;
      else if ( pdg_code == kF0ID || pdg_code == kA00ID || pdg_code == kFPrime0ID ||
                pdg_code == kF0m1500ID || pdg_code == kF2ID || pdg_code == kF0m1710ID ||
                pdg_code == kF0600ID )
        num[5]++;
      else if ( pdg_code == kRhoPlusID || pdg_code == kRho2SPlusID || pdg_code == kA1PlusID )
        num[6]++;
      else if ( pdg_code == kRhoMinusID || pdg_code == kRho2SMinusID ||
                pdg_code == kA1MinusID )
        num[7]++;
      else if ( pdg_code == kRho0ID || pdg_code == kRho2S0ID ) num[8]++;
      else if ( pdg_code == kOmegaID ) num[9]++;
      else if ( pdg_code == kPhiID ) num[10]++;
      else if ( pdg_code == kKPlusID ) num[11]++;
      else if ( pdg_code == kKMinusID ) num[12]++;
      else if ( pdg_code == kK0SID ) num[13]++;
      else if ( pdg_code == kK0LID ) num[14]++;
      else if ( pdg_code == kKStarPlusID || pdg_code == kK1PlusID ||
                pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID ||
                pdg_code == kK0StarPlusID )
        num[15]++;
      else if ( pdg_code == kKStarMinusID || pdg_code == kK1MinusID ||
                pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID ||
                pdg_code == kK0StarMinusID )
        num[16]++;
      else if ( pdg_code == kKStar0ID ) num[17]++;
      else if ( pdg_code == kKStar0BarID ) num[18]++;
      else if ( pdg_code == kK10ID || pdg_code == kK1Prime0ID ) num[19]++;
      else if ( pdg_code == kK10BarID || pdg_code == kK1Prime0BarID ) num[20]++;
      else if ( pdg_code == kEPlusID || pdg_code == kMuPlusID ) num[21]++;
      else if ( pdg_code == kEMinusID || pdg_code == kMuMinusID ) num[22]++;
      else if ( pdg_code == kNuEID || pdg_code == kNuEBarID || pdg_code == kNuMuID ||
                pdg_code == kNuMuBarID )
        num[23]++;
      else if ( pdg_code == kGammaID ) num[24]++;
      else if ( pdg_code == kGammaFSRID ) continue;
      else
      {
        num[25]++;
        cout << "Unknown particle:  " << pdg_code << endl;
      }
 
      // if (pdg_code ==  kA10ID )            num[0]++; // not present
 
      // Enter Momentums for Dalitz canidates
      // It would be prefered if we could get PreFSR Momentum
      // m_D0daughter_P4.push_back( (*it)->initialFourMomentum().vect() );
      // m_D0daughter_PID.push_back( pdg_code );
 
      if ( pdg_code == kPiPlusID )
        piPlus.setVectM( ( *it )->initialFourMomentum().vect(), xmpion );
      if ( pdg_code == kPiMinusID )
        piMinus.setVectM( ( *it )->initialFourMomentum().vect(), xmpion );
      if ( pdg_code == kK0SID ) k0.setVectM( ( *it )->initialFourMomentum().vect(), xmk0 );
      if ( pdg_code == kK0LID ) k0.setVectM( ( *it )->initialFourMomentum().vect(), xmk0 );
      HepLorentzVector daughterP4;
      if ( piPlus_4pi.size() > 2 || piMinus_4pi.size() > 2 ) continue;
      if ( kPlus_kkpipi.size() > 1 || kMinus_kkpipi.size() > 1 ) continue;
      // log << MSG::INFO << "in D0->pipipipi;" << endmsg;
      if ( pdg_code == kPiPlusID )
      {
        daughterP4.setVectM( ( *it )->initialFourMomentum().vect(), xmpion );
        piPlus_4pi.push_back( daughterP4 );
      }
      if ( pdg_code == kPiMinusID )
      {
        daughterP4.setVectM( ( *it )->initialFourMomentum().vect(), xmpion );
        piMinus_4pi.push_back( daughterP4 );
      }
      if ( pdg_code == kKPlusID )
      {
        daughterP4.setVectM( ( *it )->initialFourMomentum().vect(), xmkaon );
        kPlus_kkpipi.push_back( daughterP4 );
      }
      if ( pdg_code == kKMinusID )
      {
        daughterP4.setVectM( ( *it )->initialFourMomentum().vect(), xmkaon );
        kMinus_kkpipi.push_back( daughterP4 );
      }
      if ( pdg_code == kPi0ID )
      {
        daughterP4.setVectM( ( *it )->initialFourMomentum().vect(), xmpi0 );
        pi0_p4.push_back( daughterP4 );
      }
    }
  }
 
  // D decay daughters have been tabulated.  Now, classify decay.
  int nDaughters = 0;
  for ( int i = 0; i < 26; i++ ) { nDaughters += num[i]; }
  // if(m_debug)cout<<nDaughters<<num[ kUFVPlus ]<<num[ kPiMinus ]<< nDaughters_RhoPlus<<
  // nDaughterPiPlus_RhoPlus     << nDaughterPi0_RhoPlus<<endl;
  // if(m_debug)cout<<nDaughters<<num[ kUFVMinus ]<<num[ kPiPlus ]<< nDaughters_RhoMinus<<
  // nDaughterPiMinus_RhoMinus<<  nDaughterPi0_RhoMinus<<endl;
  // if(m_debug)cout<<nDaughters<<num[ kRho0 ]<<num[ kPi0 ]<< nDaughters_Rho0       <<
  // nDaughterPiPlus_Rho0        << nDaughterPiMinus_Rho0<<endl;
  // if(m_debug)cout<<nDaughters<< nDaughters_F0600 << nDaughterPiPlus_F0600
  // << nDaughterPiMinus_F0600<<endl; if(m_debug)cout<<nDaughters<< nDaughters_F0
  // << nDaughterPiPlus_F0           << nDaughterPiMinus_F0<<endl;
  // if(m_debug)cout<<nDaughters<< nDaughters_FPrime0<< nDaughterPiPlus_FPrime0  <<
  // nDaughterPiMinus_FPrime0<<endl; if(m_debug)cout<<nDaughters<< nDaughters_F01500    <<
  // nDaughterPiPlus_F01500  << nDaughterPiMinus_F01500<<endl;
  // if(m_debug)cout<<nDaughters<< nDaughters_F01710    << nDaughterPiPlus_F01710    <<
  // nDaughterPiMinus_F01710<<endl; if(m_debug)cout<<nDaughters<< nDaughters_F2
  // << nDaughterPiPlus_F2           << nDaughterPiMinus_F2<<endl;
  // if(m_debug)cout<<nDaughters<< nDaughters_Omega << nDaughterPiPlus_Omega
  // << nDaughterPiMinus_Omega<<endl;
 
  int nUFP0        = num[kPi0] + num[kEta] + num[kEtaPrime];
  int nUFV0        = num[kRho0] + num[kPhi] + num[kOmega] + num[kGamma];
  int nFV0         = num[kKStar0] + num[kK10];
  int nFV0Bar      = num[kKStar0Bar] + num[kK10Bar];
  int nStrange     = num[kKMinus] + num[kFVMinus] + nFV0Bar;
  int nAntiStrange = num[kKPlus] + num[kFVPlus] + nFV0;
  int nCPPlusEig =
      num[kNeutralScalar] + num[kRho0] + num[kOmega] + num[kPhi] + num[kK0S] + num[kGamma];
  int nCPMinusEig = num[kPi0] + num[kEta] + num[kEtaPrime] + num[kK0L];
  int nCPEig      = nCPPlusEig + nCPMinusEig;
  int nChargedPiK = num[kPiPlus] + num[kPiMinus] + num[kKPlus] + num[kKMinus];
  int nK0         = num[kK0S] + num[kK0L];
 
  // check Dalitz modes: KsPi+Pi- or KlPi+Pi-
  double mnm_gen     = 0.;
  double mpn_gen     = 0.;
  double mpm_gen     = 0.;
  bool   isKsPiPi    = false;
  bool   isKsKK      = false;
  bool   isKsPiPiPi0 = false;
 
  // check K3Pi modes:
 
  if ( nK0 == 1 && num[kPiPlus] == 1 && num[kPiMinus] && nDaughters == 3 )
  {
    decay_mode = kDalitz;
    // k0.boost(-0.011, 0, 0);
    // piMinus.boost(-0.011, 0, 0);
    // piPlus.boost(-0.011, 0, 0);
    // mnm_gen = (k0 + piMinus).m2();
    // mpn_gen = (k0 + piPlus).m2();
    // mpm_gen = (piPlus + piMinus).m2();
    if ( num[kK0S] == 1 ) isKsPiPi = true;
  }
  if ( num[kPiPlus] == 2 && num[kPiMinus] == 2 && nDaughters == 4 ) { decay_mode = k2Pip2Pim; }
  if ( num[kPiPlus] == 1 && num[kPiMinus] == 1 && num[kPi0] == 2 && nDaughters == 4 )
  { decay_mode = kPipPim2Pi0; }
  if ( nK0 == 1 && num[kPiPlus] == 1 && num[kPiMinus] == 1 && num[kPi0] == 1 &&
       nDaughters == 4 )
  {
    decay_mode = kK0PiPiPi0;
    if ( num[kK0S] == 1 ) isKsPiPiPi0 = true;
  }
  if ( num[kPiPlus] == 1 && num[kPiMinus] == 1 && num[kKPlus] == 1 && num[kKMinus] == 1 &&
       nDaughters == 4 )
  { decay_mode = kKKPiPi; }
  if ( num[kKPlus] == 1 && num[kKMinus] == 1 && nK0 == 1 && nDaughters == 3 )
  {
    decay_mode = kK0KK;
    if ( num[kK0S] == 1 ) isKsKK = true;
  }
 
  int nDaughters_RhoPlus( 0 );
  int nDaughterPiPlus_RhoPlus( 0 );
  int nDaughterPi0_RhoPlus( 0 );
  int nDaughters_RhoMinus( 0 );
  int nDaughterPiMinus_RhoMinus( 0 );
  int nDaughterPi0_RhoMinus( 0 );
  int nDaughters_Rho0( 0 );
  int nDaughterPiPlus_Rho0( 0 );
  int nDaughterPiMinus_Rho0( 0 );
  int nDaughters_F0600( 0 );
  int nDaughterPiPlus_F0600( 0 );
  int nDaughterPiMinus_F0600( 0 );
  int nDaughters_F0( 0 );
  int nDaughterPiPlus_F0( 0 );
  int nDaughterPiMinus_F0( 0 );
  int nDaughters_FPrime0( 0 );
  int nDaughterPiPlus_FPrime0( 0 );
  int nDaughterPiMinus_FPrime0( 0 );
  int nDaughters_F01500( 0 );
  int nDaughterPiPlus_F01500( 0 );
  int nDaughterPiMinus_F01500( 0 );
  int nDaughters_F01710( 0 );
  int nDaughterPiPlus_F01710( 0 );
  int nDaughterPiMinus_F01710( 0 );
  int nDaughters_F2( 0 );
  int nDaughterPiPlus_F2( 0 );
  int nDaughterPiMinus_F2( 0 );
  int nDaughters_Omega( 0 );
  int nDaughterPiPlus_Omega( 0 );
  int nDaughterPiMinus_Omega( 0 );
 
  for ( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end();
        it++ )
  {
    int pdg_code = ( *it )->particleProperty();
    if ( ( *it )->primaryParticle() ) continue;
    Event::McParticle it2        = ( *it )->mother();
    int               mother_pdg = 0;
 
    // finds if the particle is related to the d0 we need
    mother_pdg = it2.particleProperty();
    if ( mother_pdg == kRhoPlusID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_RhoPlus++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_RhoPlus++;
        if ( pdg_code == kPi0ID ) nDaughterPi0_RhoPlus++;
      }
    }
    if ( mother_pdg == kRhoMinusID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_RhoMinus++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_RhoMinus++;
        if ( pdg_code == kPi0ID ) nDaughterPi0_RhoMinus++;
      }
    }
    if ( mother_pdg == kRho0ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_Rho0++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_Rho0++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_Rho0++;
      }
    }
    if ( mother_pdg == kF0600ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_F0600++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_F0600++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_F0600++;
      }
    }
    if ( mother_pdg == kF0ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_F0++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_F0++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_F0++;
      }
    }
    if ( mother_pdg == kFPrime0ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_FPrime0++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_FPrime0++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_FPrime0++;
      }
    }
    if ( mother_pdg == kF0m1500ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_F01500++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_F01500++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_F01500++;
      }
    }
    if ( mother_pdg == kF0m1710ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_F01710++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_F01710++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_F01710++;
      }
    }
    if ( mother_pdg == kF2ID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_F2++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_F2++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_F2++;
      }
    }
    if ( mother_pdg == kOmegaID )
    {
      Event::McParticle it3 = it2.mother();
      if ( it3.particleProperty() == d0_pdg && pdg_code != -22 )
      {
        nDaughters_Omega++;
        if ( pdg_code == kPiPlusID ) nDaughterPiPlus_Omega++;
        if ( pdg_code == kPiMinusID ) nDaughterPiMinus_Omega++;
      }
    }
  }
 
  // The order of if-statements below matters.
  if ( decay_mode == kDalitz )
  {
    ++m_nDalitz;
 
    //  Dalitz t(8) ;
    //  TComplex D0, D0bar;
    //  if (m_dalitzModel == 0) { //Cleo model
    //      D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
    //  if (m_dalitzModel == 1) { //Babar model
    //      D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
    //  if (m_dalitzModel == 2) { //Belle model
    //      D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
 
    complex<double> D0, D0bar;
    vector<double>  k0l;
    vector<double>  pip;
    vector<double>  pim;
    k0l.clear();
    pip.clear();
    pim.clear();
    k0l.push_back( k0[0] );
    k0l.push_back( k0[1] );
    k0l.push_back( k0[2] );
    k0l.push_back( k0[3] );
    pip.push_back( piPlus[0] );
    pip.push_back( piPlus[1] );
    pip.push_back( piPlus[2] );
    pip.push_back( piPlus[3] );
    pim.push_back( piMinus[0] );
    pim.push_back( piMinus[1] );
    pim.push_back( piMinus[2] );
    pim.push_back( piMinus[3] );
 
    if ( isKsPiPi )
    {
      // D0ToKSpipi2018 tKSpipi;tKSpipi.init();
      D0ToKSpipi tKSpipi;
      tKSpipi.init();
      D0 = tKSpipi.Amp_PFT( k0l, pip, pim );
 
      vector<double> cpk0l;
      vector<double> cppip;
      vector<double> cppim;
      cpk0l.clear();
      pip.clear();
      pim.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cppim.push_back( -piPlus[0] );
      cppim.push_back( -piPlus[1] );
      cppim.push_back( -piPlus[2] );
      cppim.push_back( piPlus[3] );
      cppip.push_back( -piMinus[0] );
      cppip.push_back( -piMinus[1] );
      cppip.push_back( -piMinus[2] );
      cppip.push_back( piMinus[3] );
      D0bar = tKSpipi.Amp_PFT( cpk0l, cppip, cppim );
    }
    else
    {
      // D0ToKLpipi2018 tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx
      // is not the same as used in other analyses
      D0ToKLpipi tKLpipi;
      tKLpipi.init(); // be aware of the convention in D0ToKLpipi.cxx is not the same as used
                      // in other analyses
      D0 = tKLpipi.Amp_PFT( k0l, pip, pim );
      vector<double> cpk0l;
      vector<double> cppip;
      vector<double> cppim;
      cpk0l.clear();
      pip.clear();
      pim.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cppim.push_back( -piPlus[0] );
      cppim.push_back( -piPlus[1] );
      cppim.push_back( -piPlus[2] );
      cppim.push_back( piPlus[3] );
      cppip.push_back( -piMinus[0] );
      cppip.push_back( -piMinus[1] );
      cppip.push_back( -piMinus[2] );
      cppip.push_back( piMinus[3] );
      D0bar = tKLpipi.Amp_PFT( cpk0l, cppip, cppim );
    }
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      log << MSG::INFO << "D0 calculation " << sqrt( r2D0 ) << " " << temp << endmsg;
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsPiPi ? deltaD0 : ( deltaD0 + PI );
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      if ( mpn_gen - mnm_gen > 0. )
      {
        m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd;
        m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x;
        m_dalitzDenom += -2. * r2D0 * cosd * cosd;
      }
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsPiPi ? deltaD0 : ( deltaD0 + PI );
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      if ( mpn_gen - mnm_gen < 0. )
      {
        m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd;
        m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x;
        m_dalitzDenom += -2. * r2D0 * cosd * cosd;
      }
    }
  }
  else if ( decay_mode == k2Pip2Pim )
  {
 
    ++m_n2Pip2Pim;
 
    HepLorentzVector pip1_4pi = ( piPlus_4pi.at( 0 ) );
    HepLorentzVector pim1_4pi = ( piMinus_4pi.at( 0 ) );
    HepLorentzVector pip2_4pi = ( piPlus_4pi.at( 1 ) );
    HepLorentzVector pim2_4pi = ( piMinus_4pi.at( 1 ) );
    pip1_4pi.boost( -0.011, 0, 0 );
    pip2_4pi.boost( -0.011, 0, 0 );
    pim1_4pi.boost( -0.011, 0, 0 );
    pim2_4pi.boost( -0.011, 0, 0 );
    std::complex<double> D0, D0bar;
    // EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
    // EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
    // EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
    // EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
    std::vector<double> pip1;
    pip1.clear();
    pip1.push_back( pip1_4pi[0] );
    pip1.push_back( pip1_4pi[1] );
    pip1.push_back( pip1_4pi[2] );
    pip1.push_back( pip1_4pi[3] );
    std::vector<double> pim1;
    pim1.clear();
    pim1.push_back( pim1_4pi[0] );
    pim1.push_back( pim1_4pi[1] );
    pim1.push_back( pim1_4pi[2] );
    pim1.push_back( pim1_4pi[3] );
    std::vector<double> pip2;
    pip2.clear();
    pip2.push_back( pip2_4pi[0] );
    pip2.push_back( pip2_4pi[1] );
    pip2.push_back( pip2_4pi[2] );
    pip2.push_back( pip2_4pi[3] );
    std::vector<double> pim2;
    pim2.clear();
    pim2.push_back( pim2_4pi[0] );
    pim2.push_back( pim2_4pi[1] );
    pim2.push_back( pim2_4pi[2] );
    pim2.push_back( pim2_4pi[3] );
 
    D0To2pip2pim t2pip2pim;
    t2pip2pim.init();
    D0 = t2pip2pim.Amp( pip1, pim1, pip2, pim2 );
    // EvtVector4R cppim1(pip1_4pi[3],    -1*pip1_4pi[0],   -1*pip1_4pi[1],   -1*pip1_4pi[2]);
    // EvtVector4R cppip1(pim1_4pi[3],    -1*pim1_4pi[0],   -1*pim1_4pi[1],   -1*pim1_4pi[2]);
    // EvtVector4R cppim2(pip2_4pi[3],    -1*pip2_4pi[0],   -1*pip2_4pi[1],   -1*pip2_4pi[2]);
    // EvtVector4R cppip2(pim2_4pi[3],    -1*pim2_4pi[0],   -1*pim2_4pi[1],   -1*pim2_4pi[2]);
    std::vector<double> cppip1;
    cppip1.clear();
    cppip1.push_back( -pim1_4pi[0] );
    cppip1.push_back( -pim1_4pi[1] );
    cppip1.push_back( -pim1_4pi[2] );
    cppip1.push_back( pim1_4pi[3] );
    std::vector<double> cppim1;
    cppim1.clear();
    cppim1.push_back( -pip1_4pi[0] );
    cppim1.push_back( -pip1_4pi[1] );
    cppim1.push_back( -pip1_4pi[2] );
    cppim1.push_back( pip1_4pi[3] );
    std::vector<double> cppip2;
    cppip2.clear();
    cppip2.push_back( -pim2_4pi[0] );
    cppip2.push_back( -pim2_4pi[1] );
    cppip2.push_back( -pim2_4pi[2] );
    cppip2.push_back( pim2_4pi[3] );
    std::vector<double> cppim2;
    cppim2.clear();
    cppim2.push_back( -pip2_4pi[0] );
    cppim2.push_back( -pip2_4pi[1] );
    cppim2.push_back( -pip2_4pi[2] );
    cppim2.push_back( pip2_4pi[3] );
 
    D0bar = t2pip2pim.Amp( cppip1, cppim1, cppip2, cppim2 );
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      // if( mpn_gen - mnm_gen > 0. )
      //{
      m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd;
      m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x;
      m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd;
      //}
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      // if( mpn_gen - mnm_gen < 0. )
      //{
      m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd;
      m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x;
      m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd;
      //}
    }
  }
  else if ( decay_mode == kPipPim2Pi0 )
  {
 
    ++m_nPipPim2Pi0;
 
    HepLorentzVector pip1_4pi = ( piPlus_4pi.at( 0 ) );
    HepLorentzVector pim1_4pi = ( piMinus_4pi.at( 0 ) );
    HepLorentzVector pi01_p4  = ( pi0_p4.at( 0 ) );
    HepLorentzVector pi02_p4  = ( pi0_p4.at( 1 ) );
    pip1_4pi.boost( -0.011, 0, 0 );
    pim1_4pi.boost( -0.011, 0, 0 );
    pi01_p4.boost( -0.011, 0, 0 );
    pi02_p4.boost( -0.011, 0, 0 );
    std::complex<double> D0, D0bar;
    // EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
    // EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
    // EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
    // EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
    std::vector<double> pip1;
    pip1.clear();
    pip1.push_back( pip1_4pi[0] );
    pip1.push_back( pip1_4pi[1] );
    pip1.push_back( pip1_4pi[2] );
    pip1.push_back( pip1_4pi[3] );
    std::vector<double> pim1;
    pim1.clear();
    pim1.push_back( pim1_4pi[0] );
    pim1.push_back( pim1_4pi[1] );
    pim1.push_back( pim1_4pi[2] );
    pim1.push_back( pim1_4pi[3] );
    std::vector<double> pi01;
    pi01.clear();
    pi01.push_back( pi01_p4[0] );
    pi01.push_back( pi01_p4[1] );
    pi01.push_back( pi01_p4[2] );
    pi01.push_back( pi01_p4[3] );
    std::vector<double> pi02;
    pi02.clear();
    pi02.push_back( pi02_p4[0] );
    pi02.push_back( pi02_p4[1] );
    pi02.push_back( pi02_p4[2] );
    pi02.push_back( pi02_p4[3] );
 
    D0Topippim2pi0 tpippim2pi0;
    tpippim2pi0.init();
    D0 = tpippim2pi0.Amp( pip1, pim1, pi01, pi02 );
    // EvtVector4R cppim1(pip1_4pi[3],    -1*pip1_4pi[0],   -1*pip1_4pi[1],   -1*pip1_4pi[2]);
    // EvtVector4R cppip1(pim1_4pi[3],    -1*pim1_4pi[0],   -1*pim1_4pi[1],   -1*pim1_4pi[2]);
    // EvtVector4R cppim2(pip2_4pi[3],    -1*pip2_4pi[0],   -1*pip2_4pi[1],   -1*pip2_4pi[2]);
    // EvtVector4R cppip2(pim2_4pi[3],    -1*pim2_4pi[0],   -1*pim2_4pi[1],   -1*pim2_4pi[2]);
    std::vector<double> cppip1;
    cppip1.clear();
    cppip1.push_back( -pim1_4pi[0] );
    cppip1.push_back( -pim1_4pi[1] );
    cppip1.push_back( -pim1_4pi[2] );
    cppip1.push_back( pim1_4pi[3] );
    std::vector<double> cppim1;
    cppim1.clear();
    cppim1.push_back( -pip1_4pi[0] );
    cppim1.push_back( -pip1_4pi[1] );
    cppim1.push_back( -pip1_4pi[2] );
    cppim1.push_back( pip1_4pi[3] );
    std::vector<double> cppi01;
    cppi01.clear();
    cppi01.push_back( -pi01_p4[0] );
    cppi01.push_back( -pi01_p4[1] );
    cppi01.push_back( -pi01_p4[2] );
    cppi01.push_back( pi01_p4[3] );
    std::vector<double> cppi02;
    cppi02.clear();
    cppi02.push_back( -pi02_p4[0] );
    cppi02.push_back( -pi02_p4[1] );
    cppi02.push_back( -pi02_p4[2] );
    cppi02.push_back( pi02_p4[3] );
 
    D0bar = tpippim2pi0.Amp( cppip1, cppim1, cppi01, cppi02 );
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      // if( mpn_gen - mnm_gen > 0. )
      //{
      m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd;
      m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x;
      m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd;
      //}
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
 
      m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd;
      m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x;
      m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd;
    }
  }
  else if ( decay_mode == kK0PiPiPi0 )
  {
    ++m_nK0PiPiPi0;
 
    //  Dalitz t(8) ;
    //  TComplex D0, D0bar;
    //  if (m_dalitzModel == 0) { //Cleo model
    //      D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
    //  if (m_dalitzModel == 1) { //Babar model
    //      D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
    //  if (m_dalitzModel == 2) { //Belle model
    //      D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
    //      D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
    HepLorentzVector pip_kspipipi0 = ( piPlus_4pi.at( 0 ) );
    HepLorentzVector pim_kspipipi0 = ( piMinus_4pi.at( 0 ) );
    HepLorentzVector k0_kspipipi0  = k0;
    HepLorentzVector pi0_kspipipi0 = ( pi0_p4.at( 0 ) );
    pip_kspipipi0.boost( -0.011, 0, 0 );
    pim_kspipipi0.boost( -0.011, 0, 0 );
    k0_kspipipi0.boost( -0.011, 0, 0 );
    pi0_kspipipi0.boost( -0.011, 0, 0 );
 
    complex<double> D0, D0bar;
    vector<double>  k0l;
    vector<double>  pip;
    vector<double>  pim;
    vector<double>  pi0;
    k0l.clear();
    pip.clear();
    pim.clear();
    pi0.clear();
    k0l.push_back( k0_kspipipi0[0] );
    k0l.push_back( k0_kspipipi0[1] );
    k0l.push_back( k0_kspipipi0[2] );
    k0l.push_back( k0_kspipipi0[3] );
    pip.push_back( pip_kspipipi0[0] );
    pip.push_back( pip_kspipipi0[1] );
    pip.push_back( pip_kspipipi0[2] );
    pip.push_back( pip_kspipipi0[3] );
    pim.push_back( pim_kspipipi0[0] );
    pim.push_back( pim_kspipipi0[1] );
    pim.push_back( pim_kspipipi0[2] );
    pim.push_back( pim_kspipipi0[3] );
    pi0.push_back( pi0_kspipipi0[0] );
    pi0.push_back( pi0_kspipipi0[1] );
    pi0.push_back( pi0_kspipipi0[2] );
    pi0.push_back( pi0_kspipipi0[3] );
 
    if ( isKsPiPiPi0 )
    {
      D0ToKSpipipi0 tKSpipipi0;
      tKSpipipi0.init();
      D0 = tKSpipipi0.Amp( k0l, pip, pim, pi0 );
 
      vector<double> cpk0l;
      vector<double> cppip;
      vector<double> cppim;
      vector<double> cppi0;
      cpk0l.clear();
      cppip.clear();
      cppim.clear();
      cppi0.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cppim.push_back( -pip[0] );
      cppim.push_back( -pip[1] );
      cppim.push_back( -pip[2] );
      cppim.push_back( pip[3] );
      cppip.push_back( -pim[0] );
      cppip.push_back( -pim[1] );
      cppip.push_back( -pim[2] );
      cppip.push_back( pim[3] );
      cppi0.push_back( -pi0[0] );
      cppi0.push_back( -pi0[1] );
      cppi0.push_back( -pi0[2] );
      cppi0.push_back( pi0[3] );
      D0bar = tKSpipipi0.Amp( cpk0l, cppip, cppim, cppi0 );
    }
    else
    {
      D0ToKSpipipi0 tKSpipipi0;
      tKSpipipi0.init();
      D0 = tKSpipipi0.Amp( k0l, pip, pim, pi0 );
 
      vector<double> cpk0l;
      vector<double> cppip;
      vector<double> cppim;
      vector<double> cppi0;
      cpk0l.clear();
      cppip.clear();
      cppim.clear();
      cppi0.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cppim.push_back( -pip[0] );
      cppim.push_back( -pip[1] );
      cppim.push_back( -pip[2] );
      cppim.push_back( pip[3] );
      cppip.push_back( -pim[0] );
      cppip.push_back( -pim[1] );
      cppip.push_back( -pim[2] );
      cppip.push_back( pim[3] );
      cppi0.push_back( -pi0[0] );
      cppi0.push_back( -pi0[1] );
      cppi0.push_back( -pi0[2] );
      cppi0.push_back( pi0[3] );
      D0bar = tKSpipipi0.Amp( cpk0l, cppip, cppim, cppi0 );
      // D0ToKLpipi tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx is
      // not the same as used in other analyses D0 = tKLpipi.Amp_PFT(k0l,pip,pim);
      // vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
      // cpk0l.clear();pip.clear();pim.clear();
      // cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
      // cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
      // cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
      // D0bar = tKLpipi.Amp_PFT(cpk0l, cppip, cppim);
    }
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      log << MSG::INFO << "D0 calculation " << sqrt( r2D0 ) << " " << temp << endmsg;
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsPiPiPi0 ? ( deltaD0 + PI ) : ( deltaD0 );
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      if ( mpn_gen - mnm_gen > 0. )
      {
        m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd;
        m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x;
        m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd;
      }
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsPiPiPi0 ? ( deltaD0 + PI ) : ( deltaD0 );
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
      if ( mpn_gen - mnm_gen < 0. )
      {
        m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd;
        m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x;
        m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd;
      }
    }
  }
  else if ( decay_mode == kKKPiPi )
  {
 
    ++m_nKKPiPi;
 
    complex<double> D0, D0bar;
    vector<double>  kp;
    vector<double>  km;
    vector<double>  pip;
    vector<double>  pim;
    double          pdau[16];
    pip.clear();
    pim.clear();
    kp.clear();
    km.clear();
 
    HepLorentzVector pip_kkpipi = ( piPlus_4pi.at( 0 ) );
    HepLorentzVector pim_kkpipi = ( piMinus_4pi.at( 0 ) );
    HepLorentzVector kp_kkpipi  = ( kPlus_kkpipi.at( 0 ) );
    HepLorentzVector km_kkpipi  = ( kMinus_kkpipi.at( 0 ) );
 
    pip.push_back( pip_kkpipi[0] );
    pip.push_back( pip_kkpipi[1] );
    pip.push_back( pip_kkpipi[2] );
    pip.push_back( pip_kkpipi[3] );
    pim.push_back( pim_kkpipi[0] );
    pim.push_back( pim_kkpipi[1] );
    pim.push_back( pim_kkpipi[2] );
    pim.push_back( pim_kkpipi[3] );
    kp.push_back( kp_kkpipi[0] );
    kp.push_back( kp_kkpipi[1] );
    kp.push_back( kp_kkpipi[2] );
    kp.push_back( kp_kkpipi[3] );
    km.push_back( km_kkpipi[0] );
    km.push_back( km_kkpipi[1] );
    km.push_back( km_kkpipi[2] );
    km.push_back( km_kkpipi[3] );
 
    pdau[0]  = kp_kkpipi[0];
    pdau[1]  = kp_kkpipi[1];
    pdau[2]  = kp_kkpipi[2];
    pdau[3]  = kp_kkpipi[3];
    pdau[4]  = km_kkpipi[0];
    pdau[5]  = km_kkpipi[1];
    pdau[6]  = km_kkpipi[2];
    pdau[7]  = km_kkpipi[3];
    pdau[8]  = pip_kkpipi[0];
    pdau[9]  = pip_kkpipi[1];
    pdau[10] = pip_kkpipi[2];
    pdau[11] = pip_kkpipi[3];
    pdau[12] = pim_kkpipi[0];
    pdau[13] = pim_kkpipi[1];
    pdau[14] = pim_kkpipi[2];
    pdau[15] = pim_kkpipi[3];
 
    D0ToKKpipi tKKpipi;
    tKKpipi.init();
    D0 = tKKpipi.AMP( pdau, 1 );
 
    pdau[0]  = -km_kkpipi[0];
    pdau[1]  = -km_kkpipi[1];
    pdau[2]  = -km_kkpipi[2];
    pdau[3]  = km_kkpipi[3];
    pdau[4]  = -kp_kkpipi[0];
    pdau[5]  = -kp_kkpipi[1];
    pdau[6]  = -kp_kkpipi[2];
    pdau[7]  = kp_kkpipi[3];
    pdau[8]  = -pim_kkpipi[0];
    pdau[9]  = -pim_kkpipi[1];
    pdau[10] = -pim_kkpipi[2];
    pdau[11] = pim_kkpipi[3];
    pdau[12] = -pip_kkpipi[0];
    pdau[13] = -pip_kkpipi[1];
    pdau[14] = -pip_kkpipi[2];
    pdau[15] = pip_kkpipi[3];
    D0bar    = tKKpipi.AMP( pdau, 1 );
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      log << MSG::INFO << "D0 calculation " << sqrt( r2D0 ) << " " << temp << endmsg;
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
    }
  }
  else if ( decay_mode == kK0KK )
  {
    ++m_nK0KK;
 
    HepLorentzVector kPlus  = ( kPlus_kkpipi.at( 0 ) );
    HepLorentzVector kMinus = ( kMinus_kkpipi.at( 0 ) );
    complex<double>  D0, D0bar;
    vector<double>   k0l;
    vector<double>   kp;
    vector<double>   km;
    k0l.clear();
    kp.clear();
    km.clear();
    k0l.push_back( k0[0] );
    k0l.push_back( k0[1] );
    k0l.push_back( k0[2] );
    k0l.push_back( k0[3] );
    kp.push_back( kPlus[0] );
    kp.push_back( kPlus[1] );
    kp.push_back( kPlus[2] );
    kp.push_back( kPlus[3] );
    km.push_back( kMinus[0] );
    km.push_back( kMinus[1] );
    km.push_back( kMinus[2] );
    km.push_back( kMinus[3] );
 
    if ( isKsKK )
    {
      // D0ToKSKK tKSKK;tKSKK.init();
      // D0 = tKSKK.Amp_PFT(k0l,kp,km);
      D0ToKSLKK tKSLKK;
      tKSLKK.init( 310, 0 );
      D0 = tKSLKK.Amp( k0l, kp, km, 310, 0 );
 
      vector<double> cpk0l;
      vector<double> cpkp;
      vector<double> cpkm;
      cpk0l.clear();
      kp.clear();
      km.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cpkm.push_back( -kPlus[0] );
      cpkm.push_back( -kPlus[1] );
      cpkm.push_back( -kPlus[2] );
      cpkm.push_back( kPlus[3] );
      cpkp.push_back( -kMinus[0] );
      cpkp.push_back( -kMinus[1] );
      cpkp.push_back( -kMinus[2] );
      cpkp.push_back( kMinus[3] );
      // D0bar = tKSKK.Amp_PFT(cpk0l, cpkp, cpkm);
      D0bar = tKSLKK.Amp( cpk0l, cpkp, cpkm, 310, 0 );
    }
    else
    {
      // D0ToKSKK tKSKK;tKSKK.init();
      // D0 = tKSKK.Amp_PFT(k0l,kp,km);
      D0ToKSLKK tKSLKK;
      tKSLKK.init( 130, 1 );
      D0 = tKSLKK.Amp( k0l, kp, km, 130, 1 );
      vector<double> cpk0l;
      vector<double> cpkp;
      vector<double> cpkm;
      cpk0l.clear();
      kp.clear();
      km.clear();
      cpk0l.push_back( -k0l[0] );
      cpk0l.push_back( -k0l[1] );
      cpk0l.push_back( -k0l[2] );
      cpk0l.push_back( k0l[3] );
      cpkm.push_back( -kPlus[0] );
      cpkm.push_back( -kPlus[1] );
      cpkm.push_back( -kPlus[2] );
      cpkm.push_back( kPlus[3] );
      cpkp.push_back( -kMinus[0] );
      cpkp.push_back( -kMinus[1] );
      cpkp.push_back( -kMinus[2] );
      cpkp.push_back( kMinus[3] );
      // D0bar = tKSKK.Amp_PFT(cpk0l, cpkp, cpkm);
      D0bar = tKSLKK.Amp( cpk0l, cpkp, cpkm, 130, 1 );
    }
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      log << MSG::INFO << "D0 calculation " << sqrt( r2D0 ) << " " << temp << endmsg;
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsKK ? deltaD0 : ( deltaD0 + PI );
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
      deltaD0 = isKsKK ? deltaD0 : ( deltaD0 + PI );
    }
  }
  else if ( num[kLeptonPlus] == 1 )
  {
    decay_mode = kSLPlus;
    ++m_nSL;
 
    r2D0    = 0.;
    deltaD0 = 0.;
  }
  else if ( num[kLeptonMinus] == 1 )
  {
    decay_mode = kSLMinus;
    ++m_nSL;
 
    r2D0    = 0.;
    deltaD0 = 0.;
  }
  // check pipipi0 mode
  else if ( nDaughters == 3 && num[kPiPlus] == 1 && num[kPiMinus] == 1 && num[kPi0] == 1 )
  {
 
    if ( m_debug ) cout << "PiPiPi0 mode" << endl;
    decay_mode = kPipPimPi0;
    ++m_nPiPiPi0;
 
    HepLorentzVector pip1_4pi = ( piPlus_4pi.at( 0 ) );
    HepLorentzVector pim1_4pi = ( piMinus_4pi.at( 0 ) );
    HepLorentzVector pi01_p4  = ( pi0_p4.at( 0 ) );
    pip1_4pi.boost( -0.011, 0, 0 );
    pim1_4pi.boost( -0.011, 0, 0 );
    pi01_p4.boost( -0.011, 0, 0 );
    std::complex<double> D0, D0bar;
    // EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
    // EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
    // EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
    // EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
    std::vector<double> pip1;
    pip1.clear();
    pip1.push_back( pip1_4pi[0] );
    pip1.push_back( pip1_4pi[1] );
    pip1.push_back( pip1_4pi[2] );
    pip1.push_back( pip1_4pi[3] );
    std::vector<double> pim1;
    pim1.clear();
    pim1.push_back( pim1_4pi[0] );
    pim1.push_back( pim1_4pi[1] );
    pim1.push_back( pim1_4pi[2] );
    pim1.push_back( pim1_4pi[3] );
    std::vector<double> pi01;
    pi01.clear();
    pi01.push_back( pi01_p4[0] );
    pi01.push_back( pi01_p4[1] );
    pi01.push_back( pi01_p4[2] );
    pi01.push_back( pi01_p4[3] );
 
    D0Topipipi0 tpipipi0;
    tpipipi0.init();
    D0 = tpipipi0.Amp( pip1, pim1, pi01 );
    // EvtVector4R cppim1(pip1_4pi[3],    -1*pip1_4pi[0],   -1*pip1_4pi[1],   -1*pip1_4pi[2]);
    // EvtVector4R cppip1(pim1_4pi[3],    -1*pim1_4pi[0],   -1*pim1_4pi[1],   -1*pim1_4pi[2]);
    // EvtVector4R cppim2(pip2_4pi[3],    -1*pip2_4pi[0],   -1*pip2_4pi[1],   -1*pip2_4pi[2]);
    // EvtVector4R cppip2(pim2_4pi[3],    -1*pim2_4pi[0],   -1*pim2_4pi[1],   -1*pim2_4pi[2]);
    std::vector<double> cppip1;
    cppip1.clear();
    cppip1.push_back( -pim1_4pi[0] );
    cppip1.push_back( -pim1_4pi[1] );
    cppip1.push_back( -pim1_4pi[2] );
    cppip1.push_back( pim1_4pi[3] );
    std::vector<double> cppim1;
    cppim1.clear();
    cppim1.push_back( -pip1_4pi[0] );
    cppim1.push_back( -pip1_4pi[1] );
    cppim1.push_back( -pip1_4pi[2] );
    cppim1.push_back( pip1_4pi[3] );
    std::vector<double> cppi01;
    cppi01.clear();
    cppi01.push_back( -pi01_p4[0] );
    cppi01.push_back( -pi01_p4[1] );
    cppi01.push_back( -pi01_p4[2] );
    cppi01.push_back( pi01_p4[3] );
 
    D0bar = ( -1.0 ) * tpipipi0.Amp( cppip1, cppim1, cppi01 );
 
    if ( charm == 1 )
    {
 
      r2D0        = std::norm( D0bar ) / std::norm( D0 );
      double temp = std::arg( D0 ) - std::arg( D0bar );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
    }
    else
    {
      r2D0        = std::norm( D0 ) / std::norm( D0bar );
      double temp = std::arg( D0bar ) - std::arg( D0 );
      while ( temp < -TMath::Pi() ) { temp += 2.0 * TMath::Pi(); }
      while ( temp > TMath::Pi() ) { temp -= 2.0 * TMath::Pi(); }
      deltaD0 = temp;
 
      double cosd = cos( deltaD0 );
      double sind = sin( deltaD0 );
    }
    if ( m_debug ) cout << "deltaD0=" << deltaD0 << endl;
    if ( m_debug ) cout << "r2D0=" << r2D0 << endl;
  }
  // Check CP even modes
  else if ( ( nDaughters == 2 &&
              ( ( num[kPiPlus] == 1 && num[kPiMinus] == 1 ) || nUFP0 == 2 ||
                num[kNeutralScalar] == 2 || ( nUFP0 == 1 && nUFV0 == 1 ) ||
                ( num[kKPlus] == 1 && num[kKMinus] == 1 ) || num[kK0S] == 2 ||
                num[kK0L] == 2 || ( num[kK0L] == 1 && nUFP0 == 1 ) ||
                ( num[kK0S] == 1 && num[kNeutralScalar] == 1 ) ||
                ( num[kK0L] == 1 && nUFV0 == 1 ) ||
                ( num[kUFVPlus] == 1 && num[kPiMinus] == 1 ) ||
                ( num[kUFVMinus] == 1 && num[kPiPlus] == 1 ) ) ) ||
            ( nDaughters == 3 &&
              ( ( nCPPlusEig == 1 && num[kPi0] == 2 ) || ( num[kK0S] == 3 ) ||
                ( num[kK0S] == 1 && num[kK0L] == 2 ) ||
                ( num[kK0S] == 1 && num[kK0L] == 1 && nUFP0 == 1 ) ||
                ( num[kK0S] == 1 && nUFP0 == 2 ) ) ) ||
            ( nDaughters == 4 && ( num[kK0L] == 1 && nUFP0 == 3 ) ||
              ( ( num[kK0S] == 2 || num[kK0L] == 2 ) && nUFP0 == 2 ) ) )
  {
    decay_mode = kCPPlus;
    ++m_nCPPlus;
 
    r2D0    = 1.;
    deltaD0 = 0.;
  }
  // Check CP odd modes
  else if ( ( nDaughters == 2 && ( ( num[kK0S] == 1 && nUFP0 == 1 ) ||
                                   ( num[kK0L] == 1 && num[kNeutralScalar] == 1 ) ||
                                   ( num[kK0S] == 1 && nUFV0 == 1 ) ||
                                   ( num[kNeutralScalar] == 1 && nUFV0 == 1 ) ||
                                   ( nUFP0 == 1 && num[kNeutralScalar] == 1 ) ) ) ||
            ( nDaughters == 3 && ( ( nCPMinusEig == 3 && num[kPi0] == 2 ) || // pi0 is subset
                                                                             // of CP-
                                   ( num[kPi0] == 3 ) || ( num[kK0L] == 3 ) ||
                                   ( num[kK0S] == 2 && num[kK0L] == 1 ) ||
                                   ( ( num[kK0S] == 2 || num[kK0L] == 2 ) && nUFP0 == 1 ) ||
                                   ( num[kK0L] == 1 && nUFP0 == 2 ) ) ) ||
            ( nDaughters == 4 && ( num[kK0S] == 1 && num[kPi0] == 3 ) ||
              ( ( num[kK0S] == 1 && num[kK0L] == 1 ) && nUFP0 == 2 ) ) )
  {
    decay_mode = kCPMinus;
    ++m_nCPMinus;
 
    r2D0    = 1.;
    deltaD0 = PI;
  }
  else if ( nStrange == nAntiStrange + 1 )
  {
    if ( m_debug ) cout << nDaughters << endl;
    if ( m_debug ) cout << num[kKMinus] << endl;
    if ( m_debug ) cout << num[kPiMinus] << endl;
    if ( m_debug ) cout << num[kPiPlus] << endl;
    if ( charm == 1 )
    {
      if ( ( num[kKMinus] == 1 && num[kPiPlus] == 2 && num[kPiMinus] == 1 ) &&
           nDaughters == 4 )
      {
        if ( m_debug ) cout << "True K-3Pi " << endl;
        decay_mode = kFlavoredCF_3;
        ++m_nFlavoredCFD0_3;
 
        complex<double> D0, D0bar;
        vector<double>  kp;
        vector<double>  pi1;
        vector<double>  pi2;
        vector<double>  pi3;
        double          pdau[16];
        pi1.clear();
        pi2.clear();
        pi3.clear();
        kp.clear();
 
        HepLorentzVector kp_kpipipi  = ( kMinus_kkpipi.at( 0 ) );
        HepLorentzVector pi1_kpipipi = ( piPlus_4pi.at( 0 ) );
        HepLorentzVector pi2_kpipipi = ( piPlus_4pi.at( 1 ) );
        HepLorentzVector pi3_kpipipi = ( piMinus_4pi.at( 0 ) );
 
        kp.push_back( kp_kpipipi[0] );
        kp.push_back( kp_kpipipi[1] );
        kp.push_back( kp_kpipipi[2] );
        kp.push_back( kp_kpipipi[3] ); // k-
        pi1.push_back( pi1_kpipipi[0] );
        pi1.push_back( pi1_kpipipi[1] );
        pi1.push_back( pi1_kpipipi[2] );
        pi1.push_back( pi1_kpipipi[3] ); // pi+
        pi2.push_back( pi2_kpipipi[0] );
        pi2.push_back( pi2_kpipipi[1] );
        pi2.push_back( pi2_kpipipi[2] );
        pi2.push_back( pi2_kpipipi[3] ); // pi+
        pi3.push_back( pi3_kpipipi[0] );
        pi3.push_back( pi3_kpipipi[1] );
        pi3.push_back( pi3_kpipipi[2] );
        pi3.push_back( pi3_kpipipi[3] ); // pi-
 
        pdau[0]  = kp[0];
        pdau[1]  = kp[1];
        pdau[2]  = kp[2];
        pdau[3]  = kp[3];
        pdau[4]  = pi1[0];
        pdau[5]  = pi1[1];
        pdau[6]  = pi1[2];
        pdau[7]  = pi1[3];
        pdau[8]  = pi2[0];
        pdau[9]  = pi2[1];
        pdau[10] = pi2[2];
        pdau[11] = pi2[3];
        pdau[12] = pi3[0];
        pdau[13] = pi3[1];
        pdau[14] = pi3[2];
        pdau[15] = pi3[3];
 
        D0ToKpipipiLHCb tKpipipi;
        tKpipipi.init();
        D0 = tKpipipi.AMP( pdau, 1 );
 
        D0TopiKpipi tpiKpipi;
        tpiKpipi.init();
        std::complex<double> dcs_offset =
            0.0601387 * exp( std::complex<double>( 0, 1 ) * 1.04827 * M_PI / 180. );
        D0bar   = dcs_offset * tpiKpipi.AMP( pdau, 1 );
        deltaD0 = ( std::arg( D0 * std::conj( D0bar ) ) + m_dCF_3 );
        r2D0    = std::norm( D0bar ) / std::norm( D0 );
      }
      else if ( ( num[kKMinus] == 1 && num[kPiPlus] == 1 && num[kPi0] == 1 ) &&
                nDaughters == 3 )
      {
        decay_mode = kFlavoredCF_1;
        ++m_nFlavoredCFD0_1;
        r2D0    = pow( m_rCF_1, 2 );
        deltaD0 = m_deltaCF_1;
        RD0     = m_RCF_1;
      }
      else
      {
        // if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
        if ( m_debug ) cout << "True K-Pi " << endl;
        decay_mode = kFlavoredCF_0;
        ++m_nFlavoredCFD0_0;
        r2D0    = pow( m_rCF_0, 2 );
        deltaD0 = m_deltaCF_0;
        RD0     = m_RCF_0;
      }
    }
    else
    {
      if ( ( num[kKMinus] == 1 && num[kPiPlus] == 2 && num[kPiMinus] == 1 ) &&
           nDaughters == 4 )
      {
        decay_mode = kFlavoredCF_3;
        ++m_nFlavoredDCSD0_3;
 
        complex<double> D0, D0bar;
        vector<double>  kp;
        vector<double>  pi1;
        vector<double>  pi2;
        vector<double>  pi3;
        double          pdau[16];
        pi1.clear();
        pi2.clear();
        pi3.clear();
        kp.clear();
 
        HepLorentzVector kp_kpipipi  = ( kMinus_kkpipi.at( 0 ) );
        HepLorentzVector pi1_kpipipi = ( piPlus_4pi.at( 0 ) );
        HepLorentzVector pi2_kpipipi = ( piPlus_4pi.at( 1 ) );
        HepLorentzVector pi3_kpipipi = ( piMinus_4pi.at( 0 ) );
 
        kp.push_back( kp_kpipipi[0] );
        kp.push_back( kp_kpipipi[1] );
        kp.push_back( kp_kpipipi[2] );
        kp.push_back( kp_kpipipi[3] ); // k-
        pi1.push_back( pi1_kpipipi[0] );
        pi1.push_back( pi1_kpipipi[1] );
        pi1.push_back( pi1_kpipipi[2] );
        pi1.push_back( pi1_kpipipi[3] ); // pi+
        pi2.push_back( pi2_kpipipi[0] );
        pi2.push_back( pi2_kpipipi[1] );
        pi2.push_back( pi2_kpipipi[2] );
        pi2.push_back( pi2_kpipipi[3] ); // pi+
        pi3.push_back( pi3_kpipipi[0] );
        pi3.push_back( pi3_kpipipi[1] );
        pi3.push_back( pi3_kpipipi[2] );
        pi3.push_back( pi3_kpipipi[3] ); // pi-
 
        pdau[0]  = kp[0];
        pdau[1]  = kp[1];
        pdau[2]  = kp[2];
        pdau[3]  = kp[3];
        pdau[4]  = pi1[0];
        pdau[5]  = pi1[1];
        pdau[6]  = pi1[2];
        pdau[7]  = pi1[3];
        pdau[8]  = pi2[0];
        pdau[9]  = pi2[1];
        pdau[10] = pi2[2];
        pdau[11] = pi2[3];
        pdau[12] = pi3[0];
        pdau[13] = pi3[1];
        pdau[14] = pi3[2];
        pdau[15] = pi3[3];
 
        D0ToKpipipiLHCb tKpipipi;
        tKpipipi.init();
        D0 = tKpipipi.AMP( pdau, 1 );
 
        D0TopiKpipi tpiKpipi;
        tpiKpipi.init();
        std::complex<double> dcs_offset =
            0.0601387 * exp( std::complex<double>( 0, 1 ) * 1.04827 * M_PI / 180. );
        D0bar   = dcs_offset * tpiKpipi.AMP( pdau, 1 );
        deltaD0 = -( std::arg( D0 * std::conj( D0bar ) ) + m_dCF_3 );
        r2D0    = std::norm( D0 ) / std::norm( D0bar );
      }
      else if ( ( num[kKMinus] == 1 && num[kPiPlus] == 1 && num[kPi0] == 1 ) &&
                nDaughters == 3 )
      {
        decay_mode = kFlavoredCF_1;
        ++m_nFlavoredDCSD0_1;
        r2D0    = 1. / pow( m_rCF_1, 2 );
        deltaD0 = -m_deltaCF_1;
        RD0     = m_RCF_1;
      }
      else
      {
        // if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 ) && nDaughters == 2 )
        decay_mode = kFlavoredCF_0;
        ++m_nFlavoredDCSD0_0;
        r2D0    = 1. / pow( m_rCF_0, 2 );
        deltaD0 = -m_deltaCF_0;
        RD0     = m_RCF_0;
      }
    }
  }
  else if ( nAntiStrange == nStrange + 1 )
  {
    if ( charm == -1 )
    {
      if ( ( num[kKPlus] == 1 && num[kPiMinus] == 2 && num[kPiPlus] == 1 ) && nDaughters == 4 )
      {
        decay_mode = kFlavoredCFBar_3;
        ++m_nFlavoredCFD0_3;
 
        complex<double> D0, D0bar;
        vector<double>  kp;
        vector<double>  pi1;
        vector<double>  pi2;
        vector<double>  pi3;
        double          pdau[16];
        pi1.clear();
        pi2.clear();
        pi3.clear();
        kp.clear();
 
        HepLorentzVector kp_kpipipi  = ( kPlus_kkpipi.at( 0 ) );
        HepLorentzVector pi1_kpipipi = ( piMinus_4pi.at( 0 ) );
        HepLorentzVector pi2_kpipipi = ( piMinus_4pi.at( 1 ) );
        HepLorentzVector pi3_kpipipi = ( piPlus_4pi.at( 0 ) );
 
        kp.push_back( -kp_kpipipi[0] );
        kp.push_back( -kp_kpipipi[1] );
        kp.push_back( -kp_kpipipi[2] );
        kp.push_back( kp_kpipipi[3] ); // k-
        pi1.push_back( -pi1_kpipipi[0] );
        pi1.push_back( -pi1_kpipipi[1] );
        pi1.push_back( -pi1_kpipipi[2] );
        pi1.push_back( pi1_kpipipi[3] ); // pi+
        pi2.push_back( -pi2_kpipipi[0] );
        pi2.push_back( -pi2_kpipipi[1] );
        pi2.push_back( -pi2_kpipipi[2] );
        pi2.push_back( pi2_kpipipi[3] ); // pi+
        pi3.push_back( -pi3_kpipipi[0] );
        pi3.push_back( -pi3_kpipipi[1] );
        pi3.push_back( -pi3_kpipipi[2] );
        pi3.push_back( pi3_kpipipi[3] ); // pi-
 
        pdau[0]  = kp[0];
        pdau[1]  = kp[1];
        pdau[2]  = kp[2];
        pdau[3]  = kp[3];
        pdau[4]  = pi1[0];
        pdau[5]  = pi1[1];
        pdau[6]  = pi1[2];
        pdau[7]  = pi1[3];
        pdau[8]  = pi2[0];
        pdau[9]  = pi2[1];
        pdau[10] = pi2[2];
        pdau[11] = pi2[3];
        pdau[12] = pi3[0];
        pdau[13] = pi3[1];
        pdau[14] = pi3[2];
        pdau[15] = pi3[3];
 
        D0ToKpipipiLHCb tKpipipi;
        tKpipipi.init();
        D0 = tKpipipi.AMP( pdau, 1 );
 
        D0TopiKpipi tpiKpipi;
        tpiKpipi.init();
        std::complex<double> dcs_offset =
            0.0601387 * exp( std::complex<double>( 0, 1 ) * 1.04827 * M_PI / 180. );
        D0bar   = dcs_offset * tpiKpipi.AMP( pdau, 1 );
        deltaD0 = ( std::arg( D0 * std::conj( D0bar ) ) + m_dCF_3 );
        r2D0    = std::norm( D0bar ) / std::norm( D0 );
      }
      else if ( ( num[kKPlus] == 1 && num[kPiMinus] == 1 && num[kPi0] == 1 ) &&
                nDaughters == 3 )
      {
        decay_mode = kFlavoredCFBar_1;
        ++m_nFlavoredCFD0_1;
        r2D0    = pow( m_rCF_1, 2 );
        deltaD0 = m_deltaCF_1;
        RD0     = m_RCF_1;
      }
      else
      {
        // else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
        decay_mode = kFlavoredCFBar_0;
        ++m_nFlavoredCFD0_0;
        r2D0    = pow( m_rCF_0, 2 );
        deltaD0 = m_deltaCF_0;
        RD0     = m_RCF_0;
      }
    }
    else
    {
      if ( ( num[kKPlus] == 1 && num[kPiMinus] == 2 && num[kPiPlus] == 1 ) && nDaughters == 4 )
      {
        decay_mode = kFlavoredCFBar_3;
        ++m_nFlavoredDCSD0_3;
 
        complex<double> D0, D0bar;
        vector<double>  kp;
        vector<double>  pi1;
        vector<double>  pi2;
        vector<double>  pi3;
        double          pdau[16];
        pi1.clear();
        pi2.clear();
        pi3.clear();
        kp.clear();
 
        HepLorentzVector kp_kpipipi  = ( kPlus_kkpipi.at( 0 ) );
        HepLorentzVector pi1_kpipipi = ( piMinus_4pi.at( 0 ) );
        HepLorentzVector pi2_kpipipi = ( piMinus_4pi.at( 1 ) );
        HepLorentzVector pi3_kpipipi = ( piPlus_4pi.at( 0 ) );
 
        kp.push_back( -kp_kpipipi[0] );
        kp.push_back( -kp_kpipipi[1] );
        kp.push_back( -kp_kpipipi[2] );
        kp.push_back( kp_kpipipi[3] ); // k-
        pi1.push_back( -pi1_kpipipi[0] );
        pi1.push_back( -pi1_kpipipi[1] );
        pi1.push_back( -pi1_kpipipi[2] );
        pi1.push_back( pi1_kpipipi[3] ); // pi+
        pi2.push_back( -pi2_kpipipi[0] );
        pi2.push_back( -pi2_kpipipi[1] );
        pi2.push_back( -pi2_kpipipi[2] );
        pi2.push_back( pi2_kpipipi[3] ); // pi+
        pi3.push_back( -pi3_kpipipi[0] );
        pi3.push_back( -pi3_kpipipi[1] );
        pi3.push_back( -pi3_kpipipi[2] );
        pi3.push_back( pi3_kpipipi[3] ); // pi-
 
        pdau[0]  = kp[0];
        pdau[1]  = kp[1];
        pdau[2]  = kp[2];
        pdau[3]  = kp[3];
        pdau[4]  = pi1[0];
        pdau[5]  = pi1[1];
        pdau[6]  = pi1[2];
        pdau[7]  = pi1[3];
        pdau[8]  = pi2[0];
        pdau[9]  = pi2[1];
        pdau[10] = pi2[2];
        pdau[11] = pi2[3];
        pdau[12] = pi3[0];
        pdau[13] = pi3[1];
        pdau[14] = pi3[2];
        pdau[15] = pi3[3];
 
        D0ToKpipipiLHCb tKpipipi;
        tKpipipi.init();
        D0 = tKpipipi.AMP( pdau, 1 );
 
        D0TopiKpipi tpiKpipi;
        tpiKpipi.init();
        std::complex<double> dcs_offset =
            0.0601387 * exp( std::complex<double>( 0, 1 ) * 1.04827 * M_PI / 180. );
        D0bar   = dcs_offset * tpiKpipi.AMP( pdau, 1 );
        deltaD0 = -( std::arg( D0 * std::conj( D0bar ) ) + m_dCF_3 );
        r2D0    = std::norm( D0 ) / std::norm( D0bar );
      }
      else if ( ( num[kKPlus] == 1 && num[kPiMinus] == 1 && num[kPi0] == 1 ) &&
                nDaughters == 3 )
      {
        decay_mode = kFlavoredCFBar_1;
        ++m_nFlavoredDCSD0_1;
        r2D0    = 1. / pow( m_rCF_1, 2 );
        deltaD0 = -m_deltaCF_1;
        RD0     = m_RCF_1;
      }
      else
      {
        // if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
        decay_mode = kFlavoredCFBar_0;
        ++m_nFlavoredDCSD0_0;
        r2D0    = 1. / pow( m_rCF_0, 2 );
        deltaD0 = -m_deltaCF_0;
        RD0     = m_RCF_0;
      }
    }
  }
  else if ( nStrange == nAntiStrange ) // Unflavored or Cabibbo-suppressed
  {
    if ( ( num[kKPlus] > 0 && ( num[kKPlus] == num[kFVMinus] ||
                                num[kKPlus] == nFV0Bar ) ) || // for anti-K*0 K+ pi-
         ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0Bar ) ||
         ( num[kPiPlus] > 0 && num[kPiPlus] == num[kUFVMinus] ) )
    {
      decay_mode = kFlavoredCS;
      ++m_nFlavoredCSD0;
 
      if ( charm == 1 )
      {
        r2D0    = pow( m_rCS, 2 );
        deltaD0 = m_deltaCS;
      }
      else
      {
        r2D0    = 1. / pow( m_rCS, 2 );
        deltaD0 = -m_deltaCS;
      }
    }
    else if ( ( num[kKMinus] > 0 &&
                ( num[kKMinus] == num[kFVPlus] || num[kKMinus] == nFV0 ) ) || // for K*0 K- pi+
              ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0 ) ||
              ( num[kPiMinus] > 0 && num[kPiMinus] == num[kUFVPlus] ) )
    {
      decay_mode = kFlavoredCSBar;
      ++m_nFlavoredCSD0;
 
      if ( charm == -1 )
      {
        r2D0    = pow( m_rCS, 2 );
        deltaD0 = m_deltaCS;
      }
      else
      {
        r2D0    = 1. / pow( m_rCS, 2 );
        deltaD0 = -m_deltaCS;
      }
    }
 
    // Self-conjugate mixed-CP final states -- pick CP or non-CP at
    // random.  If CP, pick + or - at random.  If non-CP, pick
    // flavored or flavoredBar according to the appropriate value of r.
    else if ( ( nDaughters >= 3 && num[kPiPlus] == num[kPiMinus] &&
                num[kKPlus] == num[kKMinus] && nChargedPiK + nCPEig == nDaughters ) ||
              nUFV0 == nDaughters || ( num[kKStar0] > 0 && num[kKStar0] == num[kKStar0Bar] ) ||
              ( num[kK10] > 0 && num[kK10] == num[kK10Bar] ) ||
              ( num[kUFVPlus] == 1 && num[kUFVMinus] == 1 )
              // for rho+rho-; no a1+a1- and rhoa1 final states in DECAY.DEC
    )
    {
      log << MSG::DEBUG << "    [ Self-conjugate mixed-CP ]" << endmsg;
 
      if ( RandFlat::shoot() > 0.5 )
      {
        if ( RandFlat::shoot() > 0.5 )
        {
          decay_mode = kCPPlus;
          ++m_nCPPlus;
 
          r2D0    = 1.;
          deltaD0 = PI;
        }
        else
        {
          decay_mode = kCPMinus;
          ++m_nCPMinus;
 
          r2D0    = 1.;
          deltaD0 = 0.;
        }
      }
      else
      {
        // Odd # of K0S or K0L = CF; even # of K0S or K0L = SC.
        if ( nK0 % 2 )
        {
          // Nov 2007: correct for r2 -> RWS and BR != A^2
          double cutoff = m_rwsCF_0 / ( 1. + m_rwsCF_0 );
 
          if ( charm == -1 ) { cutoff = 1. - cutoff; }
 
          log << MSG::DEBUG << "    [ cutoff = " << cutoff << " ]" << endmsg;
 
          decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF_0 : kFlavoredCFBar_0;
 
          if ( ( decay_mode == kFlavoredCF_0 && charm == 1 ) ||
               ( decay_mode == kFlavoredCFBar_0 && charm == -1 ) )
          {
            ++m_nFlavoredCFD0_0;
 
            r2D0    = sqrt( m_rCF_0 );
            deltaD0 = m_deltaCF_0;
          }
          else
          {
            ++m_nFlavoredDCSD0_0;
 
            r2D0    = 1. / sqrt( m_rCF_0 );
            deltaD0 = -m_deltaCF_0;
          }
        }
        else
        {
          // Nov 2007: correct for r2 -> RWS and BR != A^2
          double cutoff = m_rwsCS / ( 1. + m_rwsCS );
 
          if ( charm == -1 ) { cutoff = 1. - cutoff; }
 
          log << MSG::DEBUG << "    [ cutoff = " << cutoff << " ]" << endmsg;
 
          decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCS : kFlavoredCSBar;
          ++m_nFlavoredCSD0;
 
          if ( ( decay_mode == kFlavoredCS && charm == 1 ) ||
               ( decay_mode == kFlavoredCSBar && charm == -1 ) )
          {
            r2D0    = sqrt( m_rCS );
            deltaD0 = m_deltaCS;
          }
          else
          {
            r2D0    = 1. / sqrt( m_rCS );
            deltaD0 = -m_deltaCS;
          }
        }
      }
    }
  }
 
  if ( num[kUnknown] >= 1 )
  {
    cout << "decay mode " << decay_mode << endl;
    cout << "D #" << charm << endl;
    cout << "pi+:     " << num[kPiPlus] << endl;
    cout << "pi-:     " << num[kPiMinus] << endl;
    cout << "pi0:     " << num[kPi0] << endl;
    cout << "eta:     " << num[kEta] << endl;
    cout << "eta':    " << num[kEtaPrime] << endl;
    cout << "f0/a0:   " << num[kNeutralScalar] << endl;
    cout << "UFV+:    " << num[kUFVPlus] << endl;
    cout << "UFV-:    " << num[kUFVMinus] << endl;
    cout << "rho0:    " << num[kRho0] << endl;
    cout << "omega:   " << num[kOmega] << endl;
    cout << "phi:     " << num[kPhi] << endl;
    cout << "K+:      " << num[kKPlus] << endl;
    cout << "K-:      " << num[kKMinus] << endl;
    cout << "K0S:     " << num[kK0S] << endl;
    cout << "K0L:     " << num[kK0L] << endl;
    cout << "FV+:     " << num[kFVPlus] << endl;
    cout << "FV-:     " << num[kFVMinus] << endl;
    cout << "K*0:     " << num[kKStar0] << endl;
    cout << "K*0b:    " << num[kKStar0Bar] << endl;
    cout << "K10:     " << num[kK10] << endl;
    cout << "K10b:    " << num[kK10Bar] << endl;
    cout << "l+:      " << num[kLeptonPlus] << endl;
    cout << "l-:      " << num[kLeptonMinus] << endl;
    cout << "nu:      " << num[kNeutrino] << endl;
    cout << "gamma:   " << num[kGamma] << endl;
    cout << "Unknown: " << num[25] << endl;
  }
 
  d0_info[0] = decay_mode;
  d0_info[1] = r2D0;
  d0_info[2] = deltaD0;
  d0_info[3] = RD0;
  d0_info[4] = FCP;
  d0_info[5] = mnm_gen;
  d0_info[6] = mpn_gen;
  d0_info[7] = double( isKsPiPi );
  d0_info[8] = mpm_gen;
  return d0_info;
}