laura is hosted by Hepforge, IPPP Durham
Laura++  v3r3
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauFlatteRes.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 "LauFlatteRes.hh"
19 #include "LauResonanceInfo.hh"
20 
21 #include "TSystem.h"
22 
24 
25 
26 LauFlatteRes::LauFlatteRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
27  LauAbsResonance(resInfo, resPairAmpInt, daughters),
28  g1_(0),
29  g2_(0),
30  mSumSq0_(0.0),
31  mSumSq1_(0.0),
32  mSumSq2_(0.0),
33  mSumSq3_(0.0),
34  useAdlerTerm_(kFALSE),
35  sA_(0.0)
36 {
37  Double_t g1Val(0.);
38  Double_t g2Val(0.);
39 
40  const TString& resName = this->getResonanceName();
41 
42  if ( resName == "f_0(980)" ) {
47  // constant factors from BES - Phys. Lett. B 607, 243 (2005)
48  // resMass should be 0.965 +/- 0.008 +/- 0.006 GeV/c^2
49  g1Val = 0.165; // +/- 0.010 +/- 0.015 GeV^2
50  g2Val = g1Val*4.21; // +/- 0.25 +/- 0.21
51 
52  // or from E791
53  //g1Val = 0.09;
54  //g2Val = 0.02;
55 
56  // or from CERN/WA76
57  //g1Val = 0.28;
58  //g2Val = 0.56;
59  } else if ( resName == "K*0_0(1430)" ) {
64  //Phys. Lett. B 572, 1 (2003)
65  // resMass should be 1.513 GeV/c^2
66  g1Val = 0.304;
67  g2Val = 0.380;
68  useAdlerTerm_ = kTRUE;
69  sA_ = 0.234;
70  } else if ( resName == "K*+_0(1430)" || resName == "K*-_0(1430)" ) {
75  //Phys. Lett. B 572, 1 (2003)
76  // resMass should be 1.513 GeV/c^2
77  g1Val = 0.304;
78  g2Val = 0.380;
79  useAdlerTerm_ = kTRUE;
80  sA_ = 0.234;
81  } else if ( resName == "a0_0(980)" ) {
86  //Phys. Rev. D 57, 3860 (1998)
87  // resMass should be 0.982 +/- 0.003 GeV/c^2
88  g1Val = 0.353;
89  g2Val = g1Val*1.03;
90  } else if ( resName == "a+_0(980)" || resName == "a-_0(980)" ) {
95  //Phys. Rev. D 57, 3860 (1998)
96  g1Val = 0.353;
97  g2Val = g1Val*1.03;
98  }
99 
100  const TString& parNameBase = this->getSanitisedName();
101 
102  TString g1Name(parNameBase);
103  g1Name += "_g1";
104  g1_ = resInfo->getExtraParameter( g1Name );
105  if ( g1_ == 0 ) {
106  g1_ = new LauParameter( g1Name, g1Val, 0.0, 10.0, kTRUE );
107  g1_->secondStage(kTRUE);
108  resInfo->addExtraParameter( g1_ );
109  }
110 
111  TString g2Name(parNameBase);
112  g2Name += "_g2";
113  g2_ = resInfo->getExtraParameter( g2Name );
114  if ( g2_ == 0 ) {
115  g2_ = new LauParameter( g2Name, g2Val, 0.0, 10.0, kTRUE );
116  g2_->secondStage(kTRUE);
117  resInfo->addExtraParameter( g2_ );
118  }
119 }
120 
122 {
123 }
124 
126 {
127  const TString& resName = this->getResonanceName();
128  if ( resName != "f_0(980)" && resName != "K*0_0(1430)" && resName != "K*+_0(1430)" && resName != "K*-_0(1430)" &&
129  resName != "a0_0(980)" && resName != "a+_0(980)" && resName != "a-_0(980)" ) {
130  std::cerr << "ERROR in LauFlatteRes::initialise : Unexpected resonance name \"" << resName << "\" for Flatte shape." << std::endl;
131  gSystem->Exit(EXIT_FAILURE);
132  }
133 
134  Int_t resSpin = this->getSpin();
135  if (resSpin != 0) {
136  std::cerr << "WARNING in LauFlatteRes::amplitude : Resonance spin is " << resSpin << "." << std::endl;
137  std::cerr << " : Flatte amplitude is only defined for scalers, resetting spin to 0." << std::endl;
138  this->changeResonance( -1.0, -1.0, 0 );
139  }
140 }
141 
142 LauComplex LauFlatteRes::resAmp(Double_t mass, Double_t spinTerm)
143 {
144  // This function returns the complex dynamical amplitude for a Flatte distribution
145  // given the invariant mass and cos(helicity) values.
146 
147  const Double_t resMass = this->getMass();
148  const Double_t resMassSq = resMass*resMass;
149  const Double_t s = mass*mass; // Invariant mass squared combination for the system
150 
151  const Double_t g1Val = this->getg1Parameter();
152  const Double_t g2Val = this->getg2Parameter();
153 
154  Double_t dMSq = resMassSq - s;
155  Double_t rho1(0.0), rho2(0.0);
156  if (s > mSumSq0_) {
157  rho1 = TMath::Sqrt(1.0 - mSumSq0_/s)/3.0;
158  if (s > mSumSq1_) {
159  rho1 += 2.0*TMath::Sqrt(1.0 - mSumSq1_/s)/3.0;
160  if (s > mSumSq2_) {
161  rho2 = 0.5*TMath::Sqrt(1.0 - mSumSq2_/s);
162  if (s > mSumSq3_) {
163  rho2 += 0.5*TMath::Sqrt(1.0 - mSumSq3_/s);
164  } else {
165  // Continue analytically below higher channel thresholds
166  // This contributes to the real part of the amplitude denominator
167  dMSq += g2Val*resMass*0.5*TMath::Sqrt(mSumSq3_/s - 1.0);
168  }
169  } else {
170  // Continue analytically below higher channel thresholds
171  // This contributes to the real part of the amplitude denominator
172  rho2 = 0.0;
173  dMSq += g2Val*resMass*(0.5*TMath::Sqrt(mSumSq2_/s - 1.0) + 0.5*TMath::Sqrt(mSumSq3_/s - 1.0));
174  }
175  } else {
176  // Continue analytically below higher channel thresholds
177  // This contributes to the real part of the amplitude denominator
178  dMSq += g1Val*resMass*2.0*TMath::Sqrt(mSumSq1_/s - 1.0)/3.0;
179  }
180  }
181 
182  Double_t massFactor = resMass;
183  if (useAdlerTerm_) {
184  massFactor *= ( s - sA_ ) / ( resMassSq - sA_ );
185  }
186  const Double_t width1 = g1Val*rho1*massFactor;
187  const Double_t width2 = g2Val*rho2*massFactor;
188  const Double_t widthTerm = width1 + width2;
189 
190  LauComplex resAmplitude(dMSq, widthTerm);
191 
192  const Double_t denomFactor = dMSq*dMSq + widthTerm*widthTerm;
193 
194  Double_t invDenomFactor = 0.0;
195  if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}
196 
197  resAmplitude.rescale(spinTerm*invDenomFactor);
198 
199  return resAmplitude;
200 }
201 
202 const std::vector<LauParameter*>& LauFlatteRes::getFloatingParameters()
203 {
204  this->clearFloatingParameters();
205 
206  if ( ! this->fixMass() ) {
207  this->addFloatingParameter( this->getMassPar() );
208  }
209 
210  // NB the width is given in terms of g1 and g2 so the normal width
211  // parameter should be ignored, hence it is not added to the list
212 
213  if ( ! this->fixg1Parameter() ) {
214  this->addFloatingParameter( g1_ );
215  }
216 
217  if ( ! this->fixg2Parameter() ) {
218  this->addFloatingParameter( g2_ );
219  }
220 
221  return this->getParameters();
222 }
223 
224 void LauFlatteRes::setResonanceParameter(const TString& name, const Double_t value)
225 {
226  if (name == "g1") {
227  this->setg1Parameter(value);
228  std::cout << "INFO in LauFlatteRes::setResonanceParameter : Setting g1 parameter to " << this->getg1Parameter() << std::endl;
229  } else if (name == "g2") {
230  this->setg2Parameter(value);
231  std::cout << "INFO in LauFlatteRes::setResonanceParameter : Setting g2 parameter to " << this->getg2Parameter() << std::endl;
232  } else {
233  std::cerr << "WARNING in LauFlatteRes::setResonanceParameter : Parameter name \"" << name << "\" not recognised." << std::endl;
234  }
235 }
236 
238 {
239  if (name == "g1") {
240  if ( g1_->fixed() ) {
241  g1_->fixed( kFALSE );
242  this->addFloatingParameter( g1_ );
243  } else {
244  std::cerr << "WARNING in LauFlatteRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
245  }
246  } else if (name == "g2") {
247  if ( g2_->fixed() ) {
248  g2_->fixed( kFALSE );
249  this->addFloatingParameter( g2_ );
250  } else {
251  std::cerr << "WARNING in LauFlatteRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
252  }
253  } else {
254  std::cerr << "WARNING in LauFlatteRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
255  }
256 }
257 
259 {
260  if (name == "g1") {
261  return g1_;
262  } else if (name == "g2") {
263  return g2_;
264  } else {
265  std::cerr << "WARNING in LauFlatteRes::getResonanceParameter: Parameter name not reconised." << std::endl;
266  return 0;
267  }
268 }
269 
270 void LauFlatteRes::setg1Parameter(const Double_t g1)
271 {
272  g1_->value( g1 );
273  g1_->genValue( g1 );
274  g1_->initValue( g1 );
275 }
276 
277 void LauFlatteRes::setg2Parameter(const Double_t g2)
278 {
279  g2_->value( g2 );
280  g2_->genValue( g2 );
281  g2_->initValue( g2 );
282 }
283 
virtual LauParameter * getResonanceParameter(const TString &name)
Access the given resonance parameter.
const Double_t mEta
Mass of eta (GeV/c^2)
Definition: LauConstants.hh:48
LauParameter * getMassPar()
Get the mass parameter of the resonance.
virtual void initialise()
Initialise the model.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
LauParameter * g2_
Channel 1 coupling parameter.
Double_t getMass() const
Get the mass of the resonance.
const TString & getResonanceName() const
Get the name of the resonance.
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.
File containing declaration of LauResonanceInfo class.
ClassImp(LauAbsCoeffSet)
const Double_t mK0
Mass of K0 (GeV/c^2)
Definition: LauConstants.hh:46
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
Class for defining the properties of a resonant particle.
const TString & name() const
The parameter name.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
Double_t mSumSq1_
Channel 1, subchannel 2 invariant mass.
Double_t fixg2Parameter() const
See if the g2 parameter is fixed or floating.
void setg1Parameter(const Double_t g1)
Set the g1 parameter.
Double_t mSumSq3_
Channel 2, subchannel 2 invariant mass.
void setg2Parameter(const Double_t g2)
Set the g2 parameter.
Double_t sA_
The Adler zero.
Bool_t fixMass() const
Get the status of resonance mass (fixed or released)
const Double_t mPi
Mass of pi+- (GeV/c^2)
Definition: LauConstants.hh:40
LauParameter * g1_
Channel 1 coupling parameter.
void addFloatingParameter(LauParameter *param)
Add parameter to the list of floating parameters.
Bool_t useAdlerTerm_
Flag to turn on Adler term in the width.
Double_t getg2Parameter() const
Get the g2 parameter.
virtual void floatResonanceParameter(const TString &name)
Allow the various parameters to float in the fit.
std::vector< LauParameter * > & getParameters()
Access the list of floating parameters.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:35
const Double_t mEtaPrime
Mass of eta&#39; (GeV/c^2)
Definition: LauConstants.hh:50
virtual ~LauFlatteRes()
Destructor.
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 initValue() const
The initial value of the parameter.
virtual const std::vector< LauParameter * > & getFloatingParameters()
Retrieve the resonance parameters, e.g. so that they can be loaded into a fit.
Double_t mSumSq2_
Channel 2, subchannel 1 invariant mass.
virtual void setResonanceParameter(const TString &name, const Double_t value)
Set value of a resonance parameter.
File containing LauConstants namespace.
Class for defining a complex number.
Definition: LauComplex.hh:47
Double_t mSumSq0_
Channel 1, subchannel 1 invariant mass.
Double_t fixg1Parameter() const
See if the g1 parameter is fixed or floating.
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm)
Complex resonant amplitude.
Double_t value() const
The value of the parameter.
Int_t getSpin() const
Get the spin of the resonance.
Class for defining the Flatte resonance model.
Definition: LauFlatteRes.hh:31
Double_t genValue() const
The value generated for the parameter.
const Double_t mK
Mass of K+- (GeV/c^2)
Definition: LauConstants.hh:44
const Double_t mPi0
Mass of pi0 (GeV/c^2)
Definition: LauConstants.hh:42
File containing declaration of LauFlatteRes class.
void clearFloatingParameters()
Clear list of floating parameters.
Double_t getg1Parameter() const
Get the g1 parameter.
Definition: LauFlatteRes.hh:97