#include <TimeDependentGenerator.h>

Classes
class	GenTimeEvent

Public Types
typedef std::pair< std::complex< double >, std::complex< double > >	AmpPair

Public Member Functions
	TimeDependentGenerator (const DalitzEventPattern &pattern, double width, double deltam, double deltagamma, double qoverp, double phi, TRandom3 rndm, TH1F h_efficiency=NULL, float resWidth=0.05, bool addExpEffects=false, Eff3piSymmetric *sEfficiency=NULL)

	TimeDependentGenerator (MINT::counted_ptr< FitAmpSum > model, MINT::counted_ptr< FitAmpSum > cpmodel, double width, double deltam, double deltagamma, double qoverp, double phi, TRandom3 , TH1F h_efficiency=NULL, float resWidth=0.05, bool addExpEffects=false, Eff3piSymmetric *sEfficiency=NULL)

std::pair< double, double >	generate_decay_time () const
	Generate a decay time, optionally including experimental effects. More...

MINT::counted_ptr< IDalitzEvent >	generate_event ()
	Generate a flavour, decay time and Dalitz event. More...

MINT::counted_ptr< IDalitzEvent >	generate_event (const double s13, const double s23)
	Generate a 3-body event at a given point in phase space. More...

bool	accept_or_reject (const int tag, const double decaytime, IDalitzEvent &)
	Decide whether the generated event should be accepted or rejected. More...

double	get_scale () const

float	get_gen_efficiency () const

double	pdf_value (int, double, IDalitzEvent &)

double	pdf_value (IDalitzEvent &)

double	envelope_value (const double, IDalitzEvent &)

AmpPair	amplitude_coefficients (const int tag, const double decaytime)

TH1F	draw_pdf_vs_time (IDalitzEvent &, unsigned, float, float, const std::string &name="pdf_vs_time")

TH1F	draw_envelope_vs_time (IDalitzEvent &, unsigned, float, float, const std::string &name="envelope_vs_time")

Static Public Member Functions
static DalitzEventPattern	anti (DalitzEventPattern pat)
	Take the CP conjugate of the head of the decay pattern. More...

Private Member Functions
void	init ()

Private Attributes
TRandom3 *	m_rndm

const DalitzEventPattern	m_pattern

const DalitzEventPattern	m_cppattern

MINT::counted_ptr< FitAmpSum >	m_model

MINT::counted_ptr< FitAmpSum >	m_cpmodel

FitAmpIncoherentSum	m_bothmodel

SignalGenerator	m_generator

const double	m_width

const double	m_deltam

const double	m_deltagamma

const std::complex< double >	m_qoverp

double	m_scale

unsigned	m_ngen

unsigned	m_naccept

TH1F *	m_h_efficiency

TSpline3	m_efficiencyFit

float	m_resWidth

bool	m_addExpEffects

Eff3piSymmetric *	m_sEfficiency

Detailed Description

Definition at line 20 of file TimeDependentGenerator.h.

Member Typedef Documentation

◆ AmpPair

typedef std::pair<std::complex<double>, std::complex<double> > TimeDependentGenerator::AmpPair

Definition at line 106 of file TimeDependentGenerator.h.

Constructor & Destructor Documentation

◆ TimeDependentGenerator() [1/2]

TimeDependentGenerator::TimeDependentGenerator	(	const DalitzEventPattern &	pattern,
		double	width,
		double	deltam,
		double	deltagamma,
		double	qoverp,
		double	phi,
		TRandom3 *	rndm,
		TH1F *	h_efficiency = `NULL`,
		float	resWidth = `0.05`,
		bool	addExpEffects = `false`,
		Eff3piSymmetric *	sEfficiency = `NULL`
	)

Constructor, takes: pattern : The event pattern to be used (the CP conjugate will automatically be added). width : the decay width in 1/ps. deltam : the delta-mass in 1/ps. deltagamma : the delta-gamma in 1/ps. qoverp : the magnitude of q/p. phi : the phase of q/p. rndm : The random number generator to use. h_efficiency : (optional) histogram to which efficiency plot will be fitted resWidth : the width of the Gaussian decay-time resolution to apply addExpEffects : whether to add efficiency and resolution to the decay time.

