diff --git a/config/particles.dat b/config/particles.dat index 18b1989..18e101b 100644 --- a/config/particles.dat +++ b/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 diff --git a/src/RapidParticleData.cc b/src/RapidParticleData.cc index 1b92544..91db7c9 100644 --- a/src/RapidParticleData.cc +++ b/src/RapidParticleData.cc @@ -5,6 +5,7 @@ #include #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 diff --git a/src/RapidParticleData.h b/src/RapidParticleData.h index 8247dfa..3811663 100644 --- a/src/RapidParticleData.h +++ b/src/RapidParticleData.h @@ -9,6 +9,7 @@ #include "RooRealVar.h" class RooRelBreitWigner; +class RooSubThreshold; class RooGounarisSakurai; class RapidParticle; @@ -16,6 +17,7 @@ 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 idToCT_; std::map idToMass_; diff --git a/src/RooSubThreshold.cc b/src/RooSubThreshold.cc new file mode 100644 index 0000000..c9f6414 --- /dev/null +++ b/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 + +//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); +} diff --git a/src/RooSubThreshold.h b/src/RooSubThreshold.h new file mode 100644 index 0000000..3a72c47 --- /dev/null +++ b/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