Differential D41 Diff 202 src/LauKMatrixProdPole.cc

Changeset View

Standalone View

src/LauKMatrixProdPole.cc

Show All 22 Lines
*/		*/

/*! \file LauKMatrixProdPole.cc		/*! \file LauKMatrixProdPole.cc
\brief File containing implementation of LauKMatrixProdPole class.		\brief File containing implementation of LauKMatrixProdPole class.
*/		*/

#include "LauKMatrixProdPole.hh"		#include "LauKMatrixProdPole.hh"
#include "LauKMatrixPropagator.hh"		#include "LauKMatrixPropagator.hh"
		#include "LauResonanceMaker.hh"

#include <iostream>		#include <iostream>

ClassImp(LauKMatrixProdPole)		ClassImp(LauKMatrixProdPole)

LauKMatrixProdPole::LauKMatrixProdPole(const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,		LauKMatrixProdPole::LauKMatrixProdPole( const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,
LauKMatrixPropagator* propagator, const LauDaughters* daughters,		LauKMatrixPropagator* propagator, const LauDaughters* daughters,
Bool_t useProdAdler) :		Bool_t useProdAdler) :
LauAbsResonance(poleName, resPairAmpInt, daughters),		LauAbsResonance( poleName, resPairAmpInt, daughters, propagator->getL(propagator->getIndex()) ),
thePropagator_(propagator),		thePropagator_(propagator),
poleIndex_(poleIndex - 1), // poleIndex goes from 1 to nPoles		poleIndex_(poleIndex - 1), // poleIndex goes from 1 to nPoles
useProdAdler_(useProdAdler)		useProdAdler_(useProdAdler)
{		{

if (useProdAdler_) {		if (useProdAdler_) {
std::cout<<"Creating K matrix production pole "<<poleName<<" with poleIndex = "		std::cout <<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
<<poleIndex<<" with s-dependent production Adler zero term"<<std::endl;		<<poleIndex<<" with s-dependent production Adler zero term"<<std::endl;
} else {		} else {
std::cout<<"Creating K matrix production pole "<<poleName<<" with poleIndex = "		std::cout <<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
<<poleIndex<<" with production Adler zero factor = 1"<<std::endl;		<<poleIndex<<" with production Adler zero factor = 1"<<std::endl;
}		}

		// `Resonance' Blatt-Weisskopf factor is handled internally, but parent must be set here. For other lineshapes, LauResonanceMaker handles this.
		this->setBarrierRadii( nullptr,LauResonanceMaker::get().getParentBWFactor( propagator->getL(propagator->getIndex()), LauBlattWeisskopfFactor::BWBarrier ) );
}		}

LauKMatrixProdPole::~LauKMatrixProdPole()		LauKMatrixProdPole::~LauKMatrixProdPole()
{		{
}		}