Definition at line 70 of file TimeDependentGenerator.cpp.

                                                                                                                  :
   m_rndm(rndm),
   m_pattern(pattern),
   m_cppattern(anti(pattern)),
   m_model(new FitAmpSum(pattern)),
   m_cpmodel(new FitAmpSum(m_cppattern)),
   m_bothmodel(*m_model),
   m_generator(m_pattern, &m_bothmodel),
   m_width(width),
   m_deltam(deltam),
   m_deltagamma(deltagamma),
   m_qoverp(polar(qoverp, phi)),
   m_scale(1.),
   m_ngen(0),
   m_naccept(0),
   m_h_efficiency(h_efficiency),
   m_efficiencyFit(),
   m_resWidth(resWidth),
   m_addExpEffects(addExpEffects),
   m_sEfficiency(sEfficiency)
 {
   init() ;
 }

◆ TimeDependentGenerator() [2/2]

TimeDependentGenerator::TimeDependentGenerator	(	MINT::counted_ptr< FitAmpSum >	model,
		MINT::counted_ptr< FitAmpSum >	cpmodel,
		double	width,
		double	deltam,
		double	deltagamma,
		double	qoverp,
		double	phi,
		TRandom3 *	rndm,
		TH1F *	h_efficiency = `NULL`,
		float	resWidth = `0.05`,
		bool	addExpEffects = `false`,
		Eff3piSymmetric *	sEfficiency = `NULL`
	)

Constructor, takes: model : the amplitude model for the decay cpmodel : the amplitude model for the CP conjugate decay width : the decay width in 1/ps. deltam : the delta-mass in 1/ps. deltagamma : the delta-gamma in 1/ps. qoverp : the magnitude of q/p. phi : the phase of q/p. rndm : The random number generator to use. h_efficiency : (optional) histogram to which efficiency plot will be fitted resWidth : the width of the Gaussian decay-time resolution to apply addExpEffects : whether to add efficiency and resolution to the decay time.

Definition at line 98 of file TimeDependentGenerator.cpp.

                                                                                                                  :
   m_rndm(rndm),
   m_pattern(model->getAmpPtr(0)->getTreePattern()),
   m_cppattern(cpmodel->getAmpPtr(0)->getTreePattern()),
   m_model(model),
   m_cpmodel(cpmodel),
   m_bothmodel(*m_model),
   m_generator(m_pattern, &m_bothmodel),
   m_width(width),
   m_deltam(deltam),
   m_deltagamma(deltagamma),
   m_qoverp(polar(qoverp, phi)),
   m_scale(1.),
   m_ngen(0),
   m_naccept(0),
   m_h_efficiency(h_efficiency),
   m_efficiencyFit(),
   m_resWidth(resWidth),
   m_addExpEffects(addExpEffects),
   m_sEfficiency(sEfficiency)
 {
   init() ;
 }

Member Function Documentation

◆ accept_or_reject()

bool TimeDependentGenerator::accept_or_reject	(	const int	tag,
		const double	decaytime,
		IDalitzEvent &	evt
	)

Decide whether the generated event should be accepted or rejected.

Definition at line 236 of file TimeDependentGenerator.cpp.

                                                                                                       {
   ++m_ngen ;
 
   double pdfval = pdf_value(tag, decaytime, evt) ;
 
   double maxval = envelope_value(decaytime, evt) ;
   if(pdfval > maxval){
     cout << "pdfval " << pdfval << " maxval " << maxval << endl ;
     //throw out_of_range("pdfval > maxval. That shouldn't happen!") ;
     m_scale *= pdfval / maxval ;
     cout << "scale " << m_scale << endl ;
     ++m_naccept ;
     return true ;
   }
   if(m_rndm->Rndm() * maxval < pdfval){
     ++m_naccept ;
     return true ;
   }
   return false ;
 }

◆ amplitude_coefficients()

TimeDependentGenerator::AmpPair TimeDependentGenerator::amplitude_coefficients	(	const int	tag,
		const double	decaytime
	)

Get the coefficients of the amplitudes for the produced flavour and the mixed flavour given the tag and decay time.

Definition at line 296 of file TimeDependentGenerator.cpp.

                                                                                     {
   double coeff = 0.5 * exp(-decaytime * 0.5 * (m_width + 0.5 * m_deltagamma)) ;
   complex<double> expterm = exp(complex<double>(0.5 * m_deltagamma * decaytime, m_deltam * decaytime)) ;
   complex<double> plusterm = 1. + expterm ;
   complex<double> minusterm = 1. - expterm ;
   complex<double> coeffprod = coeff * plusterm ;
   complex<double> coeffmix = pow(m_qoverp, tag) * coeff * minusterm ;
   return AmpPair(coeffprod, coeffmix) ;
 }

◆ anti()

DalitzEventPattern TimeDependentGenerator::anti ( DalitzEventPattern pat )

static

Take the CP conjugate of the head of the decay pattern.

Definition at line 65 of file TimeDependentGenerator.cpp.

                                                                       {
   pat[0].antiThis() ;
   return pat ;
 }

