QCMCFilter Class Reference

#include <QCMCFilter.h>

List of all members.

Public Member Functions

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

Private Attributes

IDataProviderSvc * m_evtSvc
HepPDT::ParticleDataTable * m_particleTable
int m_nUnknownEvents
int m_nUnknownDecays
int m_nD0D0barEvents
int m_nD0bar
int m_nDpDmEvents
int m_nD0D0barDiscarded
int m_nDpDmDiscarded
int m_nCPPlus
int m_nCPMinus
int m_nFlavoredCFD0
int m_nFlavoredCSD0
int m_nFlavoredDCSD0
int m_nSL
int m_nDalitz
double m_dalitzNumer1
double m_dalitzNumer2
double m_dalitzDenom
bool m_invertSelection
double m_x
double m_y
double m_rCF
double m_zCF
bool m_wCFSign
double m_rCS
double m_zCS
bool m_wCSSign
double m_dplusDminusWeight
int m_dalitzModel
bool m_useNewWeights
double m_rwsCF
double m_rwsCS
double m_deltaCF
double m_deltaCS
double m_largestWeight

Detailed Description

Definition at line 17 of file QCMCFilter.h.

Constructor & Destructor Documentation

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

Definition at line 178 of file QCMCFilter.cxx.

References m_dalitzModel, m_dplusDminusWeight, m_invertSelection, m_rCF, m_rCS, m_useNewWeights, m_wCFSign, m_wCSSign, m_x, m_y, m_zCF, and m_zCS.

00178                                                                       :
00179   Algorithm(name, pSvcLocator){
00180   //Declare the properties
00181   declareProperty("x",                 m_x=0.);
00182   declareProperty("y",                 m_y= 0.0); //For values outside of (-.11 , .06) measured y value may vary from input y
00183   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)
00184   declareProperty("zForCFModes",       m_zCF= 2.0);
00185   declareProperty("wSignForCFModes",   m_wCFSign= true);
00186   declareProperty("rForCSModes",       m_rCS= 1.0);
00187   declareProperty("zForCSModes",       m_zCS= -2.0);
00188   declareProperty("wSignForCSModes",   m_wCSSign= true);
00189   declareProperty("DplusDminusWeight", m_dplusDminusWeight= 0.5);
00190   declareProperty("dalitzModel",       m_dalitzModel= 0); // 0 CLEO-c, 1 Babar, 2 Belle
00191   declareProperty("UseNewWeights",     m_useNewWeights= true); 
00192   declareProperty("InvertSelection",   m_invertSelection= false); 
00193 }

Member Function Documentation

StatusCode QCMCFilter::execute (  ) 

Definition at line 530 of file QCMCFilter.cxx.

References abs, cos(), dalitzDenom_fil, dalitzNumer1_fil, dalitzNumer2_fil, Bes_Common::DEBUG, findD0Decay(), Bes_Common::INFO, kD0BarID, kD0ID, kDalitz, kDpID, kPsi3770ID, m_dplusDminusWeight, m_evtSvc, m_invertSelection, m_keptModeCounter(), m_largestWeight, m_modeCounter(), m_nD0bar, m_nD0D0barDiscarded, m_nD0D0barEvents, m_nDpDmDiscarded, m_nDpDmEvents, m_nUnknownEvents, m_rho_flag, m_weights(), m_x, m_y, EventModel::MC::McParticleCol, msgSvc(), Event::McParticle::particleProperty(), runNo, sin(), weight, and x.

