laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2
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 - 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 #include <iostream>
16 
17 #include "LauConstants.hh"
18 #include "LauRelBreitWignerRes.hh"
19 
21 
22 
23 LauRelBreitWignerRes::LauRelBreitWignerRes(TString resName, Double_t resMass, Double_t resWidth, Int_t resSpin,
24  Int_t resCharge, Int_t resPairAmpInt, const LauDaughters* daughters) :
25  LauAbsResonance(resName, resMass, resWidth, resSpin, resCharge, resPairAmpInt, daughters),
26  q0_(0.0),
27  p0_(0.0),
28  pstar0_(0.0),
29  resMassSq_(0.0),
30  mDaugSum_(0.0),
31  mDaugSumSq_(0.0),
32  mDaugDiff_(0.0),
33  mDaugDiffSq_(0.0),
34  mParentSq_(0.0),
35  mBachSq_(0.0),
36  FR0_(1.0),
37  FB0_(1.0)
38 {
39 }
40 
42 {
43 }
44 
46 {
47  // Set-up various constants. This must be called again if the mass/width/spin
48  // of a resonance changes...
49 
50  Int_t resSpin = this->getSpin();
51  Double_t resMass = this->getMass();
52  Double_t massDaug1 = this->getMassDaug1();
53  Double_t massDaug2 = this->getMassDaug2();
54  Double_t massBachelor = this->getMassBachelor();
55  Double_t massParent = this->getMassParent();
56 
57  // Create the mass squares, sums, differences etc.
58  resMassSq_ = resMass*resMass;
59  mDaugSum_ = massDaug1 + massDaug2;
61  mDaugDiff_ = massDaug1 - massDaug2;
63  mParentSq_ = massParent*massParent;
64  mBachSq_ = massBachelor*massBachelor;
65 
66  // Create an effective resonance pole mass to protect against resonances
67  // that are below threshold
68  Double_t effResMass = resMass;
69  if (resMassSq_ - mDaugSumSq_ < 0.0 ){
70  Double_t minMass = mDaugSum_;
71  Double_t maxMass = massParent - massBachelor;
72  Double_t tanhTerm = std::tanh( (resMass - ((minMass + maxMass)/2))/(maxMass-minMass));
73  effResMass = minMass + (maxMass-minMass)*(1+tanhTerm)/2;
74  }
75  Double_t effResMassSq = effResMass*effResMass;
76 
77  // Decay momentum of either daughter in the resonance rest frame
78  // when resonance mass = rest-mass value, m_0 (PDG value)
79  Double_t term1 = effResMassSq - mDaugSumSq_;
80  Double_t term2 = effResMassSq - mDaugDiffSq_;
81  Double_t term12 = term1*term2;
82  if (term12 > 0.0) {
83  q0_ = TMath::Sqrt(term12)/(2.0*effResMass);
84  } else {
85  q0_ = 0.0;
86  }
87 
88  // Momentum of the bachelor particle in the resonance rest frame
89  // when resonance mass = rest-mass value, m_0 (PDG value)
90  Double_t eBach = (mParentSq_ - effResMassSq - mBachSq_)/(2.0*effResMass);
91  Double_t termBach = eBach*eBach - mBachSq_;
92  if ( eBach<0.0 || termBach<0.0 ) {
93  p0_ = 0.0;
94  } else {
95  p0_ = TMath::Sqrt( termBach );
96  }
97 
98  // Momentum of the bachelor particle in the parent rest frame
99  // when resonance mass = rest-mass value, m_0 (PDG value)
100  Double_t eStarBach = (mParentSq_ + mBachSq_ - effResMassSq)/(2.0*massParent);
101  Double_t termStarBach = eStarBach*eStarBach - mBachSq_;
102  if ( eStarBach<0.0 || termStarBach<0.0 ) {
103  pstar0_ = 0.0;
104  } else {
105  pstar0_ = TMath::Sqrt( termStarBach );
106  }
107 
108  // Blatt-Weisskopf barrier factor constant: z = q^2*radius^2
109  // Calculate the Blatt-Weisskopf form factor for the case when m = m_0
110  const Double_t resR = this->getResBWRadius();
111  const Double_t parR = this->getParBWRadius();
112  const BarrierType barrierType = this->getBarrierType();
113  this->setBarrierRadii(resR, parR, barrierType);
114 
115  if (resSpin > 5) {
116  std::cerr << "WARNING in LauRelBreitWignerRes::initialise : Resonances spin is > 5. Blatt-Weisskopf form factors will be set to 1.0" << std::endl;
117  }
118 }
119 
120 void LauRelBreitWignerRes::setBarrierRadii(const Double_t resRadius, const Double_t parRadius, const BarrierType type)
121 {
122  // Reset the Blatt-Weisskopf barrier radius for the resonance and its parent
123  this->LauAbsResonance::setBarrierRadii( resRadius, parRadius, type );
124 
125  // Recalculate the Blatt-Weisskopf form factor for the case when m = m_0
126  const Double_t resR = this->getResBWRadius();
127  const Double_t parR = this->getParBWRadius();
128 
129  std::cout << "INFO in LauRelBreitWignerRes::setBarrierRadii : Recalculating barrier factor normalisations for new radii: resonance = " << resR << ", parent = " << parR << std::endl;
130 
131  Double_t zR0 = q0_*q0_*resR*resR;
132  Double_t zB0 = p0_*p0_*parR*parR;
133  if ( ( type == LauAbsResonance::BWPrimeBarrier ) || ( type == LauAbsResonance::ExpBarrier ) ) {
134  FR0_ = (resR==0.0) ? 1.0 : this->calcFFactor(zR0);
135  FB0_ = (parR==0.0) ? 1.0 : this->calcFFactor(zB0);
136  }
137 }
138 
140 {
141  // Calculate the requested form factor for the resonance, given the z value
142  Double_t fFactor(1.0);
143 
144  // For scalars the form factor is always unity
145  // TODO: and we currently don't have formulae for spin > 5
146  Int_t resSpin = this->getSpin();
147  if ( (resSpin == 0) || (resSpin>5) ) {
148  return fFactor;
149  }
150 
151  const BarrierType barrierType = this->getBarrierType();
152  if ( barrierType == LauAbsResonance::BWBarrier ) {
153  if (resSpin == 1) {
154  fFactor = TMath::Sqrt(2.0*z/(z + 1.0));
155  } else if (resSpin == 2) {
156  fFactor = TMath::Sqrt(13.0*z*z/(z*z + 3.0*z + 9.0));
157  } else if (resSpin == 3) {
158  fFactor = TMath::Sqrt(277.0*z*z*z/(z*z*z + 6.0*z*z + 45.0*z + 225.0));
159  } else if (resSpin == 4) {
160  fFactor = TMath::Sqrt(12746.0*z*z*z*z/(z*z*z*z + 10.0*z*z*z + 135.0*z*z + 1575.0*z + 11025.0));
161  } else if (resSpin == 5) {
162  fFactor = TMath::Sqrt(998881.0*z*z*z*z*z/(z*z*z*z*z + 15.0*z*z*z*z + 315.0*z*z*z + 6300.0*z*z + 99225.0*z + 893025.0));
163  }
164  } else if ( barrierType == LauAbsResonance::BWPrimeBarrier ) {
165  if (resSpin == 1) {
166  fFactor = TMath::Sqrt(1.0/(z + 1.0));
167  } else if (resSpin == 2) {
168  fFactor = TMath::Sqrt(1.0/(z*z + 3.0*z + 9.0));
169  } else if (resSpin == 3) {
170  fFactor = TMath::Sqrt(1.0/(z*z*z + 6.0*z*z + 45.0*z + 225.0));
171  } else if (resSpin == 4) {
172  fFactor = TMath::Sqrt(1.0/(z*z*z*z + 10.0*z*z*z + 135.0*z*z + 1575.0*z + 11025.0));
173  } else if (resSpin == 5) {
174  fFactor = TMath::Sqrt(1.0/(z*z*z*z*z + 15.0*z*z*z*z + 315.0*z*z*z + 6300.0*z*z + 99225.0*z + 893025.0));
175  }
176  } else if ( barrierType == LauAbsResonance::ExpBarrier ) {
177  if (resSpin == 1) {
178  fFactor = TMath::Exp( -TMath::Sqrt(z) );
179  } else if (resSpin == 2) {
180  fFactor = TMath::Exp( -z );
181  } else if (resSpin == 3) {
182  fFactor = TMath::Exp( -TMath::Sqrt(z*z*z) );
183  } else if (resSpin == 4) {
184  fFactor = TMath::Exp( -z*z );
185  } else if (resSpin == 5) {
186  fFactor = TMath::Exp( -TMath::Sqrt(z*z*z*z*z) );
187  }
188  }
189 
190  return fFactor;
191 
192 }
193 
194 LauComplex LauRelBreitWignerRes::resAmp(Double_t mass, Double_t spinTerm)
195 {
196  // This function returns the complex dynamical amplitude for a Breit-Wigner resonance,
197  // given the invariant mass and cos(helicity) values.
198 
199  LauComplex resAmplitude(0.0, 0.0);
200 
201  if (mass < 1e-10) {
202  std::cerr << "WARNING in LauRelBreitWignerRes::amplitude : mass < 1e-10." << std::endl;
203  return LauComplex(0.0, 0.0);
204  } else if (q0_ < 1e-30) {
205  return LauComplex(0.0, 0.0);
206  }
207 
208  // Calculate the width of the resonance (as a function of mass)
209  // First, calculate the various form factors.
210  // NB
211  // q is the momentum of either daughter in the resonance rest-frame,
212  // p is the momentum of the bachelor in the resonance rest-frame,
213  // pstar is the momentum of the bachelor in the parent rest-frame.
214  // These quantities have been calculate in LauAbsResonance::amplitude(...)
215 
216  Int_t resSpin = this->getSpin();
217  Double_t resMass = this->getMass();
218  Double_t resWidth = this->getWidth();
219  Double_t q = this->getQ();
220  Double_t p = this->getP();
221  //Double_t pstar = this->getPstar();
222 
223  const Double_t resR = this->getResBWRadius();
224  const Double_t parR = this->getParBWRadius();
225  Double_t zR = q*q*resR*resR;
226  Double_t zB = p*p*parR*parR;
227  Double_t fFactorR = (resR==0.0) ? 1.0 : this->calcFFactor(zR);
228  Double_t fFactorB = (parR==0.0) ? 1.0 : this->calcFFactor(zB);
229  Double_t fFactorRRatio = fFactorR/FR0_;
230  Double_t fFactorBRatio = fFactorB/FB0_;
231 
232  Double_t qRatio = q/q0_;
233  Double_t qTerm(0.0);
234  if (resSpin == 0) {
235  qTerm = qRatio;
236  } else if (resSpin == 1) {
237  qTerm = qRatio*qRatio*qRatio;
238  } else if (resSpin == 2) {
239  qTerm = TMath::Power(qRatio, 5.0);
240  } else {
241  qTerm = TMath::Power(qRatio, 2*resSpin + 1);
242  }
243 
244  Double_t totWidth = resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
245 
246  Double_t massSq = mass*mass;
247  Double_t massSqTerm = resMassSq_ - massSq;
248 
249  // Compute the complex amplitude
250  resAmplitude = LauComplex(massSqTerm, resMass*totWidth);
251 
252  // Scale by the denominator factor, as well as the spin term and Blatt-Weisskopf factors
253  resAmplitude.rescale((fFactorRRatio*fFactorBRatio*spinTerm)/(massSqTerm*massSqTerm + resMassSq_*totWidth*totWidth));
254 
255  return resAmplitude;
256 }
257 
Double_t getQ() const
Get the current value of the daughter momentum in the resonance rest frame.
Double_t getMassBachelor() const
Get the mass of the bachelor daughter.
Double_t getParBWRadius() const
Get the radius of the centrifugal barrier for the parent decay.
BarrierType getBarrierType() const
Get the form factor model.
Double_t getMass() const
Get the mass of the resonance.
Double_t getResBWRadius() const
Get the radius of the centrifugal barrier for the resonance decay.
Double_t FB0_
Value of the form factor for parent decay (at pole mass)
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 that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
Double_t calcFFactor(Double_t z)
Calculate the form factor for the resonance.
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 mDaugSum_
Sum of the two daughter masses.
Double_t getMassDaug1() const
Get the mass of daughter 1.
File containing declaration of LauRelBreitWignerRes class.
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.
Double_t p0_
Momentum of the bachelor in the resonance rest frame (at pole mass)
virtual void setBarrierRadii(const Double_t resRadius, const Double_t parRadius, const LauAbsResonance::BarrierType type)
Set the form factor model and parameters.
Double_t pstar0_
Momentum of the bachelor in the parent rest frame (at pole mass)
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 getWidth() const
Get the width of the resonance.
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.
Double_t mDaugSumSq_
Square of the sum of the two daughter masses.
virtual void initialise()
Initialise the model.
Class for defining a complex number.
Definition: LauComplex.hh:47
BarrierType
Define the allowed types of barrier factors.
Int_t getSpin() const
Get the spin of the resonance.
virtual void setBarrierRadii(const Double_t resRadius, const Double_t parRadius, const BarrierType type)
Set the form factor model and parameters.