◆ draw_envelope_vs_time()

TH1F TimeDependentGenerator::draw_envelope_vs_time	(	IDalitzEvent &	evt,
		unsigned	nbins,
		float	tmin,
		float	tmax,
		const std::string &	name = `"envelope_vs_time"`
	)

Draw the envelope value vs time at the given point in phase space.

Definition at line 319 of file TimeDependentGenerator.cpp.

                                                   {
   TH1F histo(name.c_str(), "", nbins, tmin, tmax) ;
   for(unsigned i = 1 ; i < nbins + 1 ; ++i){
     double time = histo.GetXaxis()->GetBinCenter(i) ;
     double val = envelope_value(time, evt) ;
     histo.SetBinContent(i, val) ;
     histo.SetBinError(i, 0.) ;
   }
   return histo ;
 }

◆ draw_pdf_vs_time()

TH1F TimeDependentGenerator::draw_pdf_vs_time	(	IDalitzEvent &	evt,
		unsigned	nbins,
		float	tmin,
		float	tmax,
		const std::string &	name = `"pdf_vs_time"`
	)

Draw the PDF value vs time at the given point in phase space (assumes the event has a tag).

Definition at line 306 of file TimeDependentGenerator.cpp.

                                                  {
   TH1F histo(name.c_str(), "", nbins, tmin, tmax) ;
   int tag = int(evt.getValueFromVector(GenTimeEvent::ITAG)) ;
   for(unsigned i = 1 ; i < nbins + 1 ; ++i){
     double time = histo.GetXaxis()->GetBinCenter(i) ;
     double val = pdf_value(tag, time, evt) ;
     histo.SetBinContent(i, val) ;
     histo.SetBinError(i, 0.) ;
   }
   return histo ;
 }

◆ envelope_value()

double TimeDependentGenerator::envelope_value	(	const double	decaytime,
		IDalitzEvent &	evt
	)

Get the value of the incoherent sum PDF that's used as an envelope to generate events at the given point in phase space and time.

Definition at line 232 of file TimeDependentGenerator.cpp.

                                                                                        {
   return m_bothmodel.Prob(evt) * m_scale * exp(-decaytime * m_width) * m_width ;
 }

◆ generate_decay_time()

pair< double, double > TimeDependentGenerator::generate_decay_time ( ) const

Generate a decay time, optionally including experimental effects.

Definition at line 143 of file TimeDependentGenerator.cpp.

                                                                        {
   double decaytime = m_rndm->Exp(1./m_width) ;
   double smeareddecaytime = decaytime ;
   
   if ( (m_h_efficiency != NULL) && (m_addExpEffects) ){
     int i = 0 ;
     float efficiency = 0 ;
     int maxiter = 100000 ;
     
     smeareddecaytime =  decaytime + m_rndm->Gaus(0, m_resWidth);
 
     while(true){
       if( smeareddecaytime > m_efficiencyFit.GetXmax() ){
         efficiency = m_efficiencyFit.Eval( m_efficiencyFit.GetXmax() ) ;
       } 
       else if( smeareddecaytime > m_efficiencyFit.GetXmin() ){
     efficiency = m_efficiencyFit.Eval( smeareddecaytime ) ;
       }
       else{
         efficiency = m_efficiencyFit.Eval( m_efficiencyFit.GetXmin() ) ;
       }      
 
       if(m_rndm->Rndm() < efficiency){
         break;
       }
  
       decaytime = m_rndm->Exp(1./m_width) ;
       smeareddecaytime =  decaytime + m_rndm->Gaus(0, m_resWidth);
 
       i += 1 ;
       if( i > maxiter ){ 
         cout << "WARNING: Decay time generation limit exceeded." << endl ; 
         break ;
       }
     }
   }
   return pair<double, double>(decaytime, smeareddecaytime) ;
 }

◆ generate_event() [1/2]

MINT::counted_ptr< IDalitzEvent > TimeDependentGenerator::generate_event ( )

Generate a flavour, decay time and Dalitz event.

Definition at line 207 of file TimeDependentGenerator.cpp.

                                                                      {
   int tag(0) ;
   double decaytime(0.) ;
   double smeareddecaytime(0.) ;
 
   MINT::counted_ptr<IDalitzEvent> evt(0) ;
   
   unsigned i = 0 ;
   unsigned maxattempts = 100000 ;
   // Accept/reject using the incoherent sum in m_bothmodel as an envelope to generate
   // events.
   while(true) {
     tag = m_rndm->Rndm() > 0.5 ? 1 : -1 ;
     tie(decaytime, smeareddecaytime) = generate_decay_time() ;
     evt = m_generator.newEvent() ;
     if(accept_or_reject(tag, decaytime, *evt))
       break ;
     i += 1 ;
     if( i >= maxattempts)
       throw range_error("Failed to generate an event in 100000 attempts!") ;
   }
 
   return MINT::counted_ptr<IDalitzEvent>(new GenTimeEvent(*evt, tag, decaytime, smeareddecaytime)) ;
 }