00530                                {
00531   
00532   //interface to event data service
00533   ISvcLocator* svcLocator = Gaudi::svcLocator();
00534   StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
00535   if (sc.isFailure())
00536     std::cout<<"Could not access EventDataSvc!"<<std::endl;
00537 
00538   setFilterPassed(false); 
00539 
00540   MsgStream log(msgSvc(), name());
00541   log<< MSG::INFO << "in execute()" << endmsg;
00542 
00543   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00544 
00545   if(!eventHeader)
00546     {
00547       cout<<" eventHeader "<<endl;
00548       return StatusCode::FAILURE;
00549     } 
00550   long runNo = eventHeader->runNumber();
00551   long event = eventHeader->eventNumber();
00552   log << MSG::DEBUG << "run, evtnum = "
00553       << runNo << " , "
00554       << event << endmsg;
00555 
00556   // Get McParticle collection
00557   SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
00558 
00559   //Check if event is Psi(3770)
00560   int psipp_daughter = 0;
00561   int d0_count = 0;
00562   int chd_count = 0;
00563   for( Event::McParticleCol::iterator it = mcParticles->begin();  it != mcParticles->end(); it++) 
00564     {     
00565       int pdg_code = (*it)->particleProperty();
00566       if ((*it)->primaryParticle()) continue;
00567       Event::McParticle  it2 =  (*it)->mother();
00568       int mother_pdg = 0;
00569       //finds if the particle is daughter of Psi(3770)
00570       mother_pdg = it2.particleProperty();
00571       if (mother_pdg == kPsi3770ID) 
00572         {
00573           psipp_daughter++; //Found daughter
00574           if (abs(pdg_code) == kD0ID) d0_count++; 
00575           if (abs(pdg_code) == kDpID) chd_count++;
00576           if (pdg_code == kD0BarID)   m_nD0bar++;  
00577         }
00578     }
00579 
00580   if( psipp_daughter == 0 ) //didn't find Psi(3770)
00581    {
00582       ++m_nUnknownEvents ;
00583 
00584       if( m_invertSelection )
00585         {
00586           log<< MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg ;
00587           return StatusCode::SUCCESS; //discard event
00588         }
00589       else
00590         {
00591           // -------- Write to root file
00592           setFilterPassed(true); 
00593           
00594           return StatusCode::SUCCESS; //kept event
00595         }
00596    }
00597   
00598    log<< MSG::DEBUG << "Found psi(3770) -->" << endmsg ;
00599 
00600    // Check that top particle has only two daughters + possible FSR photon
00601    if ( psipp_daughter > 3 )
00602      {
00603        ++m_nUnknownEvents ;
00604 
00605       if( m_invertSelection )
00606         {
00607           log<< MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg ;
00608           return StatusCode::SUCCESS; //discard event
00609         }
00610       else
00611         {
00612           // -------- Write to root file
00613           setFilterPassed(true);
00614            
00615           
00616           return StatusCode::SUCCESS; //kept event
00617         }
00618      }
00619 
00620    //Check if D+D- state
00621    if (chd_count == 2) {     // D+D- event
00622      ++m_nDpDmEvents ;
00623      
00624      // D+D- must be rewieghted otherwise too many D+D- events relative to D0D0bar.
00625      double random = RandFlat::shoot() ;
00626      
00627      if( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
00628          ( random > m_dplusDminusWeight &&  m_invertSelection ) )
00629        {
00630          // -------- Write to root file
00631          setFilterPassed(true);
00632           
00633          
00634          return StatusCode::SUCCESS; //event kept
00635        }
00636      else
00637        {
00638          log<< MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg ;
00639          ++m_nDpDmDiscarded ;
00640          return StatusCode::SUCCESS; //event discarded
00641        }
00642    }
00643 
00644    //Check if D0D0Bar
00645    if (d0_count != 2) { 
00646      if( !m_invertSelection )
00647        {
00648          log<< MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg ;
00649          return StatusCode::SUCCESS; //discard event
00650        }
00651      else
00652        {
00653          // -------- Write to root file
00654          setFilterPassed(true);
00655           
00656          
00657          return StatusCode::SUCCESS; //kept event
00658        }
00659    }
00660    
00661    log<< MSG::INFO << "D0-D0bar event." << endmsg ;
00662    m_nD0D0barEvents++ ;
00663 
00664    m_rho_flag = 0;
00665    
00666    //Get D0-D0bar modes and info
00667    vector<double> d0_decay = findD0Decay(1);  
00668    vector<double> d0bar_decay = findD0Decay(-1);
00669 
00670    int d0mode = d0_decay[0];
00671    int d0barmode = d0bar_decay[0];
00672    m_modeCounter[d0mode][d0barmode]++;
00673 
00674 
00675 
00676    //if any event unidentified remove
00677    if (d0_decay[0] == 9 || d0bar_decay[0] == 9 ) { 
00678      if( !m_invertSelection )
00679        {
00680          log<< MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg ;
00681          return StatusCode::SUCCESS; //discard event
00682        }
00683      else
00684        {
00685          // -------- Write to root file
00686          setFilterPassed(true);
00687           
00688          
00689          
00690          return StatusCode::SUCCESS; //kept event
00691        }
00692    }
00693    
00694    // Loop over D's
00695    double r2D0 = d0_decay[1] ;
00696    double deltaD0 = d0_decay[2] ;
00697    double r2D0bar = d0bar_decay[1] ;
00698    double deltaD0bar = d0bar_decay[2] ;
00699    
00700    
00701    //Weight based on DDbar combination
00702    double weight ;
00703    if( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz )
00704      {
00705        double r2prod = r2D0 * r2D0bar ;
00706        double x = m_x ;
00707        double y = m_y ;
00708           
00709        double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
00710        double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
00711        double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
00712        double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
00713           
00714        weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
00715        weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
00716        weight /= m_largestWeight ;
00717 
00718        //get dalitz graph information before cut
00719        int k = -10;
00720        double m2i= 0;
00721        double m2j= 0;
00722        double m2k= 0;
00723        if( d0_decay[0] == kDalitz) {
00724          k = d0bar_decay[0];
00725          m2i= d0_decay[3];
00726          m2j= d0_decay[4];
00727          m2k= d0_decay[6];}
00728        if( d0bar_decay[0] == kDalitz) {
00729          k = d0_decay[0];
00730          m2i= d0bar_decay[3];
00731          m2j= d0bar_decay[4];
00732          m2k= d0bar_decay[6];} 
00733 
00734      }
00735    else
00736      {
00737        weight = m_weights[ d0_decay[0] ][ d0bar_decay[0] ] ;
00738      }
00739    
00740    double random = RandFlat::shoot() ;
00741    log<< MSG::DEBUG << "type:  " << d0_decay[0] << " " << d0bar_decay[0] << ", weight " <<
00742      weight << " <? random " << random << endmsg ;
00743    
00744    if( ( random < weight && !m_invertSelection ) || ( random > weight && m_invertSelection ) )
00745      {
00746        // -------- Write to root file
00747        setFilterPassed(true);
00748        
00749        //get dalitz graph information after filter
00750        int k = -10;
00751        double m2i= 0;
00752        double m2j= 0;
00753        double m2k= 0;
00754        double cosd = cos( deltaD0 ) ;
00755        double sind = sin( deltaD0 ) ;
00756        if( d0_decay[0] == kDalitz) {
00757          k = d0bar_decay[0];
00758          m2i= d0_decay[3];
00759          m2j= d0_decay[4];
00760          m2k= d0_decay[6];
00761          if ( m2j - m2i > 0. ) {  //avoids doublecounting the branching fraction
00762            dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
00763            dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
00764            dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
00765        }
00766        if( d0bar_decay[0] == kDalitz) {
00767          k = d0_decay[0];
00768          m2i= d0bar_decay[3];
00769          m2j= d0bar_decay[4];
00770          m2k= d0bar_decay[6]; 
00771          cosd = cos( deltaD0bar ) ;
00772          sind = sin( deltaD0bar ) ;
00773          if( m2j - m2i < 0. ) { //avoids doublecounting the branching fraction
00774          dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
00775          dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
00776          dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
00777        } 
00778  
00779        m_keptModeCounter[d0mode][d0barmode]++;
00780        return StatusCode::SUCCESS; //event kept
00781      }
00782    else
00783      {
00784        log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg ;
00785        
00786        ++m_nD0D0barDiscarded ;
00787        return StatusCode::SUCCESS; //event discarded
00788      }
00789 }

StatusCode QCMCFilter::finalize (  ) 

Definition at line 793 of file QCMCFilter.cxx.

References dalitzDenom_fil, dalitzNumer1_fil, dalitzNumer2_fil, genRecEmupikp::i, Bes_Common::INFO, ganga-rec::j, m_dalitzDenom, m_dalitzNumer1, m_dalitzNumer2, m_keptModeCounter(), m_modeCounter(), m_nCPMinus, m_nCPPlus, m_nD0bar, m_nD0D0barDiscarded, m_nD0D0barEvents, m_nDalitz, m_nDpDmDiscarded, m_nDpDmEvents, m_nFlavoredCFD0, m_nFlavoredCSD0, m_nFlavoredDCSD0, m_nSL, m_nUnknownDecays, m_nUnknownEvents, m_rCF, m_rCS, m_wCFSign, m_wCSSign, m_weights(), m_x, m_zCF, m_zCS, msgSvc(), nD0, and num.

00793                                 {
00794   MsgStream log(msgSvc(), name());
00795   log << MSG::INFO << "in finalize()" << endmsg;
00796   log << MSG::INFO << "Number of unknown events: "  << m_nUnknownEvents << endmsg ;
00797   log << MSG::INFO << "Number of D+D- events seen: " << m_nDpDmEvents << endmsg ;
00798   log << MSG::INFO << "Number of D+D- events discarded: "  << m_nDpDmDiscarded << endmsg ;
00799   if( m_nDpDmEvents > 0 )
00800     {
00801       log << MSG::INFO <<"Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents ) << endmsg ;
00802     }
00803 
00804   log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar << endmsg ;
00805   log << MSG::INFO << "Number of D0D0bar events discarded: " << m_nD0D0barDiscarded << endmsg ;
00806   if( m_nD0D0barEvents > 0 &&  m_nCPPlus > 0 && m_nCPMinus > 0)
00807     {
00808       log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double( m_nD0D0barEvents ) << endl <<
00809         "Fraction unknown decays: " <<  0.5 * double( m_nUnknownDecays ) / double( m_nD0D0barEvents ) <<endmsg ;
00810       
00811       double nD0 = 2. * m_nD0D0barEvents ;
00812       double brSL = m_nSL / nD0 ;
00813       double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
00814       double brCF = m_nFlavoredCFD0 / nD0 ;
00815       double dBrCF = sqrt( brCF * ( 1. - brCF ) / nD0 ) ;
00816       double brCS = m_nFlavoredCSD0 / nD0 ;
00817       double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
00818       double brDCS = m_nFlavoredDCSD0 / nD0 ;
00819       double dBrDCS = sqrt( brDCS * ( 1. - brDCS ) / nD0 ) ;
00820       double brCPPlus = m_nCPPlus / nD0 ;
00821       double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
00822       double brCPMinus = m_nCPMinus / nD0 ;
00823       double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
00824       double brDalitz = m_nDalitz / nD0 ;
00825       double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
00826       double fdBrDalitz = dBrDalitz / brDalitz ;
00827       double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
00828       double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
00829       double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
00830       double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
00831       double dalitzDenomNorm = m_dalitzDenom / nD0 ;
00832       double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
00833 
00834       //after the filter was applied
00835       double fil_SL = 0, fil_FlavoredCFD0 = 0, fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0 =0, fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0;
00836       for( int i = 0 ; i < 9 ; ++i )
00837         {
00838           fil_SL += m_keptModeCounter[i][4] + m_keptModeCounter[4][i] + m_keptModeCounter[i][5] + m_keptModeCounter[5][i] ;
00839           fil_FlavoredCFD0 += m_keptModeCounter[i][1] + m_keptModeCounter[0][i];
00840           fil_FlavoredCSD0 += m_keptModeCounter[i][2] + m_keptModeCounter[2][i] + m_keptModeCounter[i][3] + m_keptModeCounter[3][i];
00841           fil_FlavoredDCSD0 += m_keptModeCounter[i][0] + m_keptModeCounter[1][i];
00842           fil_CPPlus += m_keptModeCounter[i][6] + m_keptModeCounter[6][i];
00843           fil_CPMinus += m_keptModeCounter[i][7] + m_keptModeCounter[7][i];
00844           fil_Dalitz += m_keptModeCounter[i][8] + m_keptModeCounter[8][i];
00845         }
00846       double nD0_fil = 2. * ( m_nD0D0barEvents -  m_nD0D0barDiscarded ) ;
00847       double brSL_fil = fil_SL / nD0_fil ;
00848       double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
00849       double brCF_fil = fil_FlavoredCFD0 / nD0_fil ;
00850       double dBrCF_fil = sqrt( brCF_fil * ( 1. - brCF_fil ) / nD0_fil ) ;
00851       double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
00852       double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
00853       double brDCS_fil = fil_FlavoredDCSD0 / nD0_fil ;
00854       double dBrDCS_fil = sqrt( brDCS_fil * ( 1. - brDCS_fil ) / nD0_fil ) ;
00855       double brCPPlus_fil = fil_CPPlus / nD0_fil ;
00856       double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
00857       double brCPMinus_fil = fil_CPMinus / nD0_fil ;
00858       double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
00859       double brDalitz_fil = fil_Dalitz / nD0_fil ;
00860       double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
00861       double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
00862       double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
00863       double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
00864       double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
00865       double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
00866       double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
00867       double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
00868 
00869 
00870       log << MSG::INFO <<
00871         "Original number of SL decays: " << m_nSL <<"  ==> Br(SL) = " << brSL << " +- " << dBrSL << endl <<
00872         "Filtered number of SL decays: " << fil_SL <<"  ==> Br(SL) = " << brSL_fil << " +- " << dBrSL_fil << endl <<
00873         "Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<m_nFlavoredCFD0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0 << endl <<
00874         "  ==> Br(CF) = " << brCF << " +- " << dBrCF << endl <<
00875         "  ==> Br(CS) = " << brCS << " +- " << dBrCS << endl <<
00876         "  ==> Br(DCS) = " << brDCS << " +- " << dBrDCS << endl <<
00877         "Filtered number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0 << "/" << fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0 << endl <<
00878         "  ==> Br(CF) = " << brCF_fil << " +- " << dBrCF_fil << endl <<
00879         "  ==> Br(CS) = " << brCS_fil << " +- " << dBrCS_fil << endl <<
00880         "  ==> Br(DCS) = " << brDCS_fil << " +- " << dBrDCS_fil << endl <<
00881         
00882         "Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl << 
00883         "  ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl << 
00884         "  ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl << 
00885         "Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl << 
00886         "  ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl << 
00887         "  ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl << 
00888 
00889         "Original number of Dalitz: " << m_nDalitz << endl << 
00890         "  ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
00891         "  y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
00892         "  numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
00893         "  denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
00894         "Filtered number of Dalitz: " << fil_Dalitz << endl << 
00895         "  ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
00896         "  y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
00897         "  numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
00898         "  denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg ;
00899 
00900       cout<<""<<endl;
00901 
00902       HepMatrix    differencetoWeight(10,10,0);
00903       for( int i = 0 ; i < 9 ; ++i )  { 
00904         for( int j = 0 ; j < 9 ; ++j )  {
00905           if (m_modeCounter[i][j] != 0) {
00906             differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] - m_weights[i][j]; } } }
00907 
00908       log << MSG::INFO << "Weight matrix" << m_weights <<  endmsg ;
00909       log << MSG::INFO << "D0 Modes before filter" << m_modeCounter <<  endmsg ;
00910       log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter <<  endmsg ;
00911       log << MSG::INFO << "Compare percentage difference to weights " << differencetoWeight <<  endmsg ;
00912 
00913       //Calculating and comparing Y before/after the filter
00914       // To calculate y, we want half the CS BF
00915       brCS /= 2. ;
00916       dBrCS /= 2. ;
00917       brCS_fil /= 2 ;
00918       dBrCS_fil /= 2 ;
00919 
00920       double y, dy, y_fil, dy_fil ;
00921       double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
00922       
00923       //Original MC
00924       double rCF = m_rCF ;
00925       double zCF = m_zCF ;
00926       double rCF2 = rCF * rCF ;
00927       double zCF2 = zCF * zCF ;
00928       double rCS = m_rCS ;
00929       double zCS = m_zCS ;
00930       double rCS2 = rCS * rCS ;
00931       double zCS2 = zCS * zCS ;
00932       
00933       double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
00934       double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
00935       
00936       double numFactCF =
00937         -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ; 
00938       double numFactCS =
00939         -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ; 
00940       double denFactCF = 0.5 * rCF2 * zCF2 ;  
00941       double denFactCS = 0.5 * rCS2 * zCS2 ;  
00942       double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ; 
00943       
00944       double num =
00945         brCPPlus - brCPMinus + brCF * numFactCF + brCS * numFactCS +
00946         dalitzNumerNorm ; 
00947       double den =
00948         1. - brCPPlus - brCPMinus - brCF * denFactCF - brCS * denFactCS +
00949         dalitzDenomNorm ; 
00950       y = num / den ;
00951       dy = sqrt(
00952                 ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
00953                 ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
00954                 ( numFactCF * den + denFactCF * num ) * dBrCF *
00955                 ( numFactCF * den + denFactCF * num ) * dBrCF +
00956                 ( numFactCS * den + denFactCS * num ) * dBrCS *
00957                 ( numFactCS * den + denFactCS * num ) * dBrCS +
00958                 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz *
00959                 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz
00960                 ) / ( den * den ) ;
00961       
00962       double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] + m_modeCounter[4][6] +  m_modeCounter[5][6] ;
00963       double n_cpminussl = m_modeCounter[7][4] + m_modeCounter[7][5] + m_modeCounter[4][7] +  m_modeCounter[5][7] ;     
00964       double a_cpsl = ( n_cpplussl  * m_nCPMinus ) / ( n_cpminussl * m_nCPPlus  ) ;
00965       double b_cpsl = ( n_cpminussl  * m_nCPPlus ) / ( n_cpplussl * m_nCPMinus  ) ;
00966       y_cpsl = 0.25 * ( a_cpsl  - b_cpsl ) ;
00967       dy_cpsl = 0.25 * ( a_cpsl + b_cpsl )  * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) ) ; 
00968       
00969       
00970       //Filtered MC
00971       double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
00972       double num_fil =
00973         brCPPlus_fil - brCPMinus_fil + brCF_fil * numFactCF + brCS_fil * numFactCS +
00974         dalitzNumerNorm_fil ;
00975       double den_fil =
00976         1. - brCPPlus_fil - brCPMinus_fil - brCF_fil * denFactCF - brCS_fil * denFactCS +
00977         dalitzDenomNorm_fil ;  
00978       y_fil = num_fil / den_fil ;
00979       dy_fil = sqrt(
00980                     ( num_fil + den_fil ) * ( num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil +
00981                     ( num_fil - den_fil ) * ( num_fil - den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil +
00982                     ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil *
00983                     ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil +
00984                     ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil *
00985                     ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
00986                     ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil *
00987                     ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil
00988                     ) / ( den_fil * den_fil ) ;
00989       
00990       double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] + m_keptModeCounter[4][6] +  m_keptModeCounter[5][6] ;
00991       double fil_cpminussl = m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] +  m_keptModeCounter[5][7] ;   
00992       a_cpsl = ( fil_cpplussl  * fil_CPMinus ) / ( fil_cpminussl * fil_CPPlus  ) ;
00993       b_cpsl = ( fil_cpminussl  * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus  ) ;
00994       y_fil_cpsl = 0.25 * ( a_cpsl  - b_cpsl ) ;
00995       dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ; 
00996       
00997       
00998       
00999       log << MSG::INFO <<"BR Method : Parent MC    ==> y = " << y << " +- " << dy << endmsg ;
01000       log << MSG::INFO <<"BR Method : Filtered MC  ==> y = " << y_fil << " +- " << dy_fil << endmsg ;
01001       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 corrilated"<< endmsg ;
01002       log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " << dy_fil_cpsl << endmsg ;
01003       
01004     }
01005 
01006   return StatusCode::SUCCESS;
01007 }

