WIP - Add sub- thresold pdf

config/particles.dat

@@ -34,17 +34,17 @@ ID		part		anti		mass		width		charge		spin		shape		ctau
 333		phi		---		1.019461	0.004266	 0		1.0		RBW		-1.
 380		f'_2(1525)      ---             1.5174          0.0860           0              2.0             RBW             -1.
 3122		Lambda0		Lambda0b	1.115683	0.		 0		0.5		---		78.9
-13122		Lambda(1405)    Lambda(1405)b   1.405           0.051            0              0.5             RBW             -1.
-31220		Lambda(1520)	Lambda(1520)b	1.5195		0.0156		 0		1.5		RBW		-1.
+13122		Lambda(1405)    Lambda(1405)b   1.405           0.051            0              0.5             RBW            -1.
+31220		Lambda(1520)	Lambda(1520)b	1.5195		0.0156		 0		1.5		SubTh		-1.
 23122		Lambda(1600)    Lambda(1600)b   1.600           0.200            0              0.5             RBW             -1.
-33122		Lambda(1670)    Lambda(1670)b   1.674           0.030            0              0.5             RBW             -1. 
+33122		Lambda(1670)    Lambda(1670)b   1.674           0.030            0              0.5             RBW             -1.
 13124		Lambda(1690)    Lambda(1690)b   1.690           0.070            0              1.5             RBW             -1.
-43122		Lambda(1800)    Lambda(1800)b   1.800           0.200            0              0.5             RBW             -1. 
-53122		Lambda(1810)    Lambda(1810)b   1.790           0.110            0              0.5             RBW             -1. 
+43122		Lambda(1800)    Lambda(1800)b   1.800           0.200            0              0.5             RBW             -1.
+53122		Lambda(1810)    Lambda(1810)b   1.790           0.110            0              0.5             RBW             -1.
 3126		Lambda(1820)    Lambda(1820)b   1.820           0.080            0              2.5             RBW             -1.
 13126		Lambda(1830)    Lambda(1830)b   1.825           0.090            0              2.5             RBW             -1.
 23124		Lambda(1890)    Lambda(1890)b   1.890           0.120            0              1.5             RBW             -1.
-3128		Lambda(2100)    Lambda(2100)b   2.100           0.200            0              3.5             RBW             -1.   
+3128		Lambda(2100)    Lambda(2100)b   2.100           0.200            0              3.5             RBW             -1.
 3222		Sigma+		---		1.18937		0.		 1		0.5		---		24.04
 3212		Sigma0		Sigma0b		1.192642	0.		 0		0.5		---		-1.
 3112		Sigma-		---		1.197449	0.		-1		0.5		---		44.34

src/RapidParticleData.cc

@@ -5,6 +5,7 @@
 #include <iostream>
 
 #include "RooRelBreitWigner.h"
+#include "RooSubThreshold.h"
 #include "RooGounarisSakurai.h"
 
 #include "RapidParticle.h"