◆ generate_event() [2/2]

MINT::counted_ptr< IDalitzEvent > TimeDependentGenerator::generate_event	(	const double	s13,
		const double	s23
	)

Generate a 3-body event at a given point in phase space.

Definition at line 182 of file TimeDependentGenerator.cpp.

                                                                                                        {
   int tag(0) ;
   double decaytime(0.) ;
   double smeareddecaytime(0.) ;
 
   MINT::counted_ptr<IDalitzEvent> evt(0) ;
   
   unsigned i = 0 ;
   unsigned maxattempts = 100000 ;
   // Accept/reject using the incoherent sum in m_bothmodel as an envelope to generate
   // events.
   while(true) {
     tag = m_rndm->Rndm() > 0.5 ? 1 : -1 ;
     tie(decaytime, smeareddecaytime) = generate_decay_time() ;
     evt = MINT::counted_ptr<IDalitzEvent>(new DalitzEvent(tag == 1 ? m_pattern : m_cppattern, s13, s23)) ;
     if(accept_or_reject(tag, decaytime, *evt))
       break ;
     i += 1 ;
     if( i >= maxattempts)
       throw range_error("Failed to generate an event in 100000 attempts!") ;
   }
 
   return MINT::counted_ptr<IDalitzEvent>(new GenTimeEvent(*evt, tag, decaytime, smeareddecaytime)) ;
 }

◆ get_gen_efficiency()

float TimeDependentGenerator::get_gen_efficiency ( ) const

Definition at line 261 of file TimeDependentGenerator.cpp.

                                                        {
   return m_ngen > 0 ? float(m_naccept)/m_ngen : 0. ;
 }

◆ get_scale()

double TimeDependentGenerator::get_scale ( ) const

Definition at line 257 of file TimeDependentGenerator.cpp.

                                                {
   return m_scale ;
 }

◆ init()

void TimeDependentGenerator::init ( )

private

Definition at line 127 of file TimeDependentGenerator.cpp.

                                   {
   m_bothmodel.addAsList(*m_cpmodel) ;
   if( m_h_efficiency != NULL ){
     m_efficiencyFit = TSpline3(m_h_efficiency) ;
   }
   if( m_sEfficiency != NULL ){
     m_model->setEfficiency(m_sEfficiency);
     m_cpmodel->setEfficiency(m_sEfficiency);
   }
   // Ensure the models are fully initialised.
   DalitzEvent evt(m_pattern, m_rndm) ;
   m_model->getVal(evt) ;
   m_cpmodel->getVal(evt) ;
   m_bothmodel.getVal(evt) ;
 }

◆ pdf_value() [1/2]

double TimeDependentGenerator::pdf_value	(	int	tag,
		double	decaytime,
		IDalitzEvent &	evt
	)

Get the value of the PDF given the tag, decay time and point in phase space

Definition at line 265 of file TimeDependentGenerator.cpp.

                                                                                      {
   complex<double> Ap = m_model->ComplexVal(evt) ;
   complex<double> Am = m_cpmodel->ComplexVal(evt) ;
   if(tag == -1)
     swap(Ap, Am) ;
   // This is gives the same values as the below using the amplitude coefficients.
   /*double magAp = norm(Ap) ;
   double magAm = norm(Am) ;
   complex<double> crossterm = pow(m_qoverp, tag) * conj(Ap) * Am ;
   double magqoverp = pow(norm(m_qoverp), tag) ;
   double deltamt = m_deltam * decaytime ;
   double halfdgammat = 0.5 * m_deltagamma * decaytime ;
   double pdfval = 0.5 * exp(-m_width * decaytime) 
   * ((magAp + magqoverp * magAm) * cosh(halfdgammat)
   + (magAp - magqoverp * magAm) * cos(deltamt)
   - 2 * crossterm.real() * sinh(halfdgammat)
   + 2 * crossterm.imag() * sin(deltamt)) ;
   return pdfval ;*/
   AmpPair coeffs = amplitude_coefficients(tag, decaytime) ;
   complex<double> Amp = coeffs.first * Ap + coeffs.second * Am ;
   double ampnorm = norm(Amp) ;
   return ampnorm ;
 }