vector< double > QCMCFilter::findD0Decay ( int  charm  ) 

Definition at line 1012 of file QCMCFilter.cxx.

References arg(), conj(), cos(), Bes_Common::DEBUG, false, genRecEmupikp::i, Bes_Common::INFO, kA00ID, kA0MinusID, kA0PlusID, kA1MinusID, kA1PlusID, kCPMinus, kCPPlus, kD0BarID, kD0ID, kDalitz, kEMinusID, kEPlusID, kEtaID, kEtaPrimeID, kF0ID, kF0m1500ID, kF2ID, kFlavoredCF, kFlavoredCFBar, kFlavoredCS, kFlavoredCSBar, kFPrime0ID, kGammaFSRID, kGammaID, kK0BarID, kK0ID, kK0LID, kK0SID, kK0Star0BarID, kK0Star0ID, kK0StarMinusID, kK0StarPlusID, kK10BarID, kK10ID, kK1MinusID, kK1PlusID, kK1Prime0BarID, kK1Prime0ID, kK1PrimeMinusID, kK1PrimePlusID, kK2StarMinusID, kK2StarPlusID, kKMinusID, kKPlusID, kKStar0BarID, kKStar0ID, kKStarMinusID, kKStarPlusID, kMuMinusID, kMuPlusID, kNuEBarID, kNuEID, kNuMuBarID, kNuMuID, kOmegaID, kPhiID, kPi0ID, kPiMinusID, kPiPlusID, kRho0ID, kRho2S0ID, kRho2SMinusID, kRho2SPlusID, kRhoMinusID, kRhoPlusID, kSLMinus, kSLPlus, m_dalitzDenom, m_dalitzModel, m_dalitzNumer1, m_dalitzNumer2, m_deltaCF, m_deltaCS, m_nCPMinus, m_nCPPlus, m_nDalitz, m_nFlavoredCFD0, m_nFlavoredCSD0, m_nFlavoredDCSD0, m_nSL, m_rCF, m_rCS, m_rwsCF, m_rwsCS, m_x, EventModel::MC::McParticleCol, Event::McParticle::mother(), msgSvc(), num, Event::McParticle::particleProperty(), PI, sin(), t(), true, xmk0, and xmpion.

Referenced by execute().

