laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauRelBreitWignerRes.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 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 "LauRelBreitWignerRes.hh"
19 
21 
22 
23 LauRelBreitWignerRes::LauRelBreitWignerRes(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  FR0_(1.0),
40  FP0_(1.0)
41 {
42 }
43 
45 {
46 }
47 
49 {
50  // Set-up various constants. This must be called again if the mass/width/spin
51  // of a resonance changes...
52 
53  resMass_ = this->getMass();
54  resWidth_ = this->getWidth();
55  resRadius_ = this->getResRadius();
56  parRadius_ = this->getParRadius();
57 
58  Double_t massDaug1 = this->getMassDaug1();
59  Double_t massDaug2 = this->getMassDaug2();
60  Double_t massBachelor = this->getMassBachelor();
61  Double_t massParent = this->getMassParent();
62 
63  // Create the mass squares, sums, differences etc.
65  mDaugSum_ = massDaug1 + massDaug2;
67  mDaugDiff_ = massDaug1 - massDaug2;
69  mParentSq_ = massParent*massParent;
70  mBachSq_ = massBachelor*massBachelor;
71 
72  // Create an effective resonance pole mass to protect against resonances
73  // that are below threshold
74  Double_t effResMass = resMass_;
75  Double_t effResMassSq = resMassSq_;
76  if (resMassSq_ - mDaugSumSq_ < 0.0 ){
77  Double_t minMass = mDaugSum_;
78  Double_t maxMass = massParent - massBachelor;
79  Double_t tanhTerm = std::tanh( (resMass_ - ((minMass + maxMass)/2))/(maxMass-minMass));
80  effResMass = minMass + (maxMass-minMass)*(1+tanhTerm)/2;
81  effResMassSq = effResMass*effResMass;
82  }
83 
84  // Decay momentum of either daughter in the resonance rest frame
85  // when resonance mass = rest-mass value, m_0 (PDG value)
86  Double_t term1 = effResMassSq - mDaugSumSq_;
87  Double_t term2 = effResMassSq - mDaugDiffSq_;
88  Double_t term12 = term1*term2;
89  if (term12 > 0.0) {
90  q0_ = TMath::Sqrt(term12)/(2.0*effResMass);
91  } else {
92  q0_ = 0.0;
93  }
94 
95  // Momentum of the bachelor particle in the resonance rest frame
96  // when resonance mass = rest-mass value, m_0 (PDG value)
97  Double_t eBach = (mParentSq_ - effResMassSq - mBachSq_)/(2.0*effResMass);
98  Double_t termBach = eBach*eBach - mBachSq_;
99  if ( eBach<0.0 || termBach<0.0 ) {
100  p0_ = 0.0;
101  } else {
102  p0_ = TMath::Sqrt( termBach );
103  }
104 
105  // Momentum of the bachelor particle in the parent rest frame
106  // when resonance mass = rest-mass value, m_0 (PDG value)
107  Double_t eStarBach = (mParentSq_ + mBachSq_ - effResMassSq)/(2.0*massParent);
108  Double_t termStarBach = eStarBach*eStarBach - mBachSq_;
109  if ( eStarBach<0.0 || termStarBach<0.0 ) {
110  pstar0_ = 0.0;
111  } else {
112  pstar0_ = TMath::Sqrt( termStarBach );
113  }
114 
115  // Calculate the Blatt-Weisskopf form factor for the case when m = m_0
116  const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
117  const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
118  FR0_ = (resBWFactor!=0) ? resBWFactor->calcFormFactor(q0_) : 1.0;
119  FP0_ = (parBWFactor!=0) ? parBWFactor->calcFormFactor(p0_) : 1.0;
120 }
121 
122 LauComplex LauRelBreitWignerRes::resAmp(Double_t mass, Double_t spinTerm)
123 {
124  // This function returns the complex dynamical amplitude for a Breit-Wigner resonance,
125  // given the invariant mass and cos(helicity) values.
126 
127  LauComplex resAmplitude(0.0, 0.0);
128 
129  if (mass < 1e-10) {
130  std::cerr << "WARNING in LauRelBreitWignerRes::amplitude : mass < 1e-10." << std::endl;
131  return LauComplex(0.0, 0.0);
132  } else if (q0_ < 1e-30) {
133  return LauComplex(0.0, 0.0);
134  }
135 
136  // Calculate the width of the resonance (as a function of mass)
137  // First, calculate the various form factors.
138  // NB
139  // q is the momentum of either daughter in the resonance rest-frame,
140  // p is the momentum of the bachelor in the resonance rest-frame,
141  // pstar is the momentum of the bachelor in the parent rest-frame.
142  // These quantities have been calculate in LauAbsResonance::amplitude(...)
143 
144  const Int_t resSpin = this->getSpin();
145  const Double_t resMass = this->getMass();
146  const Double_t resWidth = this->getWidth();
147  const Double_t resRadius = this->getResRadius();
148  const Double_t parRadius = this->getParRadius();
149  const Double_t q = this->getQ();
150  const Double_t p = this->getP();
151  //const Double_t pstar = this->getPstar();
152 
153  // If the mass is floating and its value has changed we need to
154  // recalculate everything that assumes that value
155  // Similarly for the BW radii
156  if ( ( (!this->fixMass()) && resMass != resMass_ ) ||
157  ( (!this->fixResRadius()) && resRadius != resRadius_ ) ||
158  ( (!this->fixParRadius()) && parRadius != parRadius_ ) ) {
159  this->initialise();
160  }
161 
162  const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
163  const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
164  const Double_t fFactorR = (resBWFactor!=0) ? resBWFactor->calcFormFactor(q) : 1.0;
165  const Double_t fFactorB = (parBWFactor!=0) ? parBWFactor->calcFormFactor(p) : 1.0;
166  const Double_t fFactorRRatio = fFactorR/FR0_;
167  const Double_t fFactorBRatio = fFactorB/FP0_;
168 
169  const Double_t qRatio = q/q0_;
170  Double_t qTerm(0.0);
171  if (resSpin == 0) {
172  qTerm = qRatio;
173  } else if (resSpin == 1) {
174  qTerm = qRatio*qRatio*qRatio;
175  } else if (resSpin == 2) {
176  qTerm = qRatio*qRatio*qRatio*qRatio*qRatio;
177  } else {
178  qTerm = TMath::Power(qRatio, 2*resSpin + 1);
179  }
180 
181  const Double_t totWidth = resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
182 
183  const Double_t massSq = mass*mass;
184  const Double_t massSqTerm = resMassSq_ - massSq;
185 
186  // Compute the complex amplitude
187  resAmplitude = LauComplex(massSqTerm, resMass*totWidth);
188 
189  // Scale by the denominator factor, as well as the spin term and Blatt-Weisskopf factors
190  resAmplitude.rescale((fFactorRRatio*fFactorBRatio*spinTerm)/(massSqTerm*massSqTerm + resMassSq_*totWidth*totWidth));
191 
192  return resAmplitude;
193 }
194 
195 const std::vector<LauParameter*>& LauRelBreitWignerRes::getFloatingParameters()
196 {
197  this->clearFloatingParameters();
198 
199  if ( ! this->fixMass() ) {
200  this->addFloatingParameter( this->getMassPar() );
201  }
202  if ( ! this->fixWidth() ) {
203  this->addFloatingParameter( this->getWidthPar() );
204  }
205  if ( ! this->fixResRadius() ) {
206  this->addFloatingParameter( this->getResBWFactor()->getRadiusParameter() );
207  }
208  if ( ! this->fixParRadius() ) {
209  this->addFloatingParameter( this->getParBWFactor()->getRadiusParameter() );
210  }
211 
212  return this->getParameters();
213 }
214 
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 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.
ClassImp(LauAbsCoeffSet)
Double_t q0_
Momentum of the daughters in the resonance rest frame (at pole mass)
Double_t getMassParent() const
Get the parent particle mass.
Double_t resMassSq_
Square of the resonance 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 resRadius_
The resonance barrier radius.
Double_t calcFormFactor(const Double_t p) const
Calculate form factor value.
virtual ~LauRelBreitWignerRes()
Destructor.
Double_t getP() const
Get the current value of the bachelor momentum in the resonance rest frame.
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm)
Complex resonant amplitude.
Double_t parRadius_
The parent barrier radius.
Double_t mDaugSum_
Sum of the two daughter masses.
Double_t resMass_
The resonance mass.
Double_t getMassDaug1() const
Get the mass of daughter 1.
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.
Double_t resWidth_
The resonance width.
File containing declaration of LauRelBreitWignerRes class.
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 FR0_
Value of the form factor for resonance decay (at pole mass)
Double_t mParentSq_
Square of the parent mass.
Double_t mBachSq_
Square of the bachelor mass.
void addFloatingParameter(LauParameter *param)
Add parameter to the list of floating parameters.
LauBlattWeisskopfFactor * getParBWFactor()
Get the centrifugal barrier for the parent decay.
Double_t p0_
Momentum of the bachelor in the resonance rest frame (at pole mass)
Double_t pstar0_
Momentum of the bachelor in the parent rest frame (at pole mass)
std::vector< LauParameter * > & getParameters()
Access the list of floating parameters.
LauParameter * getWidthPar()
Get the width parameter of the resonance.
Class for defining the relativistic Breit-Wigner resonance model.
Double_t mDaugDiff_
Difference of the two daughter masses.
Double_t mDaugDiffSq_
Square of the difference of the two daughter masses.
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)
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
File containing LauConstants namespace.
Double_t mDaugSumSq_
Square of the sum of the two daughter masses.
virtual void initialise()
Initialise the model.
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.
virtual const std::vector< LauParameter * > & getFloatingParameters()
Retrieve the resonance parameters, e.g. so that they can be loaded into a fit.
Int_t getSpin() const
Get the spin of the resonance.
void clearFloatingParameters()
Clear list of floating parameters.
Double_t FP0_
Value of the form factor for parent decay (at pole mass)