◆ pdf_value() [2/2]

double TimeDependentGenerator::pdf_value ( IDalitzEvent & evt )

Get the value of the PDF, assuming the tag & decay time are stored in the DalitzEvent like a GenTimeEvent

Definition at line 289 of file TimeDependentGenerator.cpp.

                                                           {
   return pdf_value(evt.getValueFromVector(GenTimeEvent::ITAG),
            evt.getValueFromVector(GenTimeEvent::IDECAYTIME),
            evt) ;
 }

Member Data Documentation

◆ m_addExpEffects

bool TimeDependentGenerator::m_addExpEffects

private

Definition at line 140 of file TimeDependentGenerator.h.

◆ m_bothmodel

FitAmpIncoherentSum TimeDependentGenerator::m_bothmodel

private

Definition at line 125 of file TimeDependentGenerator.h.

◆ m_cpmodel

MINT::counted_ptr<FitAmpSum> TimeDependentGenerator::m_cpmodel

private

Definition at line 124 of file TimeDependentGenerator.h.

◆ m_cppattern

const DalitzEventPattern TimeDependentGenerator::m_cppattern

private

Definition at line 122 of file TimeDependentGenerator.h.

◆ m_deltagamma

const double TimeDependentGenerator::m_deltagamma

private

Definition at line 130 of file TimeDependentGenerator.h.

◆ m_deltam

const double TimeDependentGenerator::m_deltam

private

Definition at line 129 of file TimeDependentGenerator.h.

◆ m_efficiencyFit

TSpline3 TimeDependentGenerator::m_efficiencyFit

private

Definition at line 137 of file TimeDependentGenerator.h.

◆ m_generator

SignalGenerator TimeDependentGenerator::m_generator

private

Definition at line 126 of file TimeDependentGenerator.h.

◆ m_h_efficiency

TH1F* TimeDependentGenerator::m_h_efficiency

private

Definition at line 136 of file TimeDependentGenerator.h.

◆ m_model

MINT::counted_ptr<FitAmpSum> TimeDependentGenerator::m_model

private

Definition at line 123 of file TimeDependentGenerator.h.

◆ m_naccept

unsigned TimeDependentGenerator::m_naccept

private

Definition at line 134 of file TimeDependentGenerator.h.

◆ m_ngen

unsigned TimeDependentGenerator::m_ngen

private

Definition at line 133 of file TimeDependentGenerator.h.

◆ m_pattern

const DalitzEventPattern TimeDependentGenerator::m_pattern

private

Definition at line 121 of file TimeDependentGenerator.h.

◆ m_qoverp

const std::complex<double> TimeDependentGenerator::m_qoverp

private

Definition at line 131 of file TimeDependentGenerator.h.

◆ m_resWidth

float TimeDependentGenerator::m_resWidth

private

Definition at line 139 of file TimeDependentGenerator.h.

◆ m_rndm

TRandom3* TimeDependentGenerator::m_rndm

private

Definition at line 120 of file TimeDependentGenerator.h.

◆ m_scale

double TimeDependentGenerator::m_scale

private

Definition at line 132 of file TimeDependentGenerator.h.

◆ m_sEfficiency

Eff3piSymmetric* TimeDependentGenerator::m_sEfficiency

private

Definition at line 142 of file TimeDependentGenerator.h.

◆ m_width

const double TimeDependentGenerator::m_width

private

Definition at line 128 of file TimeDependentGenerator.h.

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

Mint/TimeDependentGenerator.h
src/Utils/TimeDependentGenerator.cpp

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ AmpPair

Constructor & Destructor Documentation

◆ TimeDependentGenerator() [1/2]

◆ TimeDependentGenerator() [2/2]

Member Function Documentation

◆ accept_or_reject()

◆ amplitude_coefficients()

◆ anti()

◆ draw_envelope_vs_time()

◆ draw_pdf_vs_time()

◆ envelope_value()

◆ generate_decay_time()

◆ generate_event() [1/2]

◆ generate_event() [2/2]

◆ get_gen_efficiency()

◆ get_scale()

◆ init()

◆ pdf_value() [1/2]

◆ pdf_value() [2/2]

Member Data Documentation

◆ m_addExpEffects

◆ m_bothmodel

◆ m_cpmodel

◆ m_cppattern

◆ m_deltagamma

◆ m_deltam

◆ m_efficiencyFit

◆ m_generator

◆ m_h_efficiency

◆ m_model

◆ m_naccept

◆ m_ngen

◆ m_pattern

◆ m_qoverp

◆ m_resWidth

◆ m_rndm

◆ m_scale

◆ m_sEfficiency

◆ m_width