01013 {
01014   MsgStream log(msgSvc(), name());
01015   log << MSG::INFO <<"  In findD0Decay " << endmsg ;
01016   
01017   vector<double> d0_info(7);
01018   for(int i=0; i< 7; i++)  d0_info[i]=0;
01019   double decay_mode = 9; //d0_info[0] = decay mode;
01020   double r2D0 = -1;       //d0_info[1] = r2D0;
01021   double deltaD0 = -1;    //d0_info[2] = deltaD0;
01022   
01023   double num[26];
01024   for(int i=0; i< 26; i++)  num[i]=0;
01025   // Number of each partile type
01026   int kPiPlus        =  0;// num[0] = Pi+, a0+ 
01027   int kPiMinus       =  1;// num[1] = Pi-, a0-  
01028   int kPi0           =  2;// num[2] = Pi0  
01029   int kEta           =  3;// num[3] = Eta  
01030   int kEtaPrime      =  4;// num[4] = Eta Prime  
01031   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)
01032   int kUFVPlus       =  6;// num[6] = unflavored positive (axial) vectors: rho+, a1+, rho(2S)+  
01033   int kUFVMinus      =  7;// num[7] = unflavored negative (axial) vectors: rho-, a1-, rho(2S)-  
01034   int kRho0          =  8;// num[8] = Rho0, Rho(2S)0   
01035   int kOmega         =  9;// num[9] = Omega  
01036   int kPhi           = 10;// num[10] = Phi
01037   int kKPlus         = 11;// num[11] = K+  
01038   int kKMinus        = 12;// num[12] = K-  
01039   int kK0S           = 13;// num[13] = K short  
01040   int kK0L           = 14;// num[14] = K long  
01041   int kFVPlus        = 15;// num[15] = flavored positive (axial) vectors: K*+, K_1+, K_2*+
01042   int kFVMinus       = 16;// num[16] = flavored negative (axial) vectors: K*-, K_1-, K_2*- 
01043   int kKStar0        = 17;// num[17] = K*0 
01044   int kKStar0Bar     = 18;// num[18] = anti-K*0 
01045   int kK10           = 19;// num[19] = K1_0
01046   int kK10Bar        = 20;// num[20] = anti-K1_0
01047   int kLeptonPlus    = 21;// num[21] =  LeptonPlus: e+, mu+ 
01048   int kLeptonMinus   = 22;// num[22] =  LeptonMinus: e-, mu- 
01049   int kNeutrino      = 23;// num[23] =  Neutrino: nu, nubar 
01050   int kGamma         = 24;// num[24] =  gamma 
01051   int kUnknown       = 25;// num[25] =  Not identified particle   
01052   //int kGammaFSR      GammaFSR is ignored as we are interested in the PreFSR quantum # 
01053 
01054   // Get McParticle collection
01055   SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
01056   int d0_pdg = charm == 1 ? kD0ID : kD0BarID;
01057 
01058   HepLorentzVector piPlus, piMinus, k0;
01059 
01060   for( Event::McParticleCol::iterator it = mcParticles->begin();  it != mcParticles->end(); it++) 
01061     {     
01062       int pdg_code = (*it)->particleProperty();
01063       if ((*it)->primaryParticle()) continue;
01064       Event::McParticle  it2 =  (*it)->mother();
01065       int mother_pdg = 0;
01066 
01067       //finds if the particle is related to the d0 we need
01068       mother_pdg = it2.particleProperty();
01069 
01070       //special cases
01071       if (pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID || pdg_code == kK0BarID 
01072           || pdg_code == kK10ID || pdg_code == kK10BarID 
01073           || pdg_code == kK0Star0ID || pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID ) continue;  //don't count special intermediate states
01074 
01075       if (mother_pdg == kK0ID || mother_pdg == kK0BarID) {//Looks another step up if mother was k0 or k0bar 
01076         it2 = it2.mother();
01077         mother_pdg = it2.particleProperty(); }
01078 
01079       if (mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID) {//Looks another step up if mother was k*0 or k*0bar
01080         //if code is found to be K*0 -> K+Pi- we save  it as a KStar0 and K*0Bar -> K-Pi+ as KStar0Bar
01081         //if is it a neutral decay K*0(Bar)->K0(bar) Pi0 or Gamma.  We save them as K0(bar) and Pi0 or Gamma.
01082         if (pdg_code == kPiPlusID || pdg_code == kPiMinusID) continue;
01083         if (mother_pdg == kKStar0ID && pdg_code == kKPlusID)     pdg_code = kKStar0ID;
01084         if (mother_pdg == kKStar0BarID && pdg_code == kKMinusID) pdg_code = kKStar0BarID;
01085         it2 = it2.mother();
01086         mother_pdg = it2.particleProperty();
01087       }
01088         
01089       if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID ) { // Looks another step up if mother was K_0
01090         it2 = it2.mother();
01091         mother_pdg = it2.particleProperty();
01092       }
01093 
01094       if (mother_pdg == kK10ID || mother_pdg == kK10BarID) {//Looks another step up if mother was k_10 or k_10bar
01095         it2 = it2.mother();
01096         mother_pdg = it2.particleProperty();
01097       }
01098       
01099       
01100       
01101       //Identifies what particles are the daughters
01102       if (mother_pdg == d0_pdg)  //Found daughter 
01103         {
01104           if (pdg_code ==  kPiPlusID || pdg_code ==  kA0PlusID )                   num[0]++;
01105           else if (pdg_code ==  kPiMinusID || pdg_code ==  kA0MinusID )            num[1]++;
01106           else if (pdg_code ==  kPi0ID )                                           num[2]++;
01107           else if (pdg_code ==  kEtaID )                                           num[3]++;
01108           else if (pdg_code ==  kEtaPrimeID )                                      num[4]++;
01109           else if (pdg_code ==  kF0ID || pdg_code ==  kA00ID 
01110                    || pdg_code ==  kFPrime0ID|| pdg_code ==  kF0m1500ID
01111                    || pdg_code ==  kF2ID)                                          num[5]++;
01112           else if (pdg_code ==  kRhoPlusID || pdg_code ==  kRho2SPlusID
01113                    || pdg_code ==  kA1PlusID )                                     num[6]++;
01114           else if (pdg_code ==  kRhoMinusID || pdg_code ==  kRho2SMinusID
01115                    || pdg_code ==  kA1MinusID)                                     num[7]++;
01116           else if (pdg_code ==  kRho0ID || pdg_code ==  kRho2S0ID)                 num[8]++;
01117           else if (pdg_code ==  kOmegaID )                                         num[9]++;
01118           else if (pdg_code ==  kPhiID )                                           num[10]++;
01119           else if (pdg_code ==  kKPlusID )                                         num[11]++;
01120           else if (pdg_code ==  kKMinusID )                                        num[12]++;
01121           else if (pdg_code ==  kK0SID )                                           num[13]++;
01122           else if (pdg_code ==  kK0LID )                                           num[14]++;
01123           else if (pdg_code ==  kKStarPlusID || pdg_code == kK1PlusID 
01124                    || pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID
01125                    || pdg_code == kK0StarPlusID)                                   num[15]++;
01126           else if (pdg_code ==  kKStarMinusID || pdg_code == kK1MinusID 
01127                    || pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID
01128                    || pdg_code == kK0StarMinusID )                                 num[16]++;
01129           else if (pdg_code ==  kKStar0ID )                                        num[17]++;
01130           else if (pdg_code ==  kKStar0BarID )                                     num[18]++;
01131           else if (pdg_code ==  kK10ID || pdg_code ==  kK1Prime0ID)                num[19]++;
01132           else if (pdg_code ==  kK10BarID || pdg_code ==  kK1Prime0BarID)          num[20]++;
01133           else if (pdg_code ==  kEPlusID || pdg_code ==  kMuPlusID )               num[21]++;
01134           else if (pdg_code ==  kEMinusID || pdg_code ==  kMuMinusID )             num[22]++;
01135           else if (pdg_code ==  kNuEID || pdg_code ==  kNuEBarID
01136                    || pdg_code ==  kNuMuID || pdg_code ==  kNuMuBarID )            num[23]++;
01137           else if (pdg_code ==  kGammaID )                                         num[24]++;
01138           else if (pdg_code ==  kGammaFSRID )                                       continue;
01139           else {
01140             num[25]++;
01141             cout <<"Unknown particle:  "<< pdg_code << endl;
01142           }
01143 
01144           //if (pdg_code ==  kA10ID )            num[0]++; // not present
01145 
01146           //Enter Momentums for Dalitz canidates
01147           //It would be prefered if we could get PreFSR Momentum
01148           if (pdg_code ==  kPiPlusID)  piPlus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
01149           if (pdg_code ==  kPiMinusID) piMinus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
01150           if (pdg_code ==  kK0SID)     k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
01151           if (pdg_code ==  kK0LID)     k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);      
01152         }
01153     }
01154   
01155   // D decay daughters have been tabulated.  Now, classify decay.
01156   int nDaughters = 0 ;
01157   for( int i = 0 ; i < 26 ; i++ )
01158     {
01159       nDaughters += num[ i ] ;
01160     }
01161   
01162   int nUFP0 = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] ;
01163   int nUFV0 = num[ kRho0 ] + num[ kPhi ] + num[ kOmega ] + num[ kGamma ] ;
01164   int nFV0 = num[ kKStar0 ] + num[ kK10 ] ;
01165   int nFV0Bar = num[ kKStar0Bar ] + num[ kK10Bar ] ;
01166   int nStrange = num[ kKMinus ] + num[ kFVMinus ] + nFV0Bar ;
01167   int nAntiStrange = num[ kKPlus ] + num[ kFVPlus ] + nFV0 ;
01168   int nCPPlusEig = num[ kNeutralScalar ] + num[ kRho0 ] + num[ kOmega ] + num[ kPhi ] + num[ kK0S ] + num[ kGamma ] ;
01169   int nCPMinusEig = num[ kPi0 ] + num[ kEta ] +  num[ kEtaPrime ] + num[ kK0L ] ;
01170   int nCPEig = nCPPlusEig + nCPMinusEig ;
01171   int nChargedPiK = num[ kPiPlus ] + num[ kPiMinus ] + num[ kKPlus ] + num[ kKMinus ] ;
01172   int nK0 = num[ kK0S ] + num[ kK0L ] ;
01173   
01174  
01175   //check Dalitz modes: KsPi+Pi- or KlPi+Pi-
01176   double mnm_gen = 0. ;
01177   double mpn_gen = 0. ;
01178   double mpm_gen = 0. ;
01179   bool isKsPiPi = false ;
01180   
01181   if( nK0 == 1 && num[ kPiPlus ] == 1 && num[ kPiMinus ] && nDaughters == 3 )
01182     {
01183       decay_mode = kDalitz ;
01184       k0.boost(-0.011, 0, 0);
01185       piMinus.boost(-0.011, 0, 0);
01186       piPlus.boost(-0.011, 0, 0);
01187       mnm_gen = (k0 + piMinus).m2();
01188       mpn_gen = (k0 + piPlus).m2();
01189       mpm_gen = (piPlus + piMinus).m2();
01190       if( num[ kK0S ] == 1)  isKsPiPi = true ;
01191     }
01192     
01193   // The order of if-statements below matters.
01194   if( decay_mode == kDalitz )
01195     {
01196       ++m_nDalitz ;
01197       
01198       Dalitz t(8) ;
01199       TComplex D0, D0bar;
01200       if (m_dalitzModel == 0) { //Cleo model
01201         D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01202         D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01203       if (m_dalitzModel == 1) { //Babar model
01204         D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01205         D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01206       if (m_dalitzModel == 2) { //Belle model
01207         D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01208         D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01209       
01210       if( charm == 1){
01211         r2D0 = norm( D0bar ) / norm( D0 ) ;
01212         deltaD0 = arg( D0 * conj( D0bar ) ) + ( isKsPiPi ? PI : 0. );
01213         
01214         double cosd = cos( deltaD0 ) ;
01215         double sind = sin( deltaD0 ) ;
01216         if( mpn_gen - mnm_gen > 0. )
01217           {
01218             m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
01219             m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
01220             m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
01221           }
01222       }
01223       else {
01224         r2D0 = norm( D0 ) / norm( D0bar ) ;
01225         deltaD0 =  arg( D0bar * conj( D0 ) ) + ( isKsPiPi ? PI : 0. ) ;
01226 
01227         double cosd = cos( deltaD0 ) ;
01228         double sind = sin( deltaD0 ) ;
01229         if( mpn_gen - mnm_gen < 0. )
01230           {
01231             m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
01232             m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
01233             m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
01234           }
01235       }
01236     }
01237   else if( num[ kLeptonPlus ] == 1 )
01238     {
01239       decay_mode = kSLPlus ;
01240       ++m_nSL ;
01241       
01242       r2D0 = 0. ;
01243        deltaD0 = 0. ;
01244     }
01245   else if( num[ kLeptonMinus ] == 1 )
01246     {
01247       decay_mode = kSLMinus ;
01248       ++m_nSL ;
01249       
01250       r2D0 = 0. ;
01251       deltaD0 = 0. ;
01252     }
01253   //Check CP even modes
01254   else if(
01255           ( nDaughters == 2 &&
01256             ( ( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 ) ||
01257               nUFP0 == 2 ||
01258               num[ kNeutralScalar ] == 2 ||
01259               ( nUFP0 == 1 && nUFV0 == 1 ) ||
01260               ( num[ kNeutralScalar ] == 1 && nUFV0 == 1 ) ||
01261               ( num[ kKPlus ] == 1 && num[ kKMinus ] == 1 ) ||
01262               num[ kK0S ] == 2 || 
01263               num[ kK0L ] == 2 ||
01264               ( num[ kK0L ] == 1 && nUFP0 == 1 ) ||
01265               ( num[ kK0S ] == 1 && num[ kNeutralScalar ] == 1 ) ||
01266               ( num[ kK0L ] == 1 && nUFV0 == 1 ) ) ) ||
01267           ( nDaughters == 3 &&
01268             ( ( nCPPlusEig == 1 && num[ kPi0 ] == 2 ) ||
01269               ( num[ kK0S ] == 3 ) ) ) ||
01270           ( nDaughters == 4 &&
01271             num[ kK0L ] == 1 && num[ kPi0 ] == 3 )
01272           )
01273     {
01274       decay_mode = kCPPlus ;
01275       ++m_nCPPlus ;
01276       
01277       r2D0 = 1. ;
01278       deltaD0 = PI;
01279     }
01280   //Check CP odd modes
01281   else if(
01282           ( nDaughters == 2 && 
01283             ( ( num[ kK0S ] == 1 && nUFP0 == 1 ) ||
01284               ( num[ kK0L ] == 1 && num[ kNeutralScalar ] == 1 ) ||
01285               ( num[ kK0S ] == 1 && nUFV0 == 1 ) ||
01286               ( nUFP0 == 1 && num[ kNeutralScalar ] == 1 ) ) ) ||
01287           ( nDaughters == 3 &&
01288             ( ( nCPMinusEig == 3 && num[ kPi0 ] == 2 ) || // pi0 is subset of CP-
01289               ( num[ kPi0 ] == 3 ) ||
01290               ( num[ kK0L ] == 3 ) ) ) ||
01291           ( nDaughters == 4 &&
01292             num[ kK0S ] == 1 && num[ kPi0 ] == 3 )
01293           )
01294     {
01295       decay_mode = kCPMinus ;
01296       ++m_nCPMinus ;
01297       
01298       r2D0 = 1. ;
01299       deltaD0 = 0. ;
01300     }
01301   else if( nStrange == nAntiStrange + 1 )
01302     {
01303       decay_mode = kFlavoredCF ;
01304       
01305       if( charm == 1 )
01306         {
01307           ++m_nFlavoredCFD0 ;
01308           r2D0 = pow(m_rCF,2) ;
01309           deltaD0 = m_deltaCF ;
01310         }
01311       else
01312         {
01313           ++m_nFlavoredDCSD0 ;
01314           r2D0 = 1. / pow( m_rCF,2 ) ;
01315           deltaD0 = -m_deltaCF ;
01316         }
01317     }
01318   else if( nAntiStrange == nStrange + 1 )
01319     {
01320       decay_mode = kFlavoredCFBar ;
01321       
01322       if( charm == -1 )
01323         {
01324           ++m_nFlavoredCFD0 ;
01325           r2D0 = pow( m_rCF, 2) ;
01326           deltaD0 = m_deltaCF ;
01327         }
01328       else
01329         {
01330           ++m_nFlavoredDCSD0 ;
01331           r2D0 = 1. / pow( m_rCF, 2) ;
01332           deltaD0 = -m_deltaCF ;
01333         }
01334     }
01335   else if( nStrange == nAntiStrange ) // Unflavored or Cabibbo-suppressed
01336     {
01337       if( ( num[ kKPlus ] > 0 &&
01338             ( num[ kKPlus ] == num[ kFVMinus ] ||
01339               num[ kKPlus ] == nFV0Bar ) ) || // for anti-K*0 K+ pi-
01340           ( nK0 > 0 &&
01341             nFV0 != nFV0Bar &&
01342             nK0 == nFV0Bar ) ||
01343           ( num[ kPiPlus ] > 0 &&
01344             num[ kPiPlus ] == num[ kUFVMinus ] ) )
01345         {
01346           decay_mode = kFlavoredCS ;
01347           ++m_nFlavoredCSD0 ;
01348           
01349           if( charm == 1 )
01350             {
01351               r2D0 = pow( m_rCS, 2 ) ;
01352               deltaD0 = m_deltaCS ;
01353             }
01354           else
01355             {
01356               r2D0 = 1. / pow( m_rCS, 2 ) ;
01357               deltaD0 = -m_deltaCS ;
01358             }
01359         }
01360       else if( ( num[ kKMinus ] > 0 && ( num[ kKMinus ] == num[ kFVPlus ] || num[ kKMinus ] == nFV0 ) ) || // for K*0 K- pi+
01361                ( nK0 > 0 && nFV0 != nFV0Bar &&  nK0 == nFV0 ) ||
01362                ( num[ kPiMinus ] > 0 && num[ kPiMinus ] == num[ kUFVPlus ] ) )
01363         {
01364           decay_mode = kFlavoredCSBar ;
01365           ++m_nFlavoredCSD0 ;
01366           
01367           if( charm == -1 )
01368             {
01369               r2D0 = pow( m_rCS, 2 ) ;
01370               deltaD0 = m_deltaCS ;
01371             }
01372           else
01373             {
01374               r2D0 = 1. / pow( m_rCS, 2 ) ;
01375               deltaD0 = -m_deltaCS ;
01376             }
01377         }
01378 
01379       // Self-conjugate mixed-CP final states -- pick CP or non-CP at
01380       // random.  If CP, pick + or - at random.  If non-CP, pick
01381       // flavored or flavoredBar according to the appropriate value of r.
01382       else if( ( nDaughters >= 3 && num[ kPiPlus ] == num[ kPiMinus ] &&
01383                  num[ kKPlus ] == num[ kKMinus ] && nChargedPiK + nCPEig == nDaughters ) ||
01384                nUFV0 == nDaughters ||
01385                ( num[ kKStar0 ] > 0 && num[ kKStar0 ] == num[ kKStar0Bar ] ) ||
01386                ( num[ kK10 ] > 0 && num[ kK10 ] == num[ kK10Bar ] ) ||
01387                ( num[ kUFVPlus ] == 1 && num[ kUFVMinus ] == 1 )
01388                // for rho+rho-; no a1+a1- and rhoa1 final states in DECAY.DEC
01389                )
01390         {
01391           log << MSG::DEBUG <<"    [ Self-conjugate mixed-CP ]" << endmsg ;
01392           
01393           if( RandFlat::shoot() > 0.5 )
01394             {
01395               if( RandFlat::shoot() > 0.5 )
01396                 {
01397                   decay_mode = kCPPlus ;
01398                   ++m_nCPPlus ;
01399                   
01400                   r2D0 = 1. ;
01401                   deltaD0 =  PI;
01402                 }
01403               else
01404                 {
01405                   decay_mode = kCPMinus ;
01406                   ++m_nCPMinus ;
01407                   
01408                   r2D0 = 1. ;
01409                   deltaD0 = 0. ;
01410                 }
01411             }
01412           else
01413             {
01414               // Odd # of K0S or K0L = CF; even # of K0S or K0L = SC.
01415               if( nK0 % 2 )
01416                 {
01417                   // Nov 2007: correct for r2 -> RWS and BR != A^2
01418                   double cutoff = m_rwsCF / ( 1. + m_rwsCF ) ;
01419                   
01420                   if( charm == -1 )
01421                     {
01422                       cutoff = 1. - cutoff ;
01423                     }
01424                   
01425                   log << MSG::DEBUG <<"    [ cutoff = " << cutoff << " ]" << endmsg ;
01426                   
01427                   decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF : kFlavoredCFBar ;
01428                   
01429                   if( ( decay_mode == kFlavoredCF && charm == 1 ) ||
01430                       ( decay_mode == kFlavoredCFBar && charm == -1 ) )
01431                     {
01432                       ++m_nFlavoredCFD0 ;
01433                       
01434                       r2D0 = sqrt( m_rCF ) ;
01435                       deltaD0 = m_deltaCF ;
01436                     }
01437                   else
01438                     {
01439                       ++m_nFlavoredDCSD0 ;
01440                       
01441                       r2D0 = 1. / sqrt( m_rCF ) ;
01442                       deltaD0 = -m_deltaCF ;
01443                     }
01444                 }
01445               else
01446                 {
01447                   // Nov 2007: correct for r2 -> RWS and BR != A^2
01448                   double cutoff = m_rwsCS / ( 1. + m_rwsCS ) ;
01449                   
01450                   if( charm == -1 )
01451                     {
01452                       cutoff = 1. - cutoff ;
01453                     }
01454                   
01455                   log << MSG::DEBUG <<"    [ cutoff = " << cutoff << " ]" << endmsg ;
01456                   
01457                   decay_mode = ( RandFlat::shoot() > cutoff ) ?
01458                     kFlavoredCS : kFlavoredCSBar ;
01459                   ++m_nFlavoredCSD0 ;
01460                   
01461                   if( ( decay_mode == kFlavoredCS && charm == 1 ) ||
01462                       ( decay_mode == kFlavoredCSBar && charm == -1 ) )
01463                     {
01464                       r2D0 = sqrt( m_rCS ) ;
01465                       deltaD0 = m_deltaCS ;
01466                     }
01467                   else
01468                     {
01469                       r2D0 = 1. / sqrt( m_rCS ) ;
01470                       deltaD0 = -m_deltaCS ;
01471                     }
01472                 }
01473             }
01474         }
01475     }
01476 
01477   if (num[kUnknown] >= 1)
01478     {
01479       cout << "decay mode " << decay_mode << endl ;
01480       cout << "D #" << charm << endl ;
01481       cout << "pi+:     " <<  num[ kPiPlus ] << endl ;
01482       cout << "pi-:     " <<  num[ kPiMinus ] << endl ;
01483       cout << "pi0:     " <<  num[ kPi0 ] << endl ;
01484       cout << "eta:     " <<  num[ kEta ] << endl ;
01485       cout << "eta':    " <<  num[ kEtaPrime ] << endl ;
01486       cout << "f0/a0:   " <<  num[ kNeutralScalar ] << endl ;
01487       cout << "UFV+:    " <<  num[ kUFVPlus ] << endl ;
01488       cout << "UFV-:    " <<  num[ kUFVMinus ] << endl ;
01489       cout << "rho0:    " <<  num[ kRho0 ] << endl ;
01490       cout << "omega:   " <<  num[ kOmega ] << endl ;
01491       cout << "phi:     " <<  num[ kPhi ] << endl ;
01492       cout << "K+:      " <<  num[ kKPlus ] << endl ;
01493       cout << "K-:      " <<  num[ kKMinus ] << endl ;
01494       cout << "K0S:     " <<  num[ kK0S ] << endl ;
01495       cout << "K0L:     " <<  num[ kK0L ] << endl ;
01496       cout << "FV+:     " <<  num[ kFVPlus ] << endl ;
01497       cout << "FV-:     " <<  num[ kFVMinus ] << endl ;
01498       cout << "K*0:     " <<  num[ kKStar0 ] << endl ;
01499       cout << "K*0b:    " <<  num[ kKStar0Bar ] << endl ;
01500       cout << "K10:     " <<  num[ kK10 ] << endl ;
01501       cout << "K10b:    " <<  num[ kK10Bar ] << endl ;
01502       cout << "l+:      " <<  num[ kLeptonPlus ] << endl ;
01503       cout << "l-:      " <<  num[ kLeptonMinus ] << endl ;
01504       cout << "nu:      " <<  num[ kNeutrino ] << endl ;
01505       cout << "gamma:   " <<  num[ kGamma ] << endl ;
01506       cout << "Unknown: " <<  num[ 25 ] << endl ;  
01507     }
01508 
01509 
01510   d0_info[0]=decay_mode;
01511   d0_info[1]=r2D0;
01512   d0_info[2]=deltaD0;
01513   d0_info[3]=mnm_gen;
01514   d0_info[4]=mpn_gen;
01515   d0_info[5]= double (isKsPiPi);
01516   d0_info[6]=mpm_gen;
01517   return d0_info;
01518 }

StatusCode QCMCFilter::initialize (  ) 

Definition at line 197 of file QCMCFilter.cxx.

References bin, conj(), calibUtil::ERROR, Bes_Common::FATAL, genRecEmupikp::i, imag(), Bes_Common::INFO, ganga-rec::j, kCPMinus, kCPPlus, kDalitz, kFlavoredCF, kFlavoredCFBar, kFlavoredCS, kFlavoredCSBar, kNDecayTypes, kSLMinus, kSLPlus, m_dalitzModel, m_deltaCF, m_deltaCS, m_largestWeight, m_particleTable, m_rCF, m_rCS, m_rwsCF, m_rwsCS, m_useNewWeights, m_wCFSign, m_wCSSign, m_weights(), m_x, m_y, m_zCF, m_zCS, msgSvc(), num, s, t(), and x.

00197                                   {
00198   MsgStream log(msgSvc(), name());
00199 
00200   log << MSG::INFO << "in initialize()" << endmsg;
00201   
00202   // Get the Particle Properties Service
00203   IPartPropSvc* p_PartPropSvc;
00204   static const bool CREATEIFNOTTHERE(true);
00205   StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00206   if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) 
00207     {
00208       log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
00209       return PartPropStatus;
00210     }
00211   m_particleTable = p_PartPropSvc->PDT();
00212 
00213 
00214   //Calculates parameters here based on Parameter Input by User
00215   // (e.g. run expensive algorithms that are based on parameters
00216   //  specified by user at run-time)
00217 
00218   m_y = (m_y - 0.008)/0.292; 
00219   //corrects the y input so that we get the desired measurable y from CP:Semi-leptonic doubletag method
00220   // This equation is a result of solving equations for y dependence.
00221   //If parent MC Branching ratio's this equation will need resolved. January 2013 -DA
00222   double x = m_x ;
00223   double y = m_y ;
00224   double rCF = m_rCF ;
00225   double zCF = m_zCF ;
00226   double rCF2 = rCF * rCF ;
00227   double zCF2 = zCF * zCF ;
00228   double rCS = m_rCS ;
00229   double zCS = m_zCS ;
00230   double rCS2 = rCS * rCS ;
00231   double zCS2 = zCS * zCS ;
00232   
00233   double rmix = ( x * x + y * y ) / 2. ;
00234   double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
00235   double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
00236   double vCFCSPlus = ( zCF * zCS + wCF * wCS ) / 2. ;
00237   double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2. ;
00238   
00239   m_deltaCF = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. ) ;
00240   m_deltaCS = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. ) ;
00241 
00242   double bCF = 1. + 0.5 * rCF * ( zCF * y + wCF * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
00243   double bCS = 1. + 0.5 * rCS * ( zCS * y + wCS * x ) + 0.5 * ( y * y - x * x );
00244   double bCPPlus = 1. - y ;
00245   double bCPMinus = 1. + y ;
00246   
00247   if( !m_useNewWeights )
00248     {
00249       bCF = 1. ;
00250       m_rwsCF = rCF2 ;
00251       bCS = 1. ;
00252       m_rwsCS = rCS2 ;
00253       bCPPlus = 1. ;
00254       bCPMinus = 1. ;
00255     }
00256   
00257   m_rwsCF = ( rCF2 + 0.5 * rCF * ( zCF * y - wCF * x )  + rmix ) / bCF;//though division by b not mentioned in memo it is correct
00258   m_rwsCS = ( rCS2 + 0.5 * rCS * ( zCS * y - wCS * x )  + rmix ) / bCS;
00259 
00260   
00261   // Calculate reweighting factors, normalized to the largest factor.
00262   // Note:A fake Dalitz weight is added in to just to make initial predictions
00263   // These Dalitz weights are not used in actual code 
00264   
00265   // FlavoredCF
00266   
00267   // First calculate rate factors.
00268   m_weights[ kFlavoredCF ][ kFlavoredCFBar ] =
00269       1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
00270   m_weights[ kFlavoredCF ][ kFlavoredCF ] =
00271     rmix * m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;
00272   m_weights[ kFlavoredCF ][ kFlavoredCSBar ] =
00273     1. + rCF2 * rCS2 - rCF * rCS * vCFCSMinus ;
00274   m_weights[ kFlavoredCF ][ kFlavoredCS ] =
00275     rCF2 + rCS2 - rCF * rCS * vCFCSPlus +
00276     rmix * m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
00277   
00278   // Nov 2007: correct for r2 -> RWS and BR != A^2
00279   m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
00280   m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
00281     ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
00282   m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
00283     ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
00284   m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
00285     ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
00286   
00287   m_weights[ kFlavoredCF ][ kSLPlus ] = 1. / bCF ;
00288   m_weights[ kFlavoredCF ][ kSLMinus ] = 1. / bCF ;
00289   
00290   m_weights[ kFlavoredCF ][ kCPPlus ] =
00291     ( 1. + rCF2 + rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPPlus ;
00292   m_weights[ kFlavoredCF ][ kCPMinus ] =
00293     ( 1. + rCF2 - rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPMinus ;
00294   m_weights[ kFlavoredCF ][ kDalitz ] = 1. / bCF;
00295 
00296   // FlavoredCFBar
00297   m_weights[ kFlavoredCFBar ][ kFlavoredCFBar ] = m_weights[ kFlavoredCF ][ kFlavoredCF ] ;
00298   m_weights[ kFlavoredCFBar ][ kFlavoredCS ]    = m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
00299   m_weights[ kFlavoredCFBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF ][ kFlavoredCS ] ;
00300   m_weights[ kFlavoredCFBar ][ kSLPlus ]        = 1. / bCF ;
00301   m_weights[ kFlavoredCFBar ][ kSLMinus ]       = 1. / bCF ;
00302   m_weights[ kFlavoredCFBar ][ kCPPlus ]        = m_weights[ kFlavoredCF ][ kCPPlus ] ;
00303   m_weights[ kFlavoredCFBar ][ kCPMinus ]       = m_weights[ kFlavoredCF ][ kCPMinus ] ;
00304   m_weights[ kFlavoredCFBar ][ kDalitz ]        = 1. / bCF;
00305 
00306   // FlavoredCS
00307   
00308   // First calculate rate factors.
00309   m_weights[ kFlavoredCS ][ kFlavoredCSBar ] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2 ;
00310   m_weights[ kFlavoredCS ][ kFlavoredCS ]    = rmix * m_weights[ kFlavoredCS ][ kFlavoredCSBar ] ;
00311 
00312   // Nov 2007: correct for r2 -> RWS and BR != A^2
00313   m_weights[ kFlavoredCS ][ kFlavoredCS ]    /= m_rwsCS * bCS * bCS ;
00314   m_weights[ kFlavoredCS ][ kFlavoredCSBar ] /= ( 1. + m_rwsCS * m_rwsCS ) * bCS * bCS ;
00315   
00316   m_weights[ kFlavoredCS ][ kSLPlus ]  = 1. / bCS ;
00317   m_weights[ kFlavoredCS ][ kSLMinus ] = 1. / bCS ;
00318   
00319   m_weights[ kFlavoredCS ][ kCPPlus ]  = ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus ;
00320   m_weights[ kFlavoredCS ][ kCPMinus ] = ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus ;
00321   
00322   m_weights[ kFlavoredCS ][ kDalitz ]  = 1. / bCS;
00323  
00324   // FlavoredCSBar
00325   
00326   m_weights[ kFlavoredCSBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCS ][ kFlavoredCS ] ;
00327   m_weights[ kFlavoredCSBar ][ kSLPlus ]        = 1. / bCS ;
00328   m_weights[ kFlavoredCSBar ][ kSLMinus ]       = 1. / bCS ;
00329   m_weights[ kFlavoredCSBar ][ kCPPlus ]        = m_weights[ kFlavoredCS ][ kCPPlus ] ;
00330   m_weights[ kFlavoredCSBar ][ kCPMinus ]       = m_weights[ kFlavoredCS ][ kCPMinus ] ;
00331   m_weights[ kFlavoredCSBar ][ kDalitz ]        = 1. / bCS;  
00332 
00333   // SLPlus
00334   
00335   // SLPlus/SLPlus should have rate factor = Rmix, but none of these are
00336   // generated in uncorrelated MC.
00337   m_weights[ kSLPlus ][ kSLPlus ]  = 0. ;
00338   m_weights[ kSLPlus ][ kSLMinus ] = 1. ;
00339   m_weights[ kSLPlus ][ kCPPlus ]  = 1. / bCPPlus ;
00340   m_weights[ kSLPlus ][ kCPMinus ] = 1. / bCPMinus ;
00341   m_weights[ kSLPlus ][ kDalitz ]  = 1. ;
00342   
00343   // SLMinus
00344   
00345   m_weights[ kSLMinus ][ kSLMinus ] = 0. ;
00346   m_weights[ kSLMinus ][ kCPPlus ] = 1. / bCPPlus ;
00347   m_weights[ kSLMinus ][ kCPMinus ] = 1. / bCPMinus ;
00348   m_weights[ kSLMinus ][ kDalitz ]  = 1. ;
00349   
00350   // CPPlus
00351   
00352   m_weights[ kCPPlus ][ kCPPlus ] = 0. ;
00353   m_weights[ kCPPlus ][ kCPMinus ] = 2. / bCPPlus / bCPMinus ;
00354   m_weights[ kCPPlus ][ kDalitz ]  = 1. / bCPPlus;
00355 
00356   // CPMinus
00357   
00358   m_weights[ kCPMinus ][ kCPMinus ] = 0. ;
00359   m_weights[ kCPMinus ][ kDalitz ]  = 1. / bCPMinus;
00360 
00361   //Dalitz estimate
00362   m_weights[ kDalitz ][ kDalitz ] = 1;
00363 
00364   log << MSG::FATAL << "Weight matrix" << m_weights <<  endmsg ;
00365 
00366   //Boss 6.6.2 MC Branching fractions 
00367   HepVector fractionsD0( kNDecayTypes+1, 0 ) ; 
00368   fractionsD0[ kFlavoredCF ]    = 0.531 ;  //.510(CLEO-c BR for comparison)
00369   fractionsD0[ kFlavoredCFBar ] = 0.0082 ; //.0103
00370   fractionsD0[ kFlavoredCS ]    = 0.031 ; //.0259
00371   fractionsD0[ kFlavoredCSBar ] = 0.031 ; //.0259 
00372   fractionsD0[ kSLPlus ]        = 0.125 ;  //.129
00373   fractionsD0[ kSLMinus ]       = 0. ;       
00374   fractionsD0[ kCPPlus ]        = 0.113 ;  //.123
00375   fractionsD0[ kCPMinus ]       = 0.102 ;  //.110
00376   fractionsD0[ kDalitz ]        = 0.0588 ; //.0578
00377   HepVector scales = m_weights * fractionsD0 ;
00378   log << MSG::INFO << "Integrated scale factors " <<scales << endmsg ;
00379   
00380 
00381   HepVector fractionsD0bar( kNDecayTypes+1, 0 ) ;
00382   fractionsD0bar[ kFlavoredCF ]    = 0.0082 ;
00383   fractionsD0bar[ kFlavoredCFBar ] = 0.531 ;
00384   fractionsD0bar[ kFlavoredCS ]    = 0.031 ;
00385   fractionsD0bar[ kFlavoredCSBar ] = 0.031 ;
00386   fractionsD0bar[ kSLPlus ]        = 0. ;
00387   fractionsD0bar[ kSLMinus ]       = 0.125 ;
00388   fractionsD0bar[ kCPPlus ]        = 0.113 ;
00389   fractionsD0bar[ kCPMinus ]       = 0.102 ;
00390   fractionsD0bar[ kDalitz ]        = 0.0588 ;
00391   HepVector weight_norm = fractionsD0bar.T() * scales;
00392   log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg ;
00393  
00394 
00395 
00396   //Original MC
00397   //Branching fraction predictions after weight
00398   double brCF = ( fractionsD0bar[ kFlavoredCFBar ] * scales[ kFlavoredCFBar ] ) / weight_norm[0] ;
00399   double brCS = ( fractionsD0bar[ kFlavoredCS ] * scales[ kFlavoredCS ] +  fractionsD0bar[ kFlavoredCSBar ] * scales[ kFlavoredCSBar ] ) / weight_norm[0] ;
00400   double brDCS = ( fractionsD0bar[ kFlavoredCF ] * scales[ kFlavoredCF ] ) / weight_norm[0] ;
00401   double brCPPlus = ( fractionsD0bar[ kCPPlus ] * scales[ kCPPlus ] ) / weight_norm[0] ;
00402   double brCPMinus = ( fractionsD0bar[ kCPMinus ] * scales[ kCPMinus ] ) / weight_norm[0] ;
00403 
00404   //y=0.0 dalitz example used for prediction
00405   double dalitz_y_num = -0.000177536;
00406   double dalitz_y_den = -0.0150846;
00407 
00408   double numFactCF =
00409     -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ; 
00410   double numFactCS =
00411     -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ; 
00412   double denFactCF = 0.5 * rCF2 * zCF2 ; 
00413   double denFactCS = 0.5 * rCS2 * zCS2 ; 
00414 
00415   double num_pre =
00416     brCPPlus - brCPMinus + brCF * numFactCF + 0.5 * brCS * numFactCS + dalitz_y_num; 
00417   double den_pre =
00418     1. - brCPPlus - brCPMinus - brCF * denFactCF - 0.5 * brCS * denFactCS + dalitz_y_den ; 
00419   double y_pre = num_pre / den_pre ;
00420   double num =
00421     fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] + fractionsD0bar[ kFlavoredCFBar ] * numFactCF + fractionsD0bar[ kFlavoredCS ] * numFactCS + dalitz_y_num; 
00422   double den =
00423     1. - fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] - fractionsD0bar[ kFlavoredCFBar ] * denFactCF - fractionsD0bar[ kFlavoredCS ] * denFactCS + dalitz_y_den; 
00424   double y_prematrix = num / den ;
00425   double y_test1 = 0.25 * ( ( ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) / ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) ) -
00426                            ( ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) / ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) ) ) ;
00427 
00428   log << MSG::INFO << "An Input  Y of " << m_y << endmsg ;
00429   log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg ;
00430   log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg ;
00431   log << MSG::INFO << "The CP-SL double tag method predicts y of " <<y_test1 << endmsg ; 
00432 
00433   // Renormalize
00434   m_largestWeight = 0. ;
00435   for( int i = 0 ; i < kNDecayTypes ; ++i )
00436     {
00437       for( int j = 0 ; j < kNDecayTypes ; ++j )
00438         {
00439           if( m_weights[ i ][ j ] < 0 ) m_weights[ i ][ j ] = 0; //gets rid of negative weights in early calculations
00440           if( m_weights[ i ][ j ] > m_largestWeight )
00441             {
00442               m_largestWeight = m_weights[ i ][ j ] ;
00443             }
00444         }
00445     }
00446   m_weights /= m_largestWeight ;
00447   
00448   log << MSG::INFO <<"Renormalized weight matrix " << m_weights << endmsg ;
00449 
00450 
00451   //Set up weight parameters
00452   Dalitz t(8);
00453   double D0Sum[ 10 ], D0barSum[ 10 ] ; // last index is for sum over bins
00454   TComplex D0D0barSum[ 10 ] ;
00455   int points[ 10 ] ;
00456   
00457   for( int i = 0 ; i < 10 ; ++i )
00458     {
00459       D0Sum[ i ] = 0. ;
00460       D0barSum[ i ] = 0. ;
00461       D0D0barSum[ i ] = TComplex( 0., 0. ) ;
00462       points[ i ] = 0 ;
00463     }
00464   
00465   double m2min = 0. ;
00466   double m2max = 3. ;
00467   int nsteps = 1000 ;
00468   double stepsize = ( m2max - m2min ) / ( double ) nsteps ;
00469   
00470   for( int i = 0 ; i < nsteps ; ++i )
00471     {
00472       double m2i = m2min + ( ( double ) i + 0.5 ) * stepsize ;
00473       
00474       for( int j = 0 ; j < nsteps ; ++j ) 
00475         {
00476           double m2j = m2min + ( ( double ) j + 0.5 ) * stepsize ;
00477           
00478           if( t.Point_on_DP( m2i, m2j ) )
00479             {
00480               double m2k = 1.86450*1.86450 + 0.497648*0.497648 +
00481                 0.139570*0.139570 + 0.139570*0.139570 - m2i - m2j ;
00482               
00483               TComplex D0, D0bar;
00484               if (m_dalitzModel == 0) { //Cleo model
00485                 D0 = t.CLEO_Amplitude( m2i, m2j, m2k ) ;
00486                 D0bar = t.CLEO_Amplitude( m2j, m2i, m2k ) ;}
00487               if (m_dalitzModel == 1) { //Babar model
00488                 D0 = t.Babar_Amplitude( m2i, m2j, m2k ) ;
00489                 D0bar = t.Babar_Amplitude( m2j, m2i, m2k ) ;}
00490               if (m_dalitzModel == 2) { //Belle model
00491                 D0 = t.Amplitude( m2i, m2j, m2k ) ;
00492                 D0bar = t.Amplitude( m2j, m2i, m2k ) ;}
00493               
00494               if ( j <= i ) {// lower half only
00495                 int bin = t.getBin( m2i, m2j, m2k ) ;
00496                 if( bin == -1 ) bin = 8 ;
00497                 
00498                 ++points[ bin ] ;
00499                 D0Sum[ bin ] += norm( D0 ) ;
00500                 D0barSum[ bin ] += norm( D0bar ) ;
00501                 D0D0barSum[ bin ] += D0 * conj( D0bar ) ;
00502                 
00503                 ++points[ 9 ] ;
00504                 D0Sum[ 9 ] += norm( D0 ) ;
00505                 D0barSum[ 9 ] += norm( D0bar ) ;
00506                 D0D0barSum[ 9 ] += D0 * conj( D0bar ) ;
00507               }
00508                 
00509             }
00510         }
00511     }
00512 
00513   for( int i = 0 ; i < 10 ; ++i )
00514     {
00515       double r2 = D0barSum[ i ] / D0Sum[ i ] ;
00516       double c = real( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
00517       double s = imag( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
00518       
00519       cout << "BIN " << i << ": "
00520            << points[ i ] << " pts "
00521            << r2 << " " << c << " " << s << " " << c*c+s*s
00522            << endmsg ;
00523     }
00524   
00525   log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00526   return StatusCode::SUCCESS;
00527 }

Member Data Documentation

double QCMCFilter::m_dalitzDenom [private]

Definition at line 52 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_dalitzModel [private]

Definition at line 64 of file QCMCFilter.h.

Referenced by findD0Decay(), initialize(), and QCMCFilter().

double QCMCFilter::m_dalitzNumer1 [private]

Definition at line 50 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

double QCMCFilter::m_dalitzNumer2 [private]

Definition at line 51 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

double QCMCFilter::m_deltaCF [private]

Definition at line 69 of file QCMCFilter.h.

Referenced by findD0Decay(), and initialize().

double QCMCFilter::m_deltaCS [private]

Definition at line 70 of file QCMCFilter.h.

Referenced by findD0Decay(), and initialize().

double QCMCFilter::m_dplusDminusWeight [private]

Definition at line 63 of file QCMCFilter.h.

Referenced by execute(), and QCMCFilter().

IDataProviderSvc* QCMCFilter::m_evtSvc [private]

Definition at line 30 of file QCMCFilter.h.

Referenced by execute().

bool QCMCFilter::m_invertSelection [private]

Definition at line 54 of file QCMCFilter.h.

Referenced by execute(), and QCMCFilter().

double QCMCFilter::m_largestWeight [private]

Definition at line 72 of file QCMCFilter.h.

Referenced by execute(), and initialize().

int QCMCFilter::m_nCPMinus [private]

Definition at line 43 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nCPPlus [private]

Definition at line 42 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nD0bar [private]

Definition at line 38 of file QCMCFilter.h.

Referenced by execute(), and finalize().

int QCMCFilter::m_nD0D0barDiscarded [private]

Definition at line 40 of file QCMCFilter.h.

Referenced by execute(), and finalize().

int QCMCFilter::m_nD0D0barEvents [private]

Definition at line 37 of file QCMCFilter.h.

Referenced by execute(), and finalize().

int QCMCFilter::m_nDalitz [private]

Definition at line 49 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nDpDmDiscarded [private]

Definition at line 41 of file QCMCFilter.h.

Referenced by execute(), and finalize().

int QCMCFilter::m_nDpDmEvents [private]

Definition at line 39 of file QCMCFilter.h.

Referenced by execute(), and finalize().

int QCMCFilter::m_nFlavoredCFD0 [private]

Definition at line 44 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nFlavoredCSD0 [private]

Definition at line 45 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nFlavoredDCSD0 [private]

Definition at line 46 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nSL [private]

Definition at line 47 of file QCMCFilter.h.

Referenced by finalize(), and findD0Decay().

int QCMCFilter::m_nUnknownDecays [private]

Definition at line 36 of file QCMCFilter.h.

Referenced by finalize().

int QCMCFilter::m_nUnknownEvents [private]

Definition at line 35 of file QCMCFilter.h.

Referenced by execute(), and finalize().

HepPDT::ParticleDataTable* QCMCFilter::m_particleTable [private]

Definition at line 32 of file QCMCFilter.h.

Referenced by initialize().

double QCMCFilter::m_rCF [private]

Definition at line 57 of file QCMCFilter.h.

Referenced by finalize(), findD0Decay(), initialize(), and QCMCFilter().

double QCMCFilter::m_rCS [private]

Definition at line 60 of file QCMCFilter.h.

Referenced by finalize(), findD0Decay(), initialize(), and QCMCFilter().

double QCMCFilter::m_rwsCF [private]

Definition at line 67 of file QCMCFilter.h.

Referenced by findD0Decay(), and initialize().

double QCMCFilter::m_rwsCS [private]

Definition at line 68 of file QCMCFilter.h.

Referenced by findD0Decay(), and initialize().

bool QCMCFilter::m_useNewWeights [private]

Definition at line 65 of file QCMCFilter.h.

Referenced by initialize(), and QCMCFilter().

bool QCMCFilter::m_wCFSign [private]

Definition at line 59 of file QCMCFilter.h.

Referenced by finalize(), initialize(), and QCMCFilter().

bool QCMCFilter::m_wCSSign [private]

Definition at line 62 of file QCMCFilter.h.

Referenced by finalize(), initialize(), and QCMCFilter().

double QCMCFilter::m_x [private]

Definition at line 55 of file QCMCFilter.h.

Referenced by execute(), finalize(), findD0Decay(), initialize(), and QCMCFilter().

double QCMCFilter::m_y [private]

Definition at line 56 of file QCMCFilter.h.

Referenced by execute(), initialize(), and QCMCFilter().

double QCMCFilter::m_zCF [private]

Definition at line 58 of file QCMCFilter.h.

Referenced by finalize(), initialize(), and QCMCFilter().

double QCMCFilter::m_zCS [private]

Definition at line 61 of file QCMCFilter.h.

Referenced by finalize(), initialize(), and QCMCFilter().

Generated on Tue Nov 29 23:20:44 2016 for BOSS_7.0.2 by  doxygen 1.4.7