laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
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 - 2014.
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 #include <iostream>
16 
17 #include "LauConstants.hh"
18 #include "LauGounarisSakuraiRes.hh"
19 
21 
22 
23 LauGounarisSakuraiRes::LauGounarisSakuraiRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
24  LauAbsResonance(resInfo, resPairAmpInt, daughters),
25  q0_(0.0),
26  p0_(0.0),
27  pstar0_(0.0),
28  resMass_(0.0),
29  resMassSq_(0.0),
30  resWidth_(0.0),
31  resRadius_(0.0),
32  parRadius_(0.0),
33  mDaugSum_(0.0),
34  mDaugSumSq_(0.0),
35  mDaugDiff_(0.0),
36  mDaugDiffSq_(0.0),
37  mParentSq_(0.0),
38  mBachSq_(0.0),
39  h0_(0.0),
40  dhdm0_(0.0),
41  d_(0.0),
42  FR0_(1.0),
43  FP0_(1.0)
44 {
45 }
46 
48 {
49 }
50 
52 {
53  // Set-up various constants. This must be called again if the mass/width/spin
54  // of a resonance changes...
55 
56  resMass_ = this->getMass();
57  resWidth_ = this->getWidth();
58  resRadius_ = this->getResRadius();
59  parRadius_ = this->getParRadius();
60 
61  Int_t resSpin = this->getSpin();
62  Double_t massDaug1 = this->getMassDaug1();
63  Double_t massDaug2 = this->getMassDaug2();
64  Double_t massBachelor = this->getMassBachelor();
65  Double_t massParent = this->getMassParent();
66 
67  // Check that the spin is 1
68  if (resSpin != 1) {
69  std::cerr << "WARNING in LauGounarisSakuraiRes::initialise : Resonance spin is != 1. This lineshape is for the rho(770), setting the spin to 1." << std::endl;
70  this->changeResonance( -1.0, -1.0, 1 );
71  resSpin = this->getSpin();
72  }
73 
74  // Create the mass squares, sums, differences etc.
76  mDaugSum_ = massDaug1 + massDaug2;
78  mDaugDiff_ = massDaug1 - massDaug2;
80  mParentSq_ = massParent*massParent;
81  mBachSq_ = massBachelor*massBachelor;
82 
83  // Decay momentum of either daughter in the resonance rest frame
84  // when resonance mass = rest-mass value, m_0 (PDG value)
85  Double_t term1 = resMassSq_ - mDaugSumSq_;
86  Double_t term2 = resMassSq_ - mDaugDiffSq_;
87  Double_t term12 = term1*term2;
88  if (term12 > 0.0) {
89  q0_ = TMath::Sqrt(term12)/(2.0*resMass_);
90  } else {
91  q0_ = 0.0;
92  }
93 
94  // Momentum of the bachelor particle in the resonance rest frame
95  // when resonance mass = rest-mass value, m_0 (PDG value)
96  Double_t eBach = (mParentSq_ - resMassSq_ - mBachSq_)/(2.0*resMass_);
97  Double_t termBach = eBach*eBach - mBachSq_;
98  if ( eBach<0.0 || termBach<0.0 ) {
99  p0_ = 0.0;
100  } else {
101  p0_ = TMath::Sqrt( termBach );
102  }
103 
104  // Momentum of the bachelor particle in the parent rest frame
105  // when resonance mass = rest-mass value, m_0 (PDG value)
106  Double_t eStarBach = (mParentSq_ + mBachSq_ - resMassSq_)/(2.0*massParent);
107  Double_t termStarBach = eStarBach*eStarBach - mBachSq_;
108  if ( eStarBach<0.0 || termStarBach<0.0 ) {
109  pstar0_ = 0.0;
110  } else {
111  pstar0_ = TMath::Sqrt( termStarBach );
112  }
113 
114  // Calculate the Blatt-Weisskopf form factor for the case when m = m_0
115  const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
116  const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
117  FR0_ = (resBWFactor!=0) ? resBWFactor->calcFormFactor(q0_) : 1.0;
118  FP0_ = (parBWFactor!=0) ? parBWFactor->calcFormFactor(p0_) : 1.0;
119 
120  // Calculate the extra things needed by the G-S shape
121  h0_ = 2.0*LauConstants::invPi * q0_/resMass_ * TMath::Log((resMass_ + 2.0*q0_)/(2.0*LauConstants::mPi));
122  dhdm0_ = h0_ * (1.0/(8.0*q0_*q0_) - 1.0/(2.0*resMassSq_)) + 1.0/(LauConstants::twoPi*resMassSq_);
124  * TMath::Log((resMass_ + 2.0*q0_)/(2.0*LauConstants::mPi))
125  + resMass_/(LauConstants::twoPi*q0_)
127 }
128 
129 LauComplex LauGounarisSakuraiRes::resAmp(Double_t mass, Double_t spinTerm)
130 {
131  // This function returns the complex dynamical amplitude for a Breit-Wigner resonance,
132  // given the invariant mass and cos(helicity) values.
133 
134  LauComplex resAmplitude(0.0, 0.0);
135 
136  if (mass < 1e-10) {
137  std::cerr << "WARNING in LauGounarisSakuraiRes::amplitude : mass < 1e-10." << std::endl;
138  return LauComplex(0.0, 0.0);
139  } else if (q0_ < 1e-30) {
140  return LauComplex(0.0, 0.0);
141  }
142 
143  // Calculate the width of the resonance (as a function of mass)
144  // First, calculate the various form factors.
145  // NB
146  // q is the momentum of either daughter in the resonance rest-frame,
147  // p is the momentum of the bachelor in the resonance rest-frame,
148  // pstar is the momentum of the bachelor in the parent rest-frame.
149  // These quantities have been calculate in LauAbsResonance::amplitude(...)
150 
151  const Double_t resMass = this->getMass();
152  const Double_t resWidth = this->getWidth();
153  const Double_t resRadius = this->getResRadius();
154  const Double_t parRadius = this->getParRadius();
155 
156  // If the mass is floating and its value has changed we need to
157  // recalculate everything that assumes that value
158  // Similarly for the BW radii
159  if ( ( (!this->fixMass()) && resMass != resMass_ ) ||
160  ( (!this->fixResRadius()) && resRadius != resRadius_ ) ||
161  ( (!this->fixParRadius()) && parRadius != parRadius_ ) ) {
162  this->initialise();
163  }
164 
165  const Double_t q = this->getQ();
166  const Double_t p = this->getP();
167  //const Double_t pstar = this->getPstar();
168 
169  const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
170  const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
171  const Double_t fFactorR = (resBWFactor!=0) ? resBWFactor->calcFormFactor(q) : 1.0;
172  const Double_t fFactorB = (parBWFactor!=0) ? parBWFactor->calcFormFactor(p) : 1.0;
173  const Double_t fFactorRRatio = fFactorR/FR0_;
174  const Double_t fFactorBRatio = fFactorB/FP0_;
175 
176  const Double_t qRatio = q/q0_;
177  const Double_t qTerm = qRatio*qRatio*qRatio;
178 
179  const Double_t totWidth = resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
180 
181  const Double_t massSq = mass*mass;
182  const Double_t massSqTerm = resMassSq_ - massSq;
183 
184  const Double_t h = 2.0*LauConstants::invPi * q/mass * TMath::Log((mass + 2.0*q)/(2.0*LauConstants::mPi));
185  const Double_t f = resWidth * resMassSq_/(q0_*q0_*q0_) * (q*q * (h - h0_) + massSqTerm * q0_*q0_ * dhdm0_);
186 
187  // Compute the complex amplitude
188  resAmplitude = LauComplex(massSqTerm + f, resMass*totWidth);
189 
190  // Scale by the denominator factor, as well as the spin term and Blatt-Weisskopf factors
191  Double_t numerFactor = spinTerm*(1.0 + d_ * resWidth/resMass);
192  if (!this->ignoreBarrierScaling()) {
193  numerFactor *= fFactorRRatio*fFactorBRatio;
194  }
195  const Double_t denomFactor = (massSqTerm + f)*(massSqTerm + f) + resMassSq_*totWidth*totWidth;
196  resAmplitude.rescale(numerFactor/denomFactor);
197 
198  return resAmplitude;
199 }
200 
201 const std::vector<LauParameter*>& LauGounarisSakuraiRes::getFloatingParameters()
202 {
203  this->clearFloatingParameters();
204 
205  if ( ! this->fixMass() ) {
206  this->addFloatingParameter( this->getMassPar() );
207  }
208  if ( ! this->fixWidth() ) {
209  this->addFloatingParameter( this->getWidthPar() );
210  }
211  if ( ! this->fixResRadius() ) {
212  this->addFloatingParameter( this->getResBWFactor()->getRadiusParameter() );
213  }
214  if ( ! this->fixParRadius() ) {
215  this->addFloatingParameter( this->getParBWFactor()->getRadiusParameter() );
216  }
217 
218  return this->getParameters();
219 }
220 
Double_t getQ() const
Get the current value of the daughter momentum in the resonance rest frame.
LauParameter * getMassPar()
Get the mass parameter of the resonance.
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.
Bool_t fixWidth() const
Get the status of resonance width (fixed or released)
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.
void changeResonance(const Double_t newMass, const Double_t newWidth, const Int_t newSpin)
Allow the mass, width and spin of the resonance to be changed.
ClassImp(LauAbsCoeffSet)
Double_t getMassParent() const
Get the parent particle mass.
Class for defining the properties of a resonant particle.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
Double_t FP0_
Value of the form factor for parent decay (at pole mass)
Class for defininf the Gounaris-Sakurai resonance model.
Double_t resMass_
The resonance mass.
Double_t calcFormFactor(const Double_t p) const
Calculate form factor value.
Double_t h0_
Extra parameter required by GS shape.
Double_t getP() const
Get the current value of the bachelor momentum in the resonance rest frame.
File containing declaration of LauGounarisSakuraiRes class.
Double_t getMassDaug1() const
Get the mass of daughter 1.
virtual ~LauGounarisSakuraiRes()
Destructor.
Bool_t fixResRadius() const
Get the status of resonance barrier radius (fixed or released)
Double_t getParRadius() const
Get the radius of the parent barrier factor.
virtual void initialise()
Initialise the model.
Bool_t fixMass() const
Get the status of resonance mass (fixed or released)
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.
void addFloatingParameter(LauParameter *param)
Add parameter to the list of floating parameters.
LauBlattWeisskopfFactor * getParBWFactor()
Get the centrifugal barrier for the parent decay.
Double_t resRadius_
The resonance barrier radius.
const Double_t invPi
One over Pi.
std::vector< LauParameter * > & getParameters()
Access the list of floating parameters.
const Double_t pi
Pi.
Definition: LauConstants.hh:89
LauParameter * getWidthPar()
Get the width parameter of the resonance.
Double_t resWidth_
The resonance width.
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 getResRadius() const
Get the radius of the resonance barrier factor.
Double_t getWidth() const
Get the width of the resonance.
Bool_t fixParRadius() const
Get the status of parent barrier radius (fixed or released)
Double_t mDaugSum_
Sum of the two daughter masses.
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:285
Double_t parRadius_
The parent barrier radius.
File containing LauConstants namespace.
LauBlattWeisskopfFactor * getResBWFactor()
Get the centrifugal barrier for the resonance decay.
Class for defining a complex number.
Definition: LauComplex.hh:47
Class that implements the Blatt-Weisskopf barrier factor.
Double_t resMassSq_
Square of the resonance mass.
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm)
Complex resonant amplitude.
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)
Int_t getSpin() const
Get the spin of the resonance.
virtual const std::vector< LauParameter * > & getFloatingParameters()
Retrieve the resonance parameters, e.g. so that they can be loaded into a fit.
Double_t d_
Extra parameter required by GS shape.
Bool_t ignoreBarrierScaling() const
Get the ignore barrier factor scaling flag.
void clearFloatingParameters()
Clear list of floating parameters.