LauComplex LauKMatrixProdPole::resAmp(Double_t mass, Double_t spinTerm)		LauComplex LauKMatrixProdPole::resAmp(const Double_t mass, const Double_t spinTerm)
{		{
std::cerr << "ERROR in LauKMatrixProdPole::resAmp : This method shouldn't get called." << std::endl;
std::cerr << " Returning zero amplitude for mass = " << mass << " and spinTerm = " << spinTerm << "." << std::endl;
return LauComplex(0.0, 0.0);
}

LauComplex LauKMatrixProdPole::amplitude(const LauKinematics* kinematics)
{

// Calculate the amplitude for the K-matrix production pole.		// Calculate the amplitude for the K-matrix production pole.
// First, make sure the K-matrix propagator is up-to-date for
// the given centre-of-mass squared value ("s") from the kinematics.
LauComplex amp(0.0, 0.0);		LauComplex amp(0.0, 0.0);

if (thePropagator_ == 0) {		if (thePropagator_ == 0) {
std::cerr << "ERROR in LauKMatrixProdPole::amplitude : The propagator is null" << std::endl;		std::cerr << "ERROR in LauKMatrixProdPole::amplitude : The propagator is null" << std::endl;
return amp;		return amp;
}		}

thePropagator_->updatePropagator(kinematics);		// Get barrier factors ('resonance' factor is already accounted for internally via propagator 'Gamma' matrix)
		Double_t fFactorB(1.0);

		const Int_t resSpin = this->getSpin();
		const Double_t pstar = this->getPstar();

		if ( resSpin > 0 ) {
		const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
		if ( parBWFactor != nullptr ) {
		johndanAuthorUnsubmitted Done Inline Actions @jback, here I assume the 0th channel is the one for which the amplitude is requested. This assumption has precedence in your code, right? johndan: @jback, here I assume the 0th channel is the one for which the amplitude is requested. This…
		jbackUnsubmitted Done Inline Actions Well, you could have the case where the first channel has L = 0, but others have L > 0. So we should check if the propagator has any non-zero L channels, perhaps by returning a cached boolean that keeps track of this condition at initialisation. Do we need to include the resBWFactor? I thought this was already in the "D" matrices (denominator terms involving q, a and R). jback: Well, you could have the case where the first channel has L = 0, but others have L > 0. So we…
		johndanAuthorUnsubmitted Done Inline Actions On the first point, I don't think it matters if other channels have L>0 does it? Why do we need to track it? On the second, I think you're spot on - the resBWFactor has been dropped as you say. johndan: On the first point, I don't think it matters if other channels have L>0 does it? Why do we need…
		jbackUnsubmitted Done Inline Actions Yes, we don't need to track this since we have changed how the angular dependence is included. jback: Yes, we don't need to track this since we have changed how the angular dependence is included.
		switch ( parBWFactor->getRestFrame() ) {
		case LauBlattWeisskopfFactor::ResonanceFrame:
		fFactorB = parBWFactor->calcFormFactor(this->getP());
		break;
		case LauBlattWeisskopfFactor::ParentFrame:
		fFactorB = parBWFactor->calcFormFactor(pstar);
		break;
		case LauBlattWeisskopfFactor::Covariant:
		{
		Double_t covFactor = this->getCovFactor();
		if ( resSpin > 2 ) {
		covFactor = TMath::Power( covFactor, 1.0/resSpin );
		} else if ( resSpin == 2 ) {
		covFactor = TMath::Sqrt( covFactor );
		}
		fFactorB = parBWFactor->calcFormFactor(pstar*covFactor);
		break;
		}
		}
		}
		}

		thePropagator_->updatePropagator(mass*mass);

// Sum the pole denominator terms over all channels j, multiplying by		// Sum the pole denominator terms over all channels j, multiplying by
// the propagator terms. Note that we do not sum over poles, since we		// the propagator terms. Note that we do not sum over poles, since we
// only want one of the production pole terms.		// only want one of the production pole terms.
Int_t nChannels = thePropagator_->getNChannels();		Int_t nChannels = thePropagator_->getNChannels();
Int_t jChannel;		Int_t jChannel;

for (jChannel = 0; jChannel < nChannels; jChannel++) {		for (jChannel = 0; jChannel < nChannels; jChannel++) {
Show All 10 Lines	LauComplex LauKMatrixProdPole::resAmp(const Double_t mass, const Double_t spinTerm)
Double_t poleDenom = thePropagator_->getPoleDenomTerm(poleIndex_);		Double_t poleDenom = thePropagator_->getPoleDenomTerm(poleIndex_);

// Include Adler zero factor if requested		// Include Adler zero factor if requested
Double_t adlerZero(1.0);		Double_t adlerZero(1.0);
if (useProdAdler_) {adlerZero = thePropagator_->getAdlerZero();}		if (useProdAdler_) {adlerZero = thePropagator_->getAdlerZero();}

amp.rescale(poleDenom*adlerZero);		amp.rescale(poleDenom*adlerZero);

		// Scale by the spin term
		Double_t scale = spinTerm;

		// Include Blatt-Weisskopf barrier factor for parent
		scale *= fFactorB;

		amp.rescale(scale);

return amp;		return amp;

}		}

		const std::vector<LauParameter*>& LauKMatrixProdPole::getFloatingParameters()
		{

		this->clearFloatingParameters();

		Int_t nChannels = thePropagator_->getNChannels();

		for (int jChannel = 0 ; jChannel < nChannels ; jChannel++)
		{
		LauParameter& par_gj_ = thePropagator_->getCouplingParameter(poleIndex_, jChannel);
		if ( !par_gj_.fixed() )
		this->addFloatingParameter( &par_gj_ );
		}

		LauParameter& par_polemasssq_ = thePropagator_->getPoleMassSqParameter(poleIndex_);
		if ( !par_polemasssq_.fixed() )
		{
		this->addFloatingParameter( &par_polemasssq_ );
		}

		return this->getParameters();

		}
		tlathamUnsubmitted Done Inline Actions Braces have gone missing tlatham: Braces have gone missing