QCMCFilter Class Reference

#include <QCMCFilter.h>

Inheritance diagram for QCMCFilter:

Public Member Functions
	QCMCFilter (const std::string &name, ISvcLocator *pSvcLocator)
std::vector< double >	findD0Decay (int charm)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()

Detailed Description

Definition at line 17 of file QCMCFilter.h.

Constructor & Destructor Documentation

QCMCFilter()

QCMCFilter::QCMCFilter	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 244 of file QCMCFilter.cxx.

                                                                          :
        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); 
        }

Member Function Documentation

execute()

StatusCode QCMCFilter::execute

(

)

Definition at line 886 of file QCMCFilter.cxx.

                               {
 
    //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
    }
 
}

finalize()

StatusCode QCMCFilter::finalize

(

)

Definition at line 1234 of file QCMCFilter.cxx.

                                {
    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;
}

findD0Decay()

vector< double > QCMCFilter::findD0Decay

(

int

charm

)

Definition at line 1567 of file QCMCFilter.cxx.

{
    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;
}

Referenced by execute().

initialize()

StatusCode QCMCFilter::initialize

(

)

Definition at line 283 of file QCMCFilter.cxx.

                                  {
    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;
}  

---

The documentation for this class was generated from the following files:

• QCMCFilter.h
• QCMCFilter.cxx