@@ -167,6 +168,8 @@ void RapidParticleData::setupMass(RapidParticle* part) {
 
 	double mA = part->daughter(0)->mass();
 	double mB = part->daughter(1)->mass();
+	double mC = part->daughter(2)->mass();
+	double mPar = part->mass();
 
 	if(idA>idB) {
 		int tmpId = idB;
@@ -199,12 +202,16 @@ void RapidParticleData::setupMass(RapidParticle* part) {
 		case RapidParticleData::GS:
 			pdf = makeGS(m, mass, width, spin, mA, mB, name);
 			break;
+		case RapidParticleData::SubTh:
+			  pdf = makeSubTh(m, mass, width, spin, mA, mB, mC, mPar, name);
+				break;
 		default:
 			std::cout << "WARNING in RapidParticleData::setupMass : unknown lineshape for " << name << "." << std::endl
 				  << "                                        : using a relativistic Breit-Wigner." << std::endl;
 		/* FALLTHRU */
 		case RapidParticleData::RelBW:
 			pdf = makeRelBW(m, mass, width, spin, mA, mB, name);
+
 	}
 	RooDataSet* massdata = pdf->generate(RooArgSet(m),100000);
 	massdata->getRange(m,mmin,mmax);
@@ -226,6 +233,7 @@ void RapidParticleData::addEntry(int id, TString name, double mass, double width
 	nameToId_[name] = id;
 	sanitisedNameToId_[sanitiseName(name)] = id;
 	if(lineshape=="RBW") idToShape_[id] = RapidParticleData::RelBW;
+	else if(lineshape =="SubTh")idToShape_[id] = RapidParticleData::SubTh;
 	else if(lineshape=="GS") idToShape_[id] = RapidParticleData::GS;
 }
 
@@ -294,6 +302,24 @@ RooRelBreitWigner* RapidParticleData::makeRelBW(RooRealVar& m, double mean, doub
 	return new RooRelBreitWigner(name,name, m,*m0,*g0,*radius,*ma,*mb,*spin);
 }
 
+
+
+RooSubThreshold* RapidParticleData::makeSubTh(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, double m3, double mPar,TString name) {
+
+	double barrierFactor=3.;//TODO
+
+	RooRealVar* m0     = new RooRealVar(name+"m0",    name+"m0",     mean);
+	RooRealVar* g0     = new RooRealVar(name+"g0",    name+"g0",    gamma);
+	RooRealVar* spin   = new RooRealVar(name+"spin",  name+"spin",  thespin);
+	RooRealVar* radius = new RooRealVar(name+"radius",name+"radius",barrierFactor); // not used
+	RooRealVar* ma     = new RooRealVar(name+"ma",    name+"ma",     m1);
+	RooRealVar* mb     = new RooRealVar(name+"mb",    name+"mb",     m2);
+	RooRealVar* mc     = new RooRealVar(name+"mc",    name+"mc",     m3);
+	RooRealVar* mParent= new RooRealVar(name+"mParent",    name+"mParent",     mPar);
+	return new RooSubThreshold(name,name, m,*m0,*g0,*radius,*ma, *mb, *mc, *mParent, *spin);
+}
+
+
 RooGounarisSakurai* RapidParticleData::makeGS(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, TString name) {
 
 	double barrierFactor=3.;//TODO

src/RapidParticleData.h

@@ -9,13 +9,15 @@
 #include "RooRealVar.h"
 
 class RooRelBreitWigner;
+class RooSubThreshold;
 class RooGounarisSakurai;
 class RapidParticle;
 
 class RapidParticleData {
 	public:
 		enum ResLineShape {
 			RelBW,
+			SubTh,
 			GS
 		};
 
@@ -71,7 +73,8 @@ class RapidParticleData {
 		TString makeUniqName(TString name);
 
 		RooRelBreitWigner* makeRelBW(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, TString name);
-		RooGounarisSakurai* makeGS(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, TString name);
+		RooSubThreshold* makeSubTh(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, double m3, double mPar, TString name);
+  	RooGounarisSakurai* makeGS(RooRealVar& m, double mean, double gamma, double thespin, double m1, double m2, TString name);
 
 		std::map<int, double> idToCT_;
 		std::map<int, double> idToMass_;

src/RooSubThreshold.cc

@@ -0,0 +1,139 @@
+/*****************************************************************************
+ * Project:
+ * Package:
+ *    File: $Id: RooSubThreshold.cc,v 1.7 2020/12/10 08:53:43
+ * Authors: Y.Amhis
+ *
+ *****************************************************************************/
+
+// -- CLASS DESCRIPTION [PDF] --
+// This is an implementation of Relativistic Breit Wigner for Spin 0,1,2 particles
+// with Blatt-Weisskopf form factors and barrier functions fort RooRarFit.
+// The meson radius is set to 3.1 GeV^-1 which is similar to the EvtGen number
+// but may not be appropriate for every meson. See Phys Rev D 72, 052002 (2005).
+//////////////////////////////////////////////////////
+//
+// BEGIN_HTML
+// This is an implementation of Relativistic Breit Wigner for Spin 0,1,2 particles
+// with Blatt-Weisskopf form factors and barrier functions for RooRarFit.
+// The meson radius is set to 3.1 GeV^-1 which is similar to the EvtGen number
+// but may not be appropriate for every meson. See Phys Rev D 72, 052002 (2005).
+// END_HTML
+//
+
+
+
+#include "Riostream.h"
+
+#include "RooSubThreshold.h"
+#include "RooAbsReal.h"
+#include "RooRealVar.h"
+#include "RooAbsReal.h"
+//#include "RooComplex.h"
+#include <complex>
+
+//ClassImp(RooSubThreshold)
+
+//------------------------------------------------------------------
+RooSubThreshold::RooSubThreshold(const char *name, const char *title,
+		RooAbsReal& _x, RooAbsReal& _mean,
+		RooAbsReal& _width, RooAbsReal& _radius,
+		RooAbsReal& _mass_a, RooAbsReal& _mass_b, RooAbsReal& _mass_c, RooAbsReal& _mass_parent,
+		RooAbsReal& _spin) :
+	RooAbsPdf(name,title),
+	x("x","Dependent",this,_x),
+	mean("mean","Mean",this,_mean),
+	width("width","Width",this,_width),
+	radius("radius","Radius",this,_radius),  // meson radius e.g. 3.1 GeV^-1 //to be updated for baryons
+	mass_a("mass_a","Mass of first daughter",this,_mass_a), // e.g. proton from Lambda*
+	mass_b("mass_b","Mass of second daugher",this,_mass_b), // e.g. kaon from Lambda*
+	mass_c("mass_c","Mass of third daughter", this, _mass_c),
+	mass_parent("mass_parent","Mass of parent particle", this, _mass_parent), //e.g. Lb
+	spin("spin","Spin",this,_spin)
+{
+}
+
+//------------------------------------------------------------------
+RooSubThreshold::RooSubThreshold(const RooSubThreshold& other,
+		const char* name) :
+	RooAbsPdf(other,name),
+	x("x",this,other.x),
+	mean("mean",this,other.mean),
+	width("width",this,other.width),
+	radius("radius",this,other.radius),
+	mass_a("mass_a",this,other.mass_a),
+	mass_b("mass_b",this,other.mass_b),
+	mass_c("mass_c",this,other.mass_c),
+	mass_parent("mass_parent",this,other.mass_parent),
+	spin("spin",this,other.spin)
+{
+}
+
+//------------------------------------------------------------------
+Double_t RooSubThreshold::evaluate() const
+{
+
+	double top = getWidth();
+	double dm2 = std::pow(x,2) - std::pow(mean,2);
+	double bottom = std::pow(dm2,2) + pow(mean*getWidth(),2);
+	return top/bottom;
+}
+//------------------------------------------------------------------
+Double_t RooSubThreshold::compute_meff() const
+{
+double minMass = mass_a+mass_b;
+double maxMass = mass_parent - mass_c ; // need to substract the mass of the 3rd particle
+double tanhTerm = std::tanh( mean - ((minMass+maxMass)/2)/(maxMass-minMass));
+double mean_eff = minMass + 1./2*(maxMass-minMass)*(1.+tanhTerm)	;
+return mean_eff;
+}
+//------------------------------------------------------------------
+Double_t RooSubThreshold::getWidth() const
+{
+	Double_t q  = getQ(x);
+
+	Double_t q0 = getQ(mean_eff); //pass here m_eff for the subthrehold
+	Double_t result(0.0);
+
+	if (q>0 && q0>0 && x>0 && mean>0) {
+		result = width * getQterm(q,q0) * (mean/x)
+			* (getFF(q) / getFF(q0));
+	}
+	return (result);
+}
+
+//------------------------------------------------------------------
+Double_t RooSubThreshold::getQterm(Double_t q, Double_t q0) const
+{
+	Double_t result = pow(q/q0,2*spin+1);
+	return (result); // result = (q/q0)^(2*spin+1);
+}
+
+//------------------------------------------------------------------
+Double_t RooSubThreshold::getFF(Double_t q) const
+{
+
+	Double_t z = q * radius;
+	Double_t result(1.0);
+
+	if (spin==1) {
+		result = 1.0/(1+z*z);
+	} else if (spin==2) {
+		result = 1.0/(9 + 3*z*z + z*z*z*z);
+	}
+
+	return (result); //square of the Blatt-Weisskopf form factor
+}
+
+//-----------------------------------------------------
+Double_t RooSubThreshold::getQ(Double_t mass) const
+{
+
+	if (mass < (mass_b+mass_a)) {return(0);}
+
+	const Double_t mDaugSumSq  = (mass_b+mass_a)*(mass_b+mass_a);
+	const Double_t mDaugDiffSq = (mass_b-mass_a)*(mass_b-mass_a);
+
+	Double_t q = sqrt((mass*mass-mDaugSumSq)*(mass*mass-mDaugDiffSq))/(2*mass);
+	return(q);
+}

src/RooSubThreshold.h

@@ -0,0 +1,58 @@
+#ifndef ROO_SUBTHRESHOLD
+#define ROO_SUBTHRESHOLD
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+
+class RooRealVar;
+
+class RooSubThreshold : public RooAbsPdf {
+
+	public:
+
+		inline RooSubThreshold() {};
+
+		RooSubThreshold(const char *name,
+				const char *title,
+				RooAbsReal& _x,
+				RooAbsReal& _mean,
+				RooAbsReal& _width,
+				RooAbsReal& _radius,
+				RooAbsReal& _mass_a,
+				RooAbsReal& _mass_b,
+				RooAbsReal& _mass_c,
+				RooAbsReal& _mass_parent,
+				RooAbsReal& _spin);
+		RooSubThreshold(const RooSubThreshold& other,
+				const char* name=0) ;
+		virtual TObject* clone(const char* newname) const {
+			return new RooSubThreshold(*this,newname);
+		}
+		inline virtual ~RooSubThreshold() { }
+
+	protected:
+
+		RooRealProxy x ;      // observable (mass)
+		RooRealProxy mean ;    // mean
+		RooRealProxy width ;  // width
+		RooRealProxy radius ; // meson radius parameter
+		RooRealProxy mass_a ; // mass of first daughter from decay
+		RooRealProxy mass_b ; // mass of second daughter from decay
+		RooRealProxy mass_c ; // mass of third daughter from decay
+		RooRealProxy mass_parent ; // mass of parent particle
+		RooRealProxy spin ; // spin (= 0,1,2) //need to upgrade to baryons
+    double mean_eff; // will be needed for computation of effective mass
+		Double_t evaluate() const ;
+
+
+	private:
+		Double_t getWidth() const;
+		Double_t getQterm(Double_t q, Double_t q0) const;
+		Double_t getFF(Double_t q) const;
+		Double_t getQ(Double_t mass) const;
+	  Double_t compute_meff() const;
+
+
+};
+
+#endif