laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauModIndPartWave.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 //****************************************************************************
16 // Class for defining the model independent partial wave model
17 //****************************************************************************
18 
19 // --CLASS DESCRIPTION [MODEL] --
20 // Class for defining the model independent partial wave model
21 
22 #include <cstdlib>
23 #include <iostream>
24 
25 #include "TSpline.h"
26 
27 #include "LauConstants.hh"
28 #include "LauKinematics.hh"
29 #include "LauModIndPartWave.hh"
30 #include "LauResonanceInfo.hh"
31 
33 
34 
35 LauModIndPartWave::LauModIndPartWave(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters) :
36  LauAbsResonance(resInfo, resPairAmpInt, daughters),
37  nKnots_(0),
38  initialised_(kFALSE),
39  upperThresholdMag_(0.01),
40  upperThresholdPhase_(0.),
41  fixUpperThresholdMag_(kFALSE),
42  fixUpperThresholdPhase_(kFALSE)
43 {
44  this->addKnot(this->getMassDaug1() + this->getMassDaug2(),0.01,0.0);
45 }
46 
48 {
49  if(magSpline_) delete magSpline_;
50  if(phaseSpline_) delete phaseSpline_;
51 
52  std::vector<LauParameter*>::iterator it = magnitudePars_.begin();
53  std::vector<LauParameter*>::iterator end = magnitudePars_.end();
54 
55  for( ; it!=end; ++it) {
56  if(*it) delete *it;
57  }
58 
59  it = phasePars_.begin();
60  end = phasePars_.end();
61 
62  for( ; it!=end; ++it) {
63  if(*it) delete *it;
64  }
65 }
66 
68 {
69  if(!initialised_) {
70  //Use values for upper threshold knot that were set using setKnotAmp
72  initialised_ = kTRUE;
73  }
74  if(magSpline_!=0) {
75  delete magSpline_;
76  magSpline_=0;
77  }
78  if(phaseSpline_!=0) {
79  delete phaseSpline_;
80  phaseSpline_=0;
81  }
82 
83  for(Int_t i=0; i<nKnots_; ++i) {
84  magnitudes_[i] = magnitudePars_[i]->value();
85  phases_[i] = phasePars_[i]->value();
86  }
87 
88  Double_t* massVals = &masses_[0];
89  Double_t* magVals = &magnitudes_[0];
90  Double_t* phaseVals = &phases_[0];
91 
92  magSpline_ = new TSpline3("", massVals, magVals, nKnots_);
93  phaseSpline_ = new TSpline3("", massVals, phaseVals, nKnots_);
94 }
95 
96 void LauModIndPartWave::addKnot(Double_t mass, Double_t magVal, Double_t phaseVal, Bool_t fixMag, Bool_t fixPhase) {
97  const TString & parNameBase = this->getSanitisedName();
98 
99  if(mass < this->getMassDaug1() + this->getMassDaug2()) {
100  std::cerr << "WARNING in LauModIndPartWave::addKnot : Knot at mass " << mass << " is below the lower kinematic limit and will not be added." << std::endl;
101  std::cerr << " Lower kinematic limit is at mass " << this->getMassDaug1() + this->getMassDaug2() << std::endl;
102  return;
103  }
104 
105  if(mass > this->getMassParent() - this->getMassBachelor()) {
106  std::cerr << "WARNING in LauModIndPartWave::addKnot : Knot at mass " << mass << " is above the upper kinematic limit and will not be added." << std::endl;
107  std::cerr << " Upper kinematic limit is at mass " << this->getMassParent() - this->getMassBachelor() << std::endl;
108  return;
109  }
110 
111  if(!masses_.empty() && masses_[nKnots_-1] >= mass) {
112  std::cerr << "WARNING in LauModIndPartWave::addKnot : Knots must be added in ascending order. Knot at mass " << mass << " has not been added." << std::endl;
113  std::cerr << " Highest existing knot has mass " << masses_[nKnots_-1] << std::endl;
114  return;
115  }
116 
117  masses_.push_back(mass);
118  magnitudes_.push_back(magVal);
119  phases_.push_back(phaseVal);
120 
121  TString magName(parNameBase);
122  magName+="_A";
123  magName+=nKnots_;
124 
125  magnitudePars_.push_back(this->getResInfo()->getExtraParameter( magName ));
126  if( magnitudePars_[nKnots_] == 0) {
127  magnitudePars_[nKnots_] = new LauParameter( magName, magVal, 0.0, 10.0, fixMag);
128  magnitudePars_[nKnots_]->secondStage(kTRUE);
130  }
131 
132  TString phaseName(parNameBase);
133  phaseName+="_d";
134  phaseName+=nKnots_;
135 
136  phasePars_.push_back(this->getResInfo()->getExtraParameter( phaseName ));
137  if( phasePars_[nKnots_] == 0) {
138  phasePars_[nKnots_] = new LauParameter( phaseName, phaseVal, -6.0*LauConstants::pi, 6.0*LauConstants::pi, fixPhase);
139  phasePars_[nKnots_]->secondStage(kTRUE);
141  }
142 
143  std::cout << "INFO in LauModIndPartWave::addKnot : Knot added to resonance " << this->getResonanceName() << " at mass " << mass << std::endl;
144  if(fixMag) std::cout << " Magnitude fixed to " << magVal << std::endl;
145  else std::cout << " Magnitude set to " << magVal << std::endl;
146  if(fixPhase) std::cout << " Phase fixed to " << phaseVal << std::endl;
147  else std::cout << " Phase set to " << phaseVal << std::endl;
148 
149  ++nKnots_;
150 }
151 
152 void LauModIndPartWave::setKnotAmp(Int_t knot, Double_t magVal, Double_t phaseVal, Bool_t fixMag, Bool_t fixPhase) {
153 
154  //Out of range
155  if(knot > nKnots_ || knot < -1) {
156  std::cerr << "WARNING in LauModIndPartWave::setKnotAmp : Index " << knot << " does not correspond to an existing knot in resonance " << this->getResonanceName() << std::endl;
157  std::cerr << " Index must be in range -1 to " << nKnots_-1 << std::endl;
158  return;
159  }
160 
161  //Special value to access upper threshold knot (only added during initialisation)
162  if(knot == -1) {
163  upperThresholdMag_ = magVal;
164  upperThresholdPhase_ = phaseVal;
165  fixUpperThresholdMag_ = fixMag;
166  fixUpperThresholdPhase_ = fixPhase;
167 
168  std::cout << "INFO in LauModIndPartWave::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at upper kinematic threshold" << std::endl;
169  if(fixMag) std::cout << " Magnitude fixed to " << magVal << std::endl;
170  else std::cout << " Magnitude set to " << magVal << std::endl;
171  if(fixPhase) std::cout << " Phase fixed to " << phaseVal << std::endl;
172  else std::cout << " Phase set to " << phaseVal << std::endl;
173  }
174 
175  //Otherwise edit the values directly
176  else {
177  magnitudes_[knot] = magVal;
178  magnitudePars_[knot]->value(magVal);
179  magnitudePars_[knot]->genValue(magVal);
180  magnitudePars_[knot]->initValue(magVal);
181  magnitudePars_[knot]->fixed(fixMag);
182  phases_[knot] = phaseVal;
183  phasePars_[knot]->value(phaseVal);
184  phasePars_[knot]->genValue(phaseVal);
185  phasePars_[knot]->initValue(phaseVal);
186  phasePars_[knot]->fixed(fixPhase);
187 
188  if(knot == 0) std::cout << "INFO in LauModIndPartWave::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at lower kinematic threshold" << std::endl;
189  else std::cout << "INFO in LauModIndPartWave::setKnotAmp : Knot updated in resonance " << this->getResonanceName() << " at mass " << masses_[knot] << std::endl;
190  if(fixMag) std::cout << " Magnitude fixed to " << magVal << std::endl;
191  else std::cout << " Magnitude set to " << magVal << std::endl;
192  if(fixPhase) std::cout << " Phase fixed to " << phaseVal << std::endl;
193  else std::cout << " Phase set to " << phaseVal << std::endl;
194  }
195 }
196 
197 LauComplex LauModIndPartWave::resAmp(Double_t mass, Double_t spinTerm)
198 {
199  LauComplex amp(0.0, 0.0);
200 
201  Bool_t paramChanged(kFALSE);
202 
203  for(Int_t i=0; i<nKnots_; ++i) {
204  if(!magnitudePars_[i]->fixed() && magnitudePars_[i]->value() != magnitudes_[i]) {
205  paramChanged = kTRUE;
206  break;
207  }
208  if(!phasePars_[i]->fixed() && phasePars_[i]->value() != phases_[i]) {
209  paramChanged = kTRUE;
210  break;
211  }
212  }
213 
214  if(paramChanged) this->initialise();
215 
216  if (magSpline_ == 0 || phaseSpline_ == 0) {
217  std::cerr << "ERROR in LauModIndPartWave::amplitude : One of the splines is null" << std::endl;
218  return amp;
219  }
220 
221  Double_t mag = magSpline_->Eval(mass);
222  Double_t phase = phaseSpline_->Eval(mass);
223  LauComplex ff(mag*TMath::Cos(phase), mag*TMath::Sin(phase));
224 
225  amp = LauComplex(spinTerm,0.)*ff;
226 
227  return amp;
228 
229 }
230 
231 const std::vector<LauParameter*>& LauModIndPartWave::getFloatingParameters() {
232  this->clearFloatingParameters();
233 
234  for(Int_t i=0; i<nKnots_; ++i) {
235  if(!magnitudePars_[i]->fixed()) {
237  }
238  if(!phasePars_[i]->fixed()) {
239  this->addFloatingParameter( phasePars_[i] );
240  }
241  }
242 
243  return this->getParameters();
244 }
Class for defining a model independent partial wave component.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
Double_t getMassBachelor() const
Get the mass of the bachelor daughter.
const TString & getResonanceName() const
Get the name of the resonance.
std::vector< Double_t > magnitudes_
The values of the magnitude spline at the knots.
void addKnot(Double_t mass, Double_t magVal, Double_t phaseVal, Bool_t fixMag=kFALSE, Bool_t fixPhase=kFALSE)
Add a knot to the magnitude and phase splines.
File containing declaration of LauResonanceInfo class.
ClassImp(LauAbsCoeffSet)
Double_t upperThresholdPhase_
Phase for the knot at the upper kinematic threshold.
Double_t getMassParent() const
Get the parent particle mass.
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
Class for defining the properties of a resonant particle.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
std::vector< Double_t > masses_
The masses at which knots are defined in the magnitude and phase splines.
Int_t nKnots_
The number of knots in the magnitude and phase splines.
const TString & getSanitisedName() const
Get the name of the resonance.
void setKnotAmp(Int_t knot, Double_t magVal, Double_t phaseVal, Bool_t fixMag, Bool_t fixPhase)
Set the magnitude and phase at a given knot.
Double_t getMassDaug1() const
Get the mass of daughter 1.
Bool_t fixUpperThresholdPhase_
Whether the phase should be fixed for the knot at the upper kinematic threshold.
File containing declaration of LauKinematics class.
Bool_t initialised_
Flag to identify whether the model has been intialised yet.
Double_t getMassDaug2() const
Get the mass of daughter 2.
void addFloatingParameter(LauParameter *param)
Add parameter to the list of floating parameters.
virtual const std::vector< LauParameter * > & getFloatingParameters()
Retrieve the resonance parameters, e.g. so that they can be loaded into a fit.
Bool_t fixUpperThresholdMag_
Whether the magnitude should be fixed for the knot at the upper kinematic threshold.
std::vector< LauParameter * > & getParameters()
Access the list of floating parameters.
const Double_t pi
Pi.
Definition: LauConstants.hh:89
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm)
Complex resonant amplitude.
Double_t upperThresholdMag_
Magnitude for the knot at the upper kinematic threshold.
std::vector< LauParameter * > phasePars_
The parameters for the values at the knots in the phase spline.
Abstract class for defining type for resonance amplitude models (Breit-Wigner, Flatte etc...
virtual ~LauModIndPartWave()
Destructor.
LauResonanceInfo * getResInfo() const
Access the resonance info object.
File containing LauConstants namespace.
virtual void initialise()
Initialise the model.
void addExtraParameter(LauParameter *param)
Add an extra parameter of the resonance.
Class for defining a complex number.
Definition: LauComplex.hh:47
Double_t value() const
The value of the parameter.
std::vector< LauParameter * > magnitudePars_
The parameters for the values at the knots in the magnitude spline.
File containing declaration of LauModIndPartWave class.
TSpline3 * magSpline_
The spline used to calculate the magnitude at a given mass.
std::vector< Double_t > phases_
The values of the phase spline at the knots.
TSpline3 * phaseSpline_
The spline used to calculate the phase at a given mass.
void clearFloatingParameters()
Clear list of floating parameters.