laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauGounarisSakuraiRes.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 //****************************************************************************
16 // Class for defining the relativistic Breit-Wigner resonance model, which
17 // includes the use of Blatt-Weisskopf barrier factors.
18 //****************************************************************************
19 
20 // --CLASS DESCRIPTION [MODEL] --
21 // Class for defining the relativistic Breit-Wigner resonance model, with
22 // Blatt-Weisskopf barrier factors.
23 
24 #include <iostream>
25 
26 #include "LauConstants.hh"
27 #include "LauGounarisSakuraiRes.hh"
28 
29 ClassImp(LauGounarisSakuraiRes)
30 
31 
32 LauGounarisSakuraiRes::LauGounarisSakuraiRes(TString resName, Double_t resMass, Double_t resWidth, Int_t resSpin,
33  Int_t resCharge, Int_t resPairAmpInt, const LauDaughters* daughters) :
34  LauAbsResonance(resName, resMass, resWidth, resSpin, resCharge, resPairAmpInt, daughters),
35  q0_(0.0),
36  p0_(0.0),
37  pstar0_(0.0),
38  resMassSq_(0.0),
39  mDaugSum_(0.0),
40  mDaugSumSq_(0.0),
41  mDaugDiff_(0.0),
42  mDaugDiffSq_(0.0),
43  mParentSq_(0.0),
44  mBachSq_(0.0),
45  h0_(0.0),
46  dhdm0_(0.0),
47  d_(0.0),
48  resR_(4.0),
49  parR_(4.0),
50  resRSq_(16.0),
51  parRSq_(16.0),
52  FR0_(1.0),
53  FB0_(1.0),
54  barrierType_(LauAbsResonance::BWPrimeBarrier)
55 {
56 }
57 
59 {
60 }
61 
63 {
64  // Set-up various constants. This must be called again if the mass/width/spin
65  // of a resonance changes...
66 
67  Int_t resSpin = this->getSpin();
68  Double_t resMass = this->getMass();
69  Double_t massDaug1 = this->getMassDaug1();
70  Double_t massDaug2 = this->getMassDaug2();
71  Double_t massBachelor = this->getMassBachelor();
72  Double_t massParent = this->getMassParent();
73 
74  // Check that the spin is 1
75  if (resSpin != 1) {
76  std::cerr << "WARNING in LauGounarisSakuraiRes::initialise : Resonance spin is != 1. This lineshape is for the rho(770), setting the spin to 1." << std::endl;
77  this->changeResonance( -1.0, -1.0, 1 );
78  resSpin = this->getSpin();
79  }
80 
81  // Create the mass squares, sums, differences etc.
82  resMassSq_ = resMass*resMass;
83  mDaugSum_ = massDaug1 + massDaug2;
85  mDaugDiff_ = massDaug1 - massDaug2;
87  mParentSq_ = massParent*massParent;
88  mBachSq_ = massBachelor*massBachelor;
89 
90  // Decay momentum of either daughter in the resonance rest frame
91  // when resonance mass = rest-mass value, m_0 (PDG value)
92  Double_t term1 = resMassSq_ - mDaugSumSq_;
93  Double_t term2 = resMassSq_ - mDaugDiffSq_;
94  Double_t term12 = term1*term2;
95  if (term12 > 0.0) {
96  q0_ = TMath::Sqrt(term12)/(2.0*resMass);
97  } else {
98  q0_ = 0.0;
99  }
100 
101  // Momentum of the bachelor particle in the resonance rest frame
102  // when resonance mass = rest-mass value, m_0 (PDG value)
103  Double_t eBach = (mParentSq_ - resMassSq_ - mBachSq_)/(2.0*resMass);
104  Double_t termBach = eBach*eBach - mBachSq_;
105  if ( eBach<0.0 || termBach<0.0 ) {
106  p0_ = 0.0;
107  } else {
108  p0_ = TMath::Sqrt( termBach );
109  }
110 
111  // Momentum of the bachelor particle in the parent rest frame
112  // when resonance mass = rest-mass value, m_0 (PDG value)
113  Double_t eStarBach = (mParentSq_ + mBachSq_ - resMassSq_)/(2.0*massParent);
114  Double_t termStarBach = eStarBach*eStarBach - mBachSq_;
115  if ( eStarBach<0.0 || termStarBach<0.0 ) {
116  pstar0_ = 0.0;
117  } else {
118  pstar0_ = TMath::Sqrt( termStarBach );
119  }
120 
121  // Blatt-Weisskopf barrier factor constant: z = q^2*radius^2
122  // Calculate the Blatt-Weisskopf form factor for the case when m = m_0
124 
125  // Calculate the extra things needed by the G-S shape
126  h0_ = 2.0*LauConstants::invPi * q0_/resMass * TMath::Log((resMass + 2.0*q0_)/(2.0*LauConstants::mPi));
127  dhdm0_ = h0_ * (1.0/(8.0*q0_*q0_) - 1.0/(2.0*resMassSq_)) + 1.0/(LauConstants::twoPi*resMassSq_);
129  * TMath::Log((resMass + 2.0*q0_)/(2.0*LauConstants::mPi))
130  + resMass/(LauConstants::twoPi*q0_)
132 }
133 
134 void LauGounarisSakuraiRes::setBarrierRadii(Double_t resRadius, Double_t parRadius, LauAbsResonance::BarrierType type)
135 {
136  // Reset the Blatt-Weisskopf barrier radius for the resonance and its parent
137  resR_ = resRadius;
138  parR_ = parRadius;
139 
140  resRSq_ = resRadius*resRadius;
141  parRSq_ = parRadius*parRadius;
142 
143  barrierType_ = type;
144 
145  // Recalculate the Blatt-Weisskopf form factor for the case when m = m_0
146  Double_t zR0 = q0_*q0_*resRSq_;
147  Double_t zB0 = pstar0_*pstar0_*parRSq_;
148  if ( ( type == LauAbsResonance::BWPrimeBarrier ) || ( type == LauAbsResonance::ExpBarrier ) ) {
149  FR0_ = (resR_==0.0) ? 1.0 : this->calcFFactor(zR0);
150  FB0_ = (parR_==0.0) ? 1.0 : this->calcFFactor(zB0);
151  }
152 }
153 
155 {
156  // Calculate the requested form factor for the resonance, given the z value
157  Double_t fFactor(1.0);
159  fFactor = TMath::Sqrt(2.0*z/(z + 1.0));
161  fFactor = TMath::Sqrt(1.0/(z + 1.0));
162  } else if ( barrierType_ == LauAbsResonance::ExpBarrier ) {
163  fFactor = TMath::Exp( -TMath::Sqrt(z) );
164  }
165  return fFactor;
166 }
167 
168 LauComplex LauGounarisSakuraiRes::resAmp(Double_t mass, Double_t spinTerm)
169 {
170  // This function returns the complex dynamical amplitude for a Breit-Wigner resonance,
171  // given the invariant mass and cos(helicity) values.
172 
173  LauComplex resAmplitude(0.0, 0.0);
174 
175  if (mass < 1e-10) {
176  std::cerr << "WARNING in LauGounarisSakuraiRes::amplitude : mass < 1e-10." << std::endl;
177  return LauComplex(0.0, 0.0);
178  } else if (q0_ < 1e-30) {
179  return LauComplex(0.0, 0.0);
180  }
181 
182  // Calculate the width of the resonance (as a function of mass)
183  // First, calculate the various form factors.
184  // NB
185  // q is the momentum of either daughter in the resonance rest-frame,
186  // p is the momentum of the bachelor in the resonance rest-frame,
187  // pstar is the momentum of the bachelor in the parent rest-frame.
188  // These quantities have been calculate in LauAbsResonance::amplitude(...)
189 
190  Double_t resMass = this->getMass();
191  Double_t resWidth = this->getWidth();
192  Double_t q = this->getQ();
193  //Double_t p = this->getP();
194  Double_t pstar = this->getPstar();
195 
196  Double_t zR = q*q*resRSq_;
197  Double_t zB = pstar*pstar*parRSq_;
198  Double_t fFactorR = this->calcFFactor(zR);
199  Double_t fFactorB = this->calcFFactor(zB);
200  Double_t fFactorRRatio = fFactorR/FR0_;
201  Double_t fFactorBRatio = fFactorB/FB0_;
202 
203  Double_t qRatio = q/q0_;
204  Double_t qTerm = qRatio*qRatio*qRatio;
205 
206  Double_t totWidth = resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
207 
208  Double_t massSq = mass*mass;
209  Double_t massSqTerm = resMassSq_ - massSq;
210 
211  Double_t h = 2.0*LauConstants::invPi * q/mass * TMath::Log((mass + 2.0*q)/(2.0*LauConstants::mPi));
212  Double_t f = totWidth * resMassSq_/(q0_*q0_*q0_) * (q*q * (h - h0_) + massSqTerm * q0_*q0_ * dhdm0_);
213 
214  // Compute the complex amplitude
215  resAmplitude = LauComplex(massSqTerm + f, resMass*totWidth);
216 
217  // Scale by the denominator factor, as well as the spin term and Blatt-Weisskopf factors
218  Double_t numerFactor = fFactorRRatio*fFactorBRatio*spinTerm*(1 + d_ * resWidth/resMass);
219  Double_t denomFactor = (massSqTerm + f)*(massSqTerm + f) + resMassSq_*totWidth*totWidth;
220  resAmplitude.rescale(numerFactor/denomFactor);
221 
222  return resAmplitude;
223 }
224 
Double_t getQ() const
Get the current value of the daughter momentum in the resonance rest frame.
Double_t parRSq_
Square of the radius of the barrier for parent decay.
Double_t q0_
Momentum of the daughters in the resonance rest frame (at pole mass)
Double_t getMassBachelor() const
Get the mass of the bachelor daughter.
Double_t getMass() const
Get the mass of the resonance.
const Double_t twoPi
Two times Pi.
Definition: LauConstants.hh:93
Double_t mBachSq_
Square of the bachelor mass.
Double_t parR_
Radius of the barrier for parent decay.
Double_t getMassParent() const
Get the parent particle mass.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
Class for defininf the Gounaris-Sakurai resonance model.
Double_t calcFFactor(Double_t z)
Calculate the form factor for the resonance.
Double_t getPstar() const
Get the current value of the bachelor momentum in the parent rest frame.
Double_t h0_
Extra parameter required by GS shape.
File containing declaration of LauGounarisSakuraiRes class.
void setBarrierRadii(Double_t resRadius, Double_t parRadius, LauAbsResonance::BarrierType type)
Set the form factor model and parameters.
Double_t getMassDaug1() const
Get the mass of daughter 1.
virtual ~LauGounarisSakuraiRes()
Destructor.
virtual void initialise()
Initialise the model.
LauAbsResonance::BarrierType barrierType_
Model to be used for the form factor.
Double_t getMassDaug2() const
Get the mass of daughter 2.
Double_t dhdm0_
Extra parameter required by GS shape.
Double_t mParentSq_
Square of the parent mass.
const Double_t mPi
Mass of pi+- (GeV/c^2)
Definition: LauConstants.hh:40
Double_t mDaugDiff_
Difference between the two daughter masses.
Double_t resRSq_
Square of the radius of the barrier for resonance decay.
const Double_t invPi
One over Pi.
Double_t FB0_
Value of the form factor for parent decay (at pole mass)
const Double_t pi
Pi.
Definition: LauConstants.hh:89
Double_t mDaugDiffSq_
Square of the difference of the two daughter masses.
Double_t p0_
Momentum of the bachelor in the resonance rest frame (at pole mass)
Double_t FR0_
Value of the form factor for resonance decay (at pole mass)
Double_t getWidth() const
Get the width of the resonance.
Double_t mDaugSum_
Sum of the two daughter masses.
Double_t resR_
Radius of the barrier for resonance decay.
Abstract class for defining type for resonance amplitude models (Breit-Wigner, Flatte etc...
void rescale(Double_t scaleVal)
Scale this by a factor.
Definition: LauComplex.hh:282
File containing LauConstants namespace.
Class for defining a complex number.
Definition: LauComplex.hh:47
Double_t resMassSq_
Square of the resonance mass.
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm)
Complex resonant amplitude.
BarrierType
Define the allowed types of barrier factors.
Double_t mDaugSumSq_
Square of the sum of the two daughter masses.
Double_t pstar0_
Momentum of the bachelor in the parent rest frame (at pole mass)
void changeResonance(Double_t newMass, Double_t newWidth, Int_t newSpin)
Allow the mass, width and spin of the resonance to be changed.
Int_t getSpin() const
Get the spin of the resonance.
Double_t d_
Extra parameter required by